UNIVERSIDADE TECNOLÓGICA FEDERAL DO PARANÁ
DEPARTAMENTO DE ENGENHARIA MECÂNICA
ENGENHARIA MECÂNICA
LUAN CAVALARI LABIGALINI
PROJETO HIDRODINÂMICO DE TURBINA HIDROCINÉTICA POR MEIO DE SIMULAÇÃO NUMÉRICA COMPUTACIONAL E TÉCNICAS
DE OTIMIZAÇÃO
TRABALHO DE CONCLUSÃO DE CURSO II
LONDRINA
2021
LUAN CAVALARI LABIGALINI
PROJETO HIDRODINÂMICODE TURBINA HIDROCINÉTICA POR MEIO DE SIMULAÇÃO NUMÉRICA COMPUTACIONAL E TÉCNICAS
DE OTIMIZAÇÃO
Trabalho de Conclusão de Curso apresentado como requisito parcial à obtenção do título de Bacharel em Engenharia Mecânica, do Departamento de Engenharia Mecânica, da Universidade Tecnológica Federal do Paraná. Orientador: Prof. Dr. Ricardo de Vasconcelos Salvo Coorientador: Prof. Dr. Rafael Sene de Lima
LONDRINA
2021
TERMO DE APROVAÇÃO
PROJETO HIDRODINÂMICO DE TURBINA HIDROCINÉTICA POR MEIO DE SIMULAÇÃO NUMÉRICA COMPUTACIONAL E TÉCNICAS DE OTIMIZAÇÃO
por
LUAN CAVALARI LABIGALINI
Este(a) Trabalho de Conclusão de Curso II foi apresentado(a) em 13 de maio de 2021
como requisito parcial para a obtenção do título de Bacharel em Engenharia
Mecânica. O(a) candidato(a) foi arguido pela Banca Examinadora composta pelos
professores abaixo assinados. Após deliberação, a Banca Examinadora considerou
o trabalho aprovado.
__________________________________ Ricardo de Vasconcelos Salvo
Prof.(a) Orientador(a)
___________________________________ Aulus Roberto Romão Bineli
Membro titular
___________________________________ Ismael de Marchi Neto
Membro titular
O Termo de Aprovação assinado encontrase na Coordenação do Curso
Ministério da Educação Universidade Tecnológica Federal do Paraná
Campus Londrina
Diretoria de Graduação e Educação Profissional Departamento Acadêmico de Engenharia Mecânica DAMECLD
Todo saber é apenas a submissão da essência da vida às leis da razão.
(TOLSTÓI, Liev, 1865)
RESUMO
LABIGALINI, Luan C. Projeto hidrodinâmico de turbina hidrocinética por meio de simulação numérica computacional e técnicas de otimização. 2021. 155. Trabalho de Conclusão de Curso (Bacharelado em Engenharia Mecânica) Universidade Tecnológica Federal do Paraná. Londrina, 2021.
Mini turbinas eólicas e turbinas hidrocinéticas são formas acessíveis de descentralizar a geração de energia a partir de recursos renováveis. A fim de melhor aproveitar essas fontes disponíveis e projetar futuros equipamentos mais eficientes, este trabalho visa desenvolver uma otimização de layout de uma turbina hidrocinética. Primeiramente, o método Blade Element Momentum (BEM) é aplicado para validar os modelos aerodinâmicos da literatura com dados experimentais. É pretendido replicar o funcionamento e o desempenho da turbina. As curvas de potência e empuxo foram reproduzidas de acordo com a variação da razão de velocidades de ponta. Depois, tal previsão fora submetida a algoritmos metaheurísticos baseados no comportamento da natureza, com intuito de encontrar os melhores parâmetros de projeto para um ambiente hipotético sujeito à demanda de eletricidade de uma única pessoa. A otimização foi realizada em relação à potência e à inércia das pás, e as variáveis de entrada foram o diâmetro do rotor, o número de pás e a velocidade de rotação. A turbina desenvolvida foi simulada computacionalmente no formato de DFC, a fim de validar os resultados de torque. Em terceiro lugar, uma correção relativa ao coeficiente de arrasto e sustentação foi conduzida com base nos efeitos de solidez entre as pás. Ele foi construído a partir de parâmetros obtidos das simulações DFC de uma série de hidrofólios. Finalmente, uma protuberância foi adicionada à borda de ataque do hidrofólio com o objetivo de diminuir o arrasto e melhorar a sustentação. O efeito da cavitação foi levado em consideração e tentou ser evitado no cálculo de dimensionamento da corda. O algoritmo mais eficiente pôde encontrar uma turbina com um coeficiente de potência 18% menor que o Limite de Betz. O efeito de solidez promoveu 8,15% de diferença em comparação com o coeficiente de potência otimizado. A protuberância prejudicou o coeficiente de sustentação e arrasto e tornou possível o efeito de cavitação. Seu impacto na inércia das pás foi um aumento de 177,61%.
Palavraschave: Turbina hidrocinética. Validação de desempenho da turbina. Curva de potência. BEM. Otimização multiobjetiva. Efeito de solidez. Protuberância em hidrofólios.
ABSTRACT
LABIGALINI, Luan C. Hydrodynamical project of a hydrokinetic turbines through computational simulation and optimization techniques. 2021. 155. Trabalho de Conclusão de Curso (Bacharelado em Engenharia Mecânica) Federal University of Technology Paraná. Londrina, 2021.
Small wind turbines (SWT) and hydrokinetic turbines (HT) are affordable ways to distribute power generation from renewable resources. In order to better harness such available sources and to design future equipment as efficient as possible, this work aims to develop a layout optimization. Firstly, the Blade Element Momentum (BEM) method is applied to validate the aerodynamic models from literature with experimental data. It intends to replicate turbine operation and its performance. The power and thrust curves could be wellpredicted according to tipspeed ratio (TSR) variation. Secondly, metaheuristic algorithms based on natural behaviour are used to find the best hydrokinetic turbine design in a hypothetical environment subjected to a single person's electricity demand. The optimization was carried out concerning power and blade inertia, and the layout parameters used as input were the rotor diameter, the number of blades, and the rotational speed. The turbine developed was simulated in a CFD manner, in order to validate torque result. Thirdly, a correction concerning drag and lift coefficient was conducted based on solidity effects between blades. It was built from parameters obtained at CFD simulations of an hydrofoil array. Finally, a protuberance was added into hydrofoil’s leading edge aiming a decrease of drag and an improvement of lift. Cavitation effect was taking into account and tried to be avoided at chord calculation. The most efficient algorithm could find a turbine with a power coefficient 18% lower than the Betz Limit. The solidity effect promoted 8,15% of difference in comparison with optimized power coefficient. The protuberance harmed lift and drag coefficient, and allowed cavitation effect. Its impact on blade inertia was a growth of 177,61%. Keywords: Hydrokinetic turbine. Turbine performance validation. Power curve. BEM. Multiobjective optimization. Solidity effect. Hydrofoil’s protuberance.
LISTA DE ILUSTRAÇÕES
Figura 1. Desenvolvimento da geração de energia elétrica per capita ao longo das últimas décadas. ....................................................................................................... 14
Figura 2. Previsão de custos para fontes renováveis e fósseis. ................................ 21
Figura 3. Representação de uma turbina Pelton com fluxo tangencial (esquerda) e, à direita de cima para baixo: turbinas com fluxo radial, misto e axial. .......................... 23
Figura 4. Tipo de turbina ideal segundo altura de coluna de água (eixo vertical) e velocidade específica (eixo horizontal) ...................................................................... 24
Figura 5. Diferentes formas construtivas de turbinas hidrocinéticas horizontais de fluxo axial. ................................................................................................................. 25
Figura 6. Razões idealizadas entre diâmetro do rotor e diâmetro do cubo, segundo velocidade específica. ............................................................................................... 26
Figura 7. Representação do efeito de ângulo de ataque às forças de arrasto e sustentação do perfil aero/hidrodinâmico. ................................................................. 28
Figura 8. Alguns dos perfis NACA quatro dígitos empregados. ................................ 29
Figura 9. Exemplo esquemático da teoria formulada por Rankine e Froude. ............ 30
Figura 10. Representação esquemática da divisão do volume de controle (a); triângulo de velocidades, a representação das forças e velocidades atuantes em cada elemento da pá (b). .......................................................................................... 32
Figura 11. Comparação entre os modelos de escoamento nas pontas das pás, para um número infinito, adotado por Betz (a), e finito considerado por Prandtl (b).......... 34
Figura 12. Representação do volume de controle considerado. ............................... 36
Figura13. Nova representação do triângulo de velocidades, considerado para o BEM. ......................................................................................................................... 36
Figura 14. Representação gráfica do parâmetro de solidez. ..................................... 43
Figura 15. Distorção da linha de corrente do escoamento, ao considerar o efeito cascata. ..................................................................................................................... 44
Figura 16. Exemplo de bolhas decorrentes de cavitação em hélices marítimas. ...... 47
Figura 17. Ilustração dos efeitos do difusor nas velocidades. ................................... 49
Figura 18. Representação do perfil de velocidade na seção transversal da pá, e consequentemente, a razão de aceleração da velocidade. ...................................... 50
Figura 19. Ilustração do fluxograma de atividades. ................................................... 62
Figura 20. Resultado da validação dos modelos para fator de indução axial. ........... 64
Figura 21. Resultado da validação dos modelos de perdas na ponta da pá, para o fator de indução axial de Glauert (1935). .................................................................. 65
Figura 22. Resultado da validação dos modelos de perdas na ponta da pá, para o fator de indução axial de Wilson e Lissaman (1974). ................................................ 66
Figura 23. Resultado da validação dos modelos de fator de indução tangencial. ..... 67
Figura 24. Resultado da validação da previsão do ângulo de ataque. ...................... 69
Figura 25. Resultado da validação de previsão do coeficiente de arrasto. ............... 70
Figura 26. Resultado da validação de previsão do coeficiente de sustentação. ....... 71
Figura 27. Resultado da validação da curva de potência. ......................................... 72
Figura 28. Resultado da validação dos valores de empuxo. ..................................... 73
Figura 29. Imagens representativas das malhas criadas e do zoom em torno do perfil NACA. ....................................................................................................................... 76
Figura 30. Validação do coeficiente de arrasto. ........................................................ 78
Figura 31. Validação do coeficiente de sustentação. ................................................ 79
Figura 32. Resultado da otimização do modelo híbrido ABCSA. ............................. 85
Figura 33. Resultado da otimização do modelo híbrido FPASA............................... 85
Figura 34. Resultado da otimização do modelo híbrido PSOSA. ............................. 86
Figura 35. Distribuição da corda ao longo do comprimento da pá. ........................... 88
Figura 36. Distribuição da torção ao longo do comprimento da pá. .......................... 89
Figura 37. Distribuição da solidez ao longo do comprimento da pá. ......................... 89
Figura 38. Comparação da distribuição da corda ao longo do comprimento da pá entre os coeficientes obtidos pelo XFOIL (DRELA, 2001)e os validados na Seção 4.1.2. ......................................................................................................................... 90
Figura 39. Comparação da distribuição da torção ao longo do comprimento da pá entre os coeficientes obtidos pelo XFOIL (DRELA, 2001) e os validados na Seção 4.1.2. ......................................................................................................................... 91
Figura 40. Representação do rotor completo segundo design otimizado. ................. 93
Figura 41. Representação esquemática do domínio de simulação DFC da turbina hidrocinética. ............................................................................................................. 94
Figura 42. Comportamento temporal do torque para cada um dos passos de tempo. .................................................................................................................................. 97
Figura 43. Comportamento temporal do torque para um passo de tempo de 8E4 (0,9°/s) ao longo de três tempos de residência. ........................................................ 98
Figura 44. Campos de velocidade instantânea por componente (da esquerda para direita): vertical, horizontal e axial. ............................................................................ 99
Figura 45. Linhas de corrente decorrentes da operação da turbina, coloridas pela magnitude da velocidade média. ............................................................................. 100
Figura 46. Visualização da esteira helicoidal através da isosuperfíce de velocidade instantânea ao longo domínio. ................................................................................ 101
Figura 47. Aproximação à ponta da pá, e visualização de seu efeito pela isosuperfíce de velocidade instantânea. ...................................................................................... 102
Figura 48. Linhas de corrente contornando a pá, coloridas pela magnitude da velocidade média. ................................................................................................... 102
Figura 49. Campos de velocidade instantânea por componente em um ponto na altura da pá (da esquerda para direita): vertical, horizontal e axial. ........................ 103
Figura 50. Campo de pressão instantânea em um ponto na altura da pá. ............. 103
Figura 51. Representação esquemática da construção do domínio numérico para simulação da solidez. .............................................................................................. 105
Figura 52. Variação do coeficiente de arrasto segundo ângulo de ataque para os valores de solidez entre 0 e 1.................................................................................. 107
Figura 53. Variação do coeficiente de sustentação segundo ângulo de ataque para os valores de solidez entre 0 e 1. ............................................................................ 108
Figura 54. Variação da razão entre o coeficiente de sustentação e arrasto segundo ângulo de ataque para os valores de solidez entre 0 e 1. ....................................... 108
Figuras 55. Campos de velocidade média para o ângulo de ataque igual à 0°, seguindo a solidez de 0,1 (a) até 1,0 (j). ................................................................. 110
Figuras 56. Campos de pressão média para o ângulo de ataque igual à 5°, seguindo a solidez de 0,1 (a) até 1,0 (j). ................................................................................. 111
Figuras 57. Campos de energia cinética turbulenta para o ângulo de ataque igual à 10°, seguindo a solidez de 0,1 (a) até 1,0 (j). .......................................................... 112
Figura 58. Variação da razão ótima entre os coeficientes de sustentação e arrasto e seus ângulos de ataque para os valores de solidez entre 0 e 1. ............................. 114
Figura 59. Comparação da distribuição de corda com e sem correção de solidez . 115
Figura 60. Comparação da distribuição de torção com e sem correção de solidez.115
Figura 61. Resultado da otimização do modelo híbrido ABCSA corrigido segundo resultados de solidez. .............................................................................................. 116
Figura 62. Distribuição da corda ao longo do comprimento da pá, corrigido pelos efeitos da solidez. .................................................................................................... 118
Figura 63. Distribuição da torção ao longo do comprimento da pá, corrigido pelos efeitos da solidez. .................................................................................................... 118
Figura 64. Distribuição da solidez ao longo do comprimento da pá, corrigido pelos seus efeitos no escoamento. ................................................................................... 119
Figura 65. Comparação do coeficiente de arrasto do perfil com e sem ressalto. .... 121
Figura 66. Comparação do coeficiente de sustentação do perfil com e sem ressalto. ................................................................................................................................ 121
Figuras 67. Campos de energia cinética turbulenta para o ângulo de ataque igual à 10°(a) sem e (b) com o ressalto. ............................................................................. 123
Figuras 68. Campos de velocidade média para o ângulo de ataque igual à 10°(a) sem e (b) com o ressalto. ........................................................................................ 124
Figuras 69. Campos de energia cinética turbulenta para o ângulo de ataque igual à 18°(a) sem e (b) com o ressalto. ............................................................................. 125
Figuras 70. Campos de pressão média para o ângulo de ataque igual à 18°(a) sem e (b) com o ressalto. ................................................................................................... 126
Figuras 71. Campos de velocidade média para o ângulo de ataque igual à 18°(a) sem e (b) com o ressalto. ........................................................................................ 127
Figura 72. Variação do coeficiente de arrasto para o perfil com ressalto, segundo ângulo de ataque para os valores de solidez entre 0 e 1. ....................................... 128
Figura 73. Variação do coeficiente de sustentação para o perfil com ressalto, segundo ângulo de ataque para os valores de solidez entre 0 e 1. ........................ 129
Figura 74. Variação da razão entre o coeficiente de sustentação e arrasto para o perfil com ressalto, segundo ângulo de ataque para os valores de solidez entre 0 e 1. ................................................................................................................................ 129
Figuras 75. Campos de energia cinética turbulenta para o ângulo de ataque igual à 10° e uma solidez de 1,0 de perfis (a) sem e (b) com o ressalto............................. 131
Figuras 76. Campos de velocidade média para o ângulo de ataque igual à 10° e uma solidez de 1,0 de perfis (a) sem e (b) com o ressalto. ............................................. 132
Figura 77. Variação da razão ótima entre os coeficientes de sustentação e arrasto para o perfil com ressalto e seus ângulos de ataque, para os valores de solidez entre 0 e 1. ....................................................................................................................... 133
Figura 78. Comparação da distribuição de corda com e sem ressalto. ................... 134
Figura 79. Comparação da distribuição de torção com e ressalto. .......................... 134
Figura 80. Resultado da otimização do modelo híbrido ABCSA para um perfil com ressalto, e também corrigido segundo resultados de solidez. ................................. 135
Figura 81. Distribuição da corda ao longo do comprimento da pá, segundo coeficientes de perfil com ressalto e corrigido pelos efeitos da solidez. .................. 137
Figura 82. Distribuição da torção ao longo do comprimento da pá, segundo coeficientes de perfil com ressalto e corrigido pelos efeitos da solidez. .................. 137
Figura 83. Distribuição da solidez ao longo do comprimento da pá, segundo coeficientes de perfil com ressalto e corrigido pelos efeitos da solidez. .................. 138
Figura 84. Comparação da distribuição da corda ao longo do comprimento da pá com e sem correção que evita a cavitação, segundo coeficientes de perfil com ressalto e corrigido pelos efeitos da solidez. ........................................................... 139
LISTA DE TABELAS
Tabela 1 Validação dos modelos de coeficiente de potência no ponto experimental. .................................................................................................................................. 68
Tabela 2 Resumo dos modelos validados em replicar funcionamento da turbina. .. 74
Tabela 3 Detalhes das malhas construídas para a validação do perfil isolado. ...... 75
Tabela 4 Condições iniciais para as simulações de perfil isolado. .......................... 77
Tabela 5 Constantes de projeto e operação. .......................................................... 81
Tabela 6 Delimitação da região de domínio. ........................................................... 82
Tabela 7 Constantes dos modelos de otimização. .................................................. 84
Tabela 8 Resultado dos designs ótimos para turbina. ............................................ 86
Tabela 9 Resultado da geometria otimizada para turbina, através do modelo ABCSA segundo os coeficientes validados na Seção 4.1.2. ............................................ 92
Tabela 10 Resultado da geometria otimizada para turbina, através do modelo ABCSA segundo diferentes números de pessoas como demanda. ................................. 92
Tabela 11 Detalhes das malhas construídas para a simulação da turbina. ............. 95
Tabela 12 Condições iniciais para as simulações da turbina. ................................. 95
Tabela 13 Resultado das simulações de DFC para independência de passo de tempo. ....................................................................................................................... 97
Tabela 14 Quantidade de elementos para cada uma das malhas de solidez. ...... 105
Tabela 15 Condições iniciais para as simulações de solidez. ............................... 106
Tabela 16 Resultado dos coeficientes de pressão mínimo segundo ângulo ataque para os valores de solidez entre 0 e 1. .................................................................... 109
Tabela 17 Comparação dos resultados de design ótimos para turbina com e sem correção de solidez. ................................................................................................ 117
Tabela 18 Resultado dos coeficientes de pressão mínimo para o perfil com ressalto, segundo ângulo ataque para os valores de solidez entre 0 e 1. ............... 130
Tabela 19 Comparação dos resultados de design ótimos para turbina de perfil: com e sem ressalto, e também com e sem correção de solidez. .................................... 136
LISTA DE SIGLAS, ACRÔNIMOS E SÍMBOLOS
LISTA DE SIGLAS
ABC Artificial Bee Colony DFC Dinâmica dos Fluidos Computacional FPA Flower Pollination Algorithm GCI Grid Convergence Index IBGE Instituto Brasileiro de Geografia e Estatísticas MRF Moving Reference Frame ONS Operador Nacional do Sistema elétrico OWD Our World in Data PSO Particle Swarm Optimization SA Simulated Annealing SDG Sustainable Development Goals TTR Teorema do Transporte de Reynolds UNDP United Nations Development Programme
LISTA DE ACRÔNIMOS
ANA Agência Nacional das Águas BEM Blade Element Momentum BET Blade Element Theory EPE Empresa de Pesquisa Energética GREENS Ground Renewable Expeditionary Energy Systems IEA International Energy Association NACA U.S. National Advisory Committee on Aeronautics NASA National Aeronautics and Space Administration ONG Organização Não Governamental ONU Organização das Nações Unidas RANS ReynoldsAveraged NavierStokes
LISTA DE SÍMBOLOS
a Fator de indução axial no rotor/disco a’ Fator de indução tangencial A Área [m²]
b Fator de indução axial na esteira B Número de pás c Corda do perfil [m] CD Coeficiente de Arrasto do perfil CL Coeficiente de Sustentação do perfil Cn Coeficiente de Força Normal CP Coeficiente de Potência Cp Coeficiente de Pressão Ct Coeficiente de Força Tangencial CT Coeficiente de Empuxo D Diâmetro do rotor [m] DHub Diâmetro do hub [m] F Fator de correção nas pontas fs Fator de segurança g Gravidade [m/s²] h Profundidade local [m] H Altura de coluna d’água [m] J Inércia total [kg*m2] KA Coeficiente de área L Comprimento característico [m] m Massa [kg] Nspec Velocidade específica da turbina P Potência [W] p Pressão [kPa] p0 Pressão de referência [kPa] patm Pressão atmosférica [kPa] pv Pressão de vapor do líquido [kPa] q Ordem de convergência Q Torque [N*m]
r Posição radial [m] R Raio do rotor [m] RHub Raio do hub [m] T Empuxo [Nm] t Máxima espessura do perfil hidrodinâmico U0 Velocidade de corrente livre [m/s] ui Componente do vetor velocidade ui Média temporal da velocidade ui’ Flutuação temporal da velocidade VCAV Velocidade de cavitação W Velocidade relativa [m/s] α Ângulo de Ataque [°] θ Ângulo de Torção da Pá [°] λ Razão entre a velocidade tangencial da ponta da pá com a velocidade
de corrente livre λr Razão local da velocidade tangencial da pá com a velocidade de
corrente livre μ Viscosidade dinâmica [kg/m*s]
μt Viscosidade dinâmica turbulenta [kg/m*s]
ρ Densidade absoluta do fluido [kg/m3]
ρm Densidade absoluta do material [kg/m3]
σ Solidez (ou Solidity) σ CAV Índice de cavitação Φ Ângulo de Escoamento [°] Ω Velocidade angular do fluido no plano do disco [rad/s] ε Coeficiente de peso γ Coeficiente parabólico da demanda de potência
ω Velocidade angular do disco/rotor [rad/s] δ Erro relativo dos parâmetros em análise de malhas τ Razão entre malhas φ Coeficiente em análise de malhas
SUMÁRIO
1 INTRODUÇÃO .....................................................................................................13
1.1 OBJETIVO ........................................................................................................18
1.2 JUSTIFICATIVA ................................................................................................19
2 REVISÃO BIBLIOGRÁFICA ................................................................................22
2.1 TURBINAS HIDRÁULICAS E HIDROCINÉTICAS ............................................22
2.2 SIMULAÇÃO NUMÉRICA COMPUTACIONAL .................................................30
2.2.1 Principais desenvolvedores de teorias relacionadas aos rotores ...................30
2.2.1.1 RankineFroude ..........................................................................................30
2.2.1.2 FroudeDrzewiecki (Teoria do Elemento da Pá) .........................................32
2.2.1.3 Prandtl e Betz..............................................................................................33
2.2.1.4 Glauert (Teoria do Momento da Pá) ............................................................35
2.2.1.5 Wilson e Lissaman ......................................................................................37
2.2.1.6 Spera ..........................................................................................................38
2.2.1.7 Shen ............................................................................................................39
2.2.1.8 Demais correções .......................................................................................40
2.2.1.8.1 Fatores de indução axial e tangencial .......................................................40
2.2.1.8.2 Fator de perda nas pontas das pás ..........................................................41
2.2.1.8.3 Coeficientes de empuxo e potência ..........................................................41
2.2.1.9 Efeito cascata das pás ................................................................................42
2.2.1.10 Cavitação ...................................................................................................45
2.2.1.11 Difusores ...................................................................................................48
2.2.2 Dinâmica dos fluidos computacional ...............................................................51
2.2.2.1 Modelos de turbulência ...............................................................................52
2.2.2.2 Métodos numéricos .....................................................................................54
2.3 OTIMIZAÇÃO ....................................................................................................55
2.3.1 Particle Swarm Optimization ...........................................................................56
2.3.2 Artificial Bee Colony ........................................................................................57
2.3.3 Flower Pollination Algorithm............................................................................58
2.3.4 Simulated Annealing .......................................................................................59
3 METODOLOGIA ...................................................................................................61
4 RESULTADOS .....................................................................................................63
4.1 VALIDAÇÃO .....................................................................................................63
4.1.1 MODELOS TEÓRICOS ..................................................................................63
4.1.2 PARÂMETROS DE SIMULAÇÃO DFC ..........................................................74
4.2 OTIMIZAÇÃO ....................................................................................................80
4.2.1 DESIGN ..........................................................................................................80
4.2.1.1 Algoritmo .....................................................................................................80
4.2.1.2 Simulação DFC da turbina hidrocinética .....................................................93
4.2.2 SOLIDEZ ........................................................................................................104
4.2.2.1 Simulações DFC dos perfis em série ..........................................................104
4.2.2.2 Algoritmo .....................................................................................................113
4.2.3 RESSALTO .....................................................................................................120
4.2.3.1 Simulações DFC dos perfis com ressalto ...................................................120
4.2.3.2 Algoritmo .....................................................................................................132
5 CONCLUSÃO .......................................................................................................140
REFERÊNCIAS .......................................................................................................142
13
1 INTRODUÇÃO
A luta contra as mudanças climáticas, agravada pelo crescimento da
população mundial, é uma das maiores preocupações da humanidade no século 21.
Em 2011, as 7,01 bilhões de pessoas ao redor do mundo (ROSER; RITCHIE; ORTIZ
OSPINA, 2019) foram abastecidas por 550 EJ de energia (SIIROLA, 2014). A título de
comparação, em 2016, o suprimento de energia cresceu 4,73% (IEA, 2018), enquanto
que no mesmo período a humanidade já era composta por 7,43 bilhões de pessoas –
6% a mais em relação a 2011.
Significa dizer que, atualmente, a taxa de crescimento da população mundial
é maior que a taxa de crescimento do fornecimento da energia. Se considerarmos
apenas a eletricidade, os números de geração de energia elétrica percapita seguem
a mesma tendência – por mais que apresente taxa relativamente maior. Como pode
ser visto na Figura 1, ao considerar o mesmo período de 2011 à 2016, em escala
global, esse valor médio foi de aproximadamente 5%; ao passo que países em
desenvolvimento como a Índia, China e Brasil os números são muito mais expressivos:
o primeiro teve quase 28% de aumento, o segundo mais de 25% e o terceiro em torno
de 4% (ROSER; RITCHIE; ORTIZOSPINA, 2019; RITCHIE, 2021). Além disso, um
fator crítico é o fato que esse suprimento não fora evoluído de maneira “limpa”, ou
seja, livre de emissão de gases poluentes, uma vez que os níveis de dióxido de
carbono apenas se intensificaram (IEA, 2018).
14
Figura 1. Desenvolvimento da geração de energia elétrica per capita ao longo das últimas décadas.
Fonte: Próprio autor. Dados de Ritchie (2021) e Roser; Ritchie; OrtizOspina (2019).
Por mais que o carbono siga seu ciclo natural, ele está desbalanceado
(SIIROLA, 2014). A taxa de emissão de carbono está excedendo a razão de
dissolução deste nos oceanos, e quanto mais moléculas de dióxido de carbono
existirem na atmosfera, piores as consequências para o clima mundial. Por isso,
políticas de estado, regulamentações e tratados são feitos não apenas com intuito de
atender uma crescente demanda de energia global de maneira sustentável, mas
também para reduzir as emissões dos gases de origem fóssil (BRACKEN;
BULKELEY; MAYNARD, 2014).
Em 1972, sob o contexto histórico da conturbada década de 60 e
curiosamente sob a comoção mundial da primeira foto registrada da Terra, foi
convocada pela Organização das Nações Unidas (ONU) (“Meio ambiente”, 2014) a
primeira Conferência das Nações Unidas sobre o Ambiente Humano, na cidade de
Estocolmo, na Suécia. Neste evento, foram estabelecidos 19 princípios no “Manifesto
Ambiental”, buscando inspirar as nações quanto a preservação do meio ambiente.
Passados vinte anos, motivados principalmente pela descoberta do buraco da
camada de ozônio e do crescente desmatamento na floresta Amazônica (CARVALHO,
2012), mais de 100 chefes de estado se encontraram no evento conhecido como ECO
15
92 ou “Cúpula da Terra” – na cidade do Rio de Janeiro. Nele, os representantes se
comprometeram em adotar a chamada “Agenda 21”, que incentiva o desenvolvimento
econômico sustentável tanto no âmbito de proteção ambiental, quanto na redução de
pobreza extrema.
Entretanto, o compromisso fora ratificado apenas cinco anos depois, em
Kyoto, no Japão. No conhecido “Protocolo de Kyoto”, delimitouse legalmente que as
nações reduziriam as emissões de gases do efeito estufa, com metas específicas para
cada país – essas mais rigorosas para os países desenvolvidos, considerados
responsáveis pelo alto nível de poluição dos últimos 150 anos (QUADROS, 2017).
Já em 2015, após inúmeras constatações das mudanças climáticas, a
comunidade científica global alertou no sentido de que, nos anos que viriam, o
aumento da temperatura média global deveria ser limitado em 2ºC. Assim sendo, foi
assinado em Paris, na França, o mais recente acordo climático, denominado “Acordo
de Paris”. Dessa vez, todos os países que ratificassem o referido acordo deveriam
reduzir os níveis de emissão dos gases de efeito estufa.
Além disso, em Nova York, os mesmos países se comprometeram a seguir a
chamada “Agenda 2030”, que inclui os Objetivos para o Desenvolvimento Sustentável
(ODS, ou como mais usual, em inglês Sustainable Development Goals SDGs): como
o sétimo objetivo da Agenda, está a “Energia Limpa e Acessível”, cujos principais
tópicos são: 7.1 em relação ao acesso à energia; 7.2 pelo aumento de
compartilhamento e consumo de energia provida de fontes renováveis e 7.3 que
objetiva o desenvolvimento de melhoras relacionadas a eficiência energética global.
O SDG7 aponta que quase 13% da população mundial (aproximadamente 1
bilhão de pessoas) vivem sem acesso à energia elétrica, sendo que a maioria vive em
países subdesenvolvidos na África Subsaariana e central e na América Central e do
Sul. É importante salientar que mesmo os países em que o acesso universal à energia
elétrica foi atingido entre 2010 e 2016 (SDG7), como por exemplo o Brasil, ainda
possuem extensas áreas cuja população não tem tal acesso – isso porque o território
nacional de 8.515 milhões de m² (IBGE, 2012) e os locais de difícil acesso tornam a
transmissão de energia elétrica um desafio. A título de exemplo, ainda em 2018 o
estado de Roraima, no Brasil, é o único não interligado ao SIN (Sistema Interligado
Nacional), a rede elétrica interligada brasileira (ELETROBRAS, 2018), sendo
necessário importar a energia elétrica provida pela Venezuela; infelizmente, em 2019,
a Venezuela cortou o fornecimento de energia ao estado de Roraima, sendo
16
necessária a reativação de diversas termoelétricas para o abastecimento da
população local (OLIVEIRA, 2020). A previsão é que essa situação não se altere ao
menos até o ano de 2024 (ONS).
É por isso que muitas comunidades locais apelam para uma produção
regional, em que na maioria das vezes se faz o uso de motogeradores a Diesel.
Nesse sentido, o tópico 7.2 do SDG7 explicitamente demonstra quão delicada é a
questão da mudança climática e a emissão de gases ligados ao efeito estufa.
Dessa forma, os recursos de energia renovável se apresentam como uma
solução viável e têm se consolidado à matriz energética mundial, uma vez que as
percentagens de produção a partir de biomassa, água, solar e vento aumentaram.
Assim, mesmo que representem quase 20% do fornecimento em 2016, é notável o
crescimento desse percentual em relação aos 5 anos anteriores: quase 420% para as
placas fotovoltaicas (63.170 para 328.038 GWh), aproximadamente 120%
proveniente do vento (436.010 para 957.694 GWh) e 16% para a mais estabelecida e
com números expressivos, energia provida pela água (3.597.861 para 4.170.035
GWh), na qual a principal produtora é a China (IEA, 2018).
Em 2012 essa última fonte era essencial em mais de 150 países, sendo que
em cerca de 25 destes houve dependência para suprimir 90% de sua demanda
elétrica, e destes, ainda, 12 eram 100% dependentes (IRENA, 2012). No cenário
nacional, em 2019, esse valor ultrapassava os 63%, em que 60,19% provinham de
usinas hidrelétricas e 3,59% de pequenas e micro centrais hidrelétricas (ANEEL,
2019). A energia provida pela água é amplamente utilizada na forma de hidroelétricas,
cujas barragens permitem a sua utilização na forma de energia potencial, e assim,
produzem uma grande quantidade de energia elétrica. Entretanto, seu dano ambiental
e social pode ser incalculável por conta das largas áreas alagadas (KAUNDA;
KIMAMBO; NIELSEN, 2014). Do ponto de vista ecológico, por exemplo, há
impedimento do movimento natural dos peixes, alteração da qualidade da água (em
termos de temperatura, dissolução de gases e sedimentação) e alteração de habitat
aquático natural (LEHR; KEELEY; KINGERY, 2016).
Entretanto, atualmente a estratégia de utilização da água como fonte de
energia não é restrita e tem se desenvolvido com base no uso de sua energia cinética,
abundante por exemplo em correntes marítimas e em rios. Além disso, a rede elétrica
não possui infraestrutura suficiente para se estender por todo o território brasileiro, e
consequentemente os locais isolados são prejudicados. Tendo isso em vista, turbinas
17
hidrocinéticas têm sido foco de pesquisas, uma vez que as baixas velocidades de
corrente e a alta densidade de energia tornam mais acessível o atendimento à
demanda local de algumas famílias em áreas remotas (VERMAAK; KUSAKANA;
KOKO, 2014).
Para facilitar o desenvolvimento deste tipo de projeto, um estudo mais
aprofundado é necessário e isso pode ser traduzido em prototipagem. Mas antes de
qualquer coisa, uma maneira menos custosa e mais otimizada para isso é uma análise
preliminar através de simulações numéricas, seja por algoritmos de performance, seja
pela Dinâmica dos Fluidos Computacional (DFC, traduzido do inglês Computational
Fluid Dynamics, CFD). A primeira trata de uma análise genérica, em que o resultado
informa parâmetros tangíveis e importantes para um dimensionamento prévio. Para
uma investigação mais minuciosa, utilizase a segunda.
A DFC tem como base a construção de uma cópia computacional do modelo
físico, e após aplicar a técnica de discretização às equações governantes, a qual
transforma as equações diferenciais governantes em relações algébricas discretas.
Estas são resolvidas por meio de métodos numéricos (LEHR; KEELEY; KINGERY,
2016). Como tratado por Ardizzon; Cavazzini; Pavesi (2014), a DFC é uma ferramenta
útil para avaliações precisas sobre a previsão da operação de turbomáquinas, se
tornando o estado da arte para esse tipo de projeto.
Ambos tratamentos são importantes meios para fornecer informações para
tomadas de decisões no desenvolvimento de projetos, ou mesmo para proporcionar
dados para que essa tarefa seja cada vez mais precisa. Aliado à evolução de
componentes eletrônicos, linguagens de programação e processadores
computacionais, a investigação dos problemas se tornou mais sofisticada. Como
consequência, diversas técnicas puderam ser desenvolvidas, como o aprendizado de
máquinas e a inteligência artificial, que têm como base a mimetização da natureza
procurando fazer com que as máquinas tenham “inteligência suficiente para tomar
decisões de um assunto específico” aliando a velocidade cada vez maior dos
processadores computacionais com a complexa ideia da formação do conhecimento
humano.
Nesse sentido, citase também os modelos de otimização, cujo objetivo é
maximizar ou minimizar uma equação sob um domínio de soluções factíveis ao
problema, segundo uma série de critérios (SIERKSMA, 2002). E, dessa mesma forma,
alguns dos algoritmos mais utilizados foram construídos com base na teoria genética
18
(algoritmos genéticos), movimento dos pássaros em busca de alimento (PSO) e
colônias de abelhas (ABC).
Um grande exemplo da aplicação que traduz a otimização de processos,
aprendizado de máquinas e DFC, são os “Digital Twins”. Como o próprio nome diz,
são cópias digitais dos modelos reais, cujo funcionamento é monitorado por sensores.
Esses sensores providenciam as informações necessárias para a previsão da
operação no modelo virtual pela DFC. Desta antevisão, os algoritmos são capazes de
otimizar o desempenho, e assim, “a máquina” tomar a melhor decisão segundo uma
série de condições.
A partir desses cenários é que o presente trabalho se desenvolve. No primeiro
capítulo, nas seções seguintes, são apresentados os objetivos e a motivação para o
estudo. Já o segundo capítulo, trata de todo conhecimento teórico e específico
necessário para o desenvolvimento dos objetivos. No terceiro é apresentada a
metodologia de desenvolvimento em que os objetivos são abordados. No quarto
capítulo os principais resultados são apresentados, assim como sua análise. O quinto,
conta com a conclusão do trabalho desenvolvido.
1.1 OBJETIVO
Sob as perspectivas apresentadas, têmse como objetivo principal o
desenvolvimento do projeto hidrodinâmico de uma turbina hidrocinética, de forma que
sua performance e seu design sejam otimizados a partir de dados obtidos por DFC e
por modelos analíticos desenvolvidos ao longo da história.
Devese salientar, que a densidade energética por pessoa que essa deve
gerar seja o suficiente segundo previsões nacionais no âmbito residencial cujo valor
representa pouco mais de 250 W. Ainda, assumese uma produção isolada, ou seja,
que a potência de projeto será baseada na média diária e que para suprir os picos de
demanda necessite a utilização de uma bateria.
19
1.2 JUSTIFICATIVA
O contexto apresentado na primeira seção quanto às energias renováveis, em
especial no estudo de fontes hídricas e do desenvolvimento tecnológico, evidencia um
campo com vastas oportunidades a serem exploradas e crescimento da performance
tangível. Ao ampliar a possibilidade de aproveitar os recursos hídricos para uma
geração localizada, o potencial hidráulico total dos países é expandido: 2.474 TWh da
China, 1.670 TWh da Rússia, 1.339 TWh dos EUA e 1.250 TWh do Brasil, anualmente,
por exemplo. Ao todo, a capacidade mundial é de aproximadamente 15.955 TWh/ano,
ou 4,8 vezes a mais que a produção dessa fonte em relação a 2012 (IRENA, 2012).
Entretanto, dos países citados o Brasil é o que lidera no quesito de
aproveitamento de tais recursos: 84% do disponível já é aproveitado; enquanto os
outros ainda possuem larga margem de aproveitamento: apenas 19% dos recursos
russos e 16% dos chineses são usufruídos (IRENA, 2012). A política de globalizar o
alcance da energia elétrica, tratado no SDG7, também é diretamente beneficiada: das
hoje aproximadamente 1 bilhão de pessoas sem acesso à eletricidade (HART;
SALING, 2012), 1/3 por outro lado têm acesso à água em movimento (VERMAAK;
KUSAKANA; KOKO, 2014).
Se tratando especificamente do quadro hidrocinético – que pode incluir as
correntes fluviais e/ou as marítimas , o potencial possível a ser explorado, por
exemplo, nos Estados Unidos é de 120 TWh/ano (JACOBSON, 2012), enquanto que
na Inglaterra é 22 TWh/ano (BLACK & VEATCH, 2005), e na Irlanda, 10,46 TWh/ano
(O’ROURKE; BOYLE; REYNOLDS, 2010). Estudos desse tipo, de alcance nacional,
até o momento, não foram realizados no Brasil. Entretanto, alguns potenciais em
específico foram avaliados: por meio de dados geográficos e com base no
funcionamento de turbinas já existentes, o potencial anual de potência localizada na
jusante do rio Iguaçu, cuja velocidade do rio pode variar entre 1,7 e 2,5 m/s, é de
159,89 MWh, enquanto que no rio Paraná esse valor é de 680,26 MWh – sujeito a
uma velocidade que pode chegar até quase 1,6 m/s (ARAÚJO, 2016); por meio de
modelagem numérica computacional e de aproveitamento viável do arranjo de
turbinas: 815,3 MWh/ano para o rio Jamari localizado no estado de Rondônia, cuja
velocidade de corrente gira em torno de 1,5 m/s, e 258,1 MWh/ano para o rio Curuá
Una, no estado do Pará, que atinge uma velocidade média de correntezas superior a
20
1,2 m/s (SANTOS, 2019). Isso sem contar os outros rios que juntos compõem os 12%
da totalidade de água livre do mundo, segundo a Agência Nacional das Águas (ANA),
e as correntes marítimas, que não foram até então consideradas.
Diversos esforços têm sido realizados com intuito de promover esse tipo de
tecnologia: seja militar, pelo projeto GREENS (SCHLEICHER; RIGLIN; OZTEKIN,
2015); seja por meio de Organizações Não Governamentais (ONGs), pela chamada
Practical Action; seja por programas vinculados à ONU, pelo Programa de
Desenvolvimento das Nações Unidas (PDNU, ou do inglês United Nations
Development Programme – UNDP) (KAUNDA; KIMAMBO; NIELSEN, 2014). Além
disso, também têm na versatilidade um destaque, podendo ser utilizado em uma
geração híbrida. Nesta configuração, diversas fontes renováveis são usufruídas em
conjunto, e isso traz algumas vantagens diretas, como a descentralização da rede e
sua invulnerabilidade às oscilações de demanda energética (ARDIZZON; CAVAZZINI;
PAVESI, 2014).
Diferentes classificações podem ser dadas quanto a escala da potência
elétrica gerada, variando de acordo com o autor e região. Micro unidades geradoras
podem ser classificadas entre 0,5 e 100 kW (SCHLEICHER; RIGLIN; OZTEKIN,
2015), as mini unidades entre 100 e 1.000 kW, e as pequenas, de 1.000 a 30.000 kW
(VERMAAK; KUSAKANA; KOKO, 2014). Isso afeta diretamente os custos de
instalação: por conta da grande quantidade de energia gerada, as hidrelétricas tendem
a ter uma densidade de custo por energia menor, diferente do apresentado pelas
unidades de pequena escala (IRENA, 2012). Uma comparação entre os gastos
previstos da geração de energia elétrica pode ser vista na Figura 2, em que nela se
vê quão vantajoso a produção por fonte hídrica é perante às demais fontes renováveis
e fósseis.
21
Figura 2. Previsão de custos para fontes renováveis e fósseis.
Fonte: WILLIAMS; SIMPSON, (2009).
Por outro lado, as consequências ambientais e sociais são incalculáveis,
como já citado; e quando comparado às outras fontes renováveis, por mais que exija
um investimento inicial maior (algo em torno de 3125 U$/kW, contra 2000 U$/kW da
eólica, por exemplo (CORRÊA DA SILVA; DE MARCHI NETO; SILVA SEIFERT,
2016)), seu tempo de amortização de investimento (payback) pode ser facilmente
atingido, uma vez que a vida útil das unidades micro geradoras, por exemplo, gira em
torno dos 50 anos (KAUNDA; KIMAMBO; NIELSEN, 2014). Considerando a geração
hidrocinética, ressaltase ainda que o projeto pode ser feito com certa flexibilidade, já
que a cada aumento da velocidade de corrente livre da água – onde sua densidade
por si só já possibilita uma geração elétrica maior –, há um acréscimo cúbico de
potência gerada (VERMAAK; KUSAKANA; KOKO, 2014).
22
2 REVISÃO BIBLIOGRÁFICA
Este tópico é dividido em duas subseções. Na primeira é realizada uma
contextualização das turbinas hidráulicas e hidrocinéticas. A segunda subseção
discute os detalhes acerca do projeto hidrodinâmico das turbinas hidrocinéticas e os
métodos de solução das simulações numéricas computacionais, seja por algoritmos,
seja dinâmica dos fluidos computacionais.
2.1 TURBINAS HIDRÁULICAS E HIDROCINÉTICAS
O princípio de funcionamento das turbinas hidráulicas consiste na conversão
da energia potencial contida na água, em energia mecânica por meio da velocidade
rotativa do rotor. De acordo com Brekke(2001), as turbinas hidráulicas podem ser
divididas em dois grandes grupos:
i. Turbinas de reação, onde existe uma diferença de pressão entre a
entrada e a saída do rotor. Nesse tipo de turbina enquanto o fluido
percorre seu caminho através do rotor ocorre variação tanto em sua
energia cinética quanto em sua energia de pressão.
ii. Turbinas de impulsão, onde não existe diferença de pressão entre a
entrada e a saída do rotor. Nesse tipo de turbina ocorre variação
apenas na energia cinética do fluido durante sua passagem pelo rotor
da máquina. Ou seja, o rotor de turbinas de impulsão trabalha a
pressão constante.
Considerando aplicações comerciais, os principais tipos de turbinas
hidráulicas utilizados atualmente são as turbinas Pelton (turbinas de impulsão),
Francis e Kaplan (ambas de reação). Como o tema turbinas hidráulicas é muito amplo,
diversos outros parâmetros podem ser utilizados para subclassificar os tipos de
turbinas. Dentre estes é possível destacar a direção principal do escoamento através
do rotor. Considerando esse quesito as turbinas hidráulicas podem ser
subclassificadas em turbinas de fluxo tangencial, fluxo radial, fluxo misto e fluxo axial,
como representado na Figura 3.
23
Figura 3. Representação de uma turbina Pelton com fluxo tangencial (esquerda) e, à direita de cima para baixo: turbinas com fluxo radial, misto e axial.
Fonte: Adaptado de White (2011).
A alteração na direção principal do fluxo ao longo do rotor promove a
adaptação da turbina às diferentes condições de projeto. Turbinas de fluxo tangencial
são indicadas para altura de carga elevada (desnível entre o reservatório superior e
inferior) e pequenas vazões. Turbinas de fluxo radial operam com máxima eficiência
para alturas de carga e vazões intermediárias, ao passo que turbinas de fluxo axial
são indicadas para vazões mais altas e pequeno desnível geométrico entre os
reservatórios.
A primeira turbina hidráulica de alta eficiência foi construída por Benoit
Fourneyron entre 1824 e 1827. Essa turbina possuía fluxo radial na saída e uma
eficiência máxima da ordem de 85% (HEINLOTH, 2006). Ao longo de quase dois
séculos diversos tipos diferentes de turbinas foram construídos, no entanto, a maior
parte das turbinas testadas deixaram de ser utilizadas, inclusive a própria turbina de
Fourneyron. Ainda de acordo com Heinloth (2006), de um ponto de vista comercial
atualmente apenas três tipos de turbinas hidráulicas são largamente empregados
devido à sua superioridade em determinadas alturas de carga, sendo essas:
i. Turbina Pelton turbinas hidráulicas de fluxo tangencial;
ii. Turbina Francis turbinas hidráulicas de fluxo radial (compreendendo
desde turbinas radiais puras, Francis lenta, até turbinas praticamente
de fluxo misto, Francis ultrarrápida);
iii. Turbina Kaplan turbinas de fluxo axial, as quais podem ser de eixo
vertical ou de eixo horizontal (Kaplan bulbo).
24
Uma representação gráfica do uso apropriado dos principais tipos de turbinas
como função da velocidade específica pode ser vista na Figura 4.
Figura 4. Tipo de turbina ideal segundo altura de coluna de água (eixo vertical) e velocidade
específica (eixo horizontal)
Fonte: Adaptado de Menny (2006).
A forma construtiva da maior parte das turbinas hidrocinéticas desenvolvidas
até o momento se aproxima bastante da forma das turbinas hidráulicas de fluxo axial
convencionais. Assim, é nesta classificação que as turbinas hidrocinéticas se
encaixam (KAUNDA; KIMAMBO; NIELSEN, 2014).
Essa tecnologia é voltada à aplicações de pequena escala, com geração
próxima de 10 kW, e utiliza principalmente as velocidades de corrente de rios, ondas
e correntes marítimas (VERMAAK; KUSAKANA; KOKO, 2014). Algumas formas
construtivas podem ser vistas na Figura 5. Elas são compostas basicamente pelo
rotor, a caixa de câmbio, a nacele e o gerador, enquanto em algumas aplicações se
usa difusores (OBLAS, 2016).
25
Figura 5. Diferentes formas construtivas de turbinas hidrocinéticas horizontais de fluxo axial.
Fonte: Adaptado de (KHAN et al., 2009; VERMAAK; KUSAKANA; KOKO, 2014).
O rotor é o formado pelo cubo e pelo conjunto de pás, e tem como função a
geração da velocidade rotativa e a transmissão dessa para o eixo. O diâmetro do cubo
pode ser determinado tanto com base experimental, como visto na Figura 6, quanto
por formulação teórica, segundo as Equações (1) e (2) (PENCHE; ESHA, 2004). As
pás são construídas para que absorvam a maior quantidade de energia possível do
escoamento (HALDER; RHEE; SAMAD, 2017). Para isso, devem combinar a maior
relação entre sustentação e arrasto, e ao mesmo tempo não exceder os limites
estruturais e de desgaste, diretamente relacionados com a cavitação.
26
Figura 6. Razões idealizadas entre diâmetro do rotor e diâmetro do cubo, segundo velocidade específica.
Fonte: Menny (2006).
𝐷𝐻𝑢𝑏 = (0.25 + 0.0951
𝑁𝑠𝑝𝑒𝑐)𝐷 (1)
𝑁𝑠𝑝𝑒𝑐 = 𝜔√
(2𝑔𝐻)3/4 (2)
Onde Nspec é a velocidade específica adimensional, D o diâmetro do rotor, ω
a velocidade de rotação da turbina, a vazão volumétrica, g a gravidade e H a coluna
de água disponível; ainda quanto a Figura 6, δ representa a coluna de fluido
adimensional. Nesse sentido, a sua seção transversal são perfis aerodinâmicos, os
chamados aerofólios – ou, neste caso, hidrofólios. A caixa de câmbio é responsável
por sincronizar a velocidade de rotação do eixo com a do gerador. Este, por sua vez,
é quem converte a energia mecânica em energia elétrica. Por fim, a nacele é
incumbida de garantir a fixação da turbina e/ou do difusor.
27
Ao considerar o movimento de um fluido sobre a superfície de um aerofólio,
ou hidrofólio, é esperado que a distribuição de pressão se aproxime da distribuição de
pressão em fluido ideal (invíscido). Para números de Reynolds, eq. (3), elevados os
efeitos da viscosidade são confinados a uma região muito fina adjacente à fronteira
sólida. O fluido se adere a parede e as forças de fricção retardam o movimento do
fluido em uma fina camada próxima à parede. Nessa fina camada a velocidade do
fluido aumenta de zero na parede (condição de não deslizamento) até o valor
correspondente à velocidade da corrente livre. A camada em questão é chamada de
camada limite (SCHLICHTING; KESTIN, 1979).
𝑅𝑒 = 𝐿𝑈0𝜌
𝜇 (3)
Onde L é o comprimento característico, U0 a velocidade de corrente livre, ρ a
massa específica do fluido e μ a viscosidade dinâmica.
O escoamento na camada limite, sobretudo na região adjacente à parede,
possuí baixa quantidade de movimento. Em alguns casos a camada limite aumenta
consideravelmente sua espessura a jusante e o escoamento no interior da mesma
passa a se mover na direção reversa. Isso faz com que essa parcela de fluido com
baixa quantidade de movimento seja forçada para longe da parede. Esse fenômeno é
chamado de separação da camada limite e está sempre associado à geração de
vórtices e consequentemente a uma grande perda de energia na esteira viscosa
formada (SCHLICHTING; KESTIN, 1979). Essa perda de energia pode ser
relacionada à variação na distribuição de pressão em relação ao caso do escoamento
de um fluido invíscido o chamado arrasto de pressão ou de forma.
A separação da camada limite pode acontecer sempre que existir um
gradiente adverso de pressão. Isso normalmente ocorre no escoamento ao redor de
corpos rombudos como, por exemplo, esferas ou cilindros, no entanto, também é
esperado em escoamentos ao redor de aero/hidrofólios quando o ângulo de ataque
(este inclusive representado por α na Figura 7) se torna suficientemente elevado.
Inicialmente o incremento do ângulo de ataque é benéfico, promovendo o incremento
da sustentação. Porém, em determinado ângulo de ataque, a perda de contato devido
à separação é tamanha que gera o chamado estol, que é a queda brusca da
sustentação.
28
Figura 7. Representação do efeito de ângulo de ataque às forças de arrasto e sustentação do
perfil aero/hidrodinâmico.
Fonte: Adaptado de Hansen (2013).
Por outro lado, a medida que o número de Reynolds do escoamento aumenta
pode ocorrer a transição para a turbulência. A camada limite turbulenta possuí mais
energia e assim pode resistir a um gradiente adverso de pressão de maior intensidade
do que uma camada limite laminar. No entanto, a transição também promove um
aumento no arrasto de fricção devido a um gradiente de velocidade elevado.
Os perfis aerodinâmicos são projetados com intuito de promover a maior
sustentação (positiva ou negativa) e o menor arrasto possível para um dado
escoamento. Normalmente as informações a respeito do arrasto e sustentação de um
determinado perfil são apresentadas como coeficientes adimensionais de arrasto e de
sustentação, segundo as equações (4) e (5), onde 𝐴𝑠 é a área superficial do aerofólio.
O comportamento dos coeficientes de arrasto e sustentação é próprio de cada perfil,
e depende tanto do número de Reynolds como do ângulo de ataque do escoamento:
conforme aumenta, os coeficientes aumentam, já que o perfil aerodinâmico permite a
criação de maiores gradientes de velocidade e pressão entre as superfícies superior
e inferior, tendo como limite para sustentação o ângulo de estol.
𝐶𝐿 = 2 ∗ 𝐹𝑆𝑢𝑠𝑡𝑒𝑛𝑡𝑎çã𝑜
𝜌𝐴𝑠𝑈02 (4)
𝐶𝐷 = 2 ∗ 𝐹𝐴𝑟𝑟𝑎𝑠𝑡𝑜
𝜌𝐴𝑠𝑈02 (5)
29
A família de perfis NACA (do inglês U.S. National Advisory Committee on
Aeronautics) de quatro dígitos foi desenvolvida pela agência hoje conhecida por NASA
(do inglês National Aeronautics and Space Administration) no início dos anos 1920, e
contribuiu para o aprimoramento de asas e hélices da indústria aeronáutica após a
Primeira Guerra Mundial (WOOD, 2011a). Segundo pesquisas anteriores, alguns dos
perfis desse tipo mais indicados para turbinas hidrocinéticas são: NACA 0015
(HALDER; RHEE; SAMAD, 2017), NACA 4412 (WANG; YIN; YAN, 2019),NACA 4415,
NACA 4418 e NACA 4424 (ABUAN; HOWELL, 2019) e por fim NACA
6412(HOGHOOGHI; DURALI; KASHEF, 2018). A forma desses aerofólios pode ser
vista na Figura 8.
A diferença da sua nomenclatura é intrínseca a sua construção: o primeiro e
segundo dígitos referemse a percentagem de abaulamento (cambagem) e a sua
posição máxima, respectivamente e ambos com relação à corda; os dois últimos
dígitos se referem à espessura máxima, também tendo a corda como referência.
Figura 8. Alguns dos perfis NACA quatro dígitos empregados.
Fonte: Próprio autor.
30
2.2 SIMULAÇÃO NUMÉRICA COMPUTACIONAL
2.2.1 Principais desenvolvedores de teorias relacionadas aos rotores
O principal intuito do desenvolvimento das teorias relacionadas aos rotores, é
quanto sua performance, o plano do rotor, e seu efeito no escoamento, as chamadas
esteiras.
2.2.1.1 RankineFroude
Rankine (1865) e Froude(1889) foram os pioneiros em publicar trabalhos
referentes à compreensão de turbinas e hélices. Suas teorias consistiam na análise
de momento axial e unidimensional do escoamento contra um “obstáculo”, chamado
de disco atuador, como pode ser visto na Figura 9. Ou seja, não investigava os efeitos
de esteira e de velocidades rotativas, além de considerar a queda de pressão e a
velocidade axial uniformes ao longo da área (BRANLARD, 2017).
Figura 9. Exemplo esquemático da teoria formulada por Rankine e Froude.
Fonte: Adaptado de Branlard (2017).
Considerando dois volumes de controle, antes e depois do disco, e utilizando
o Teorema de Transporte de Reynolds (TTR) para conservação da massa, momento
31
e energia, encontrase as seguintes relações (6) e (7) para potência e empuxo,
respectivamente:
𝑃 = 1
2𝜌𝐴𝑈0
3𝐶𝑃 (6)
𝑇 = 1
2𝜌𝐴𝑈0
2𝐶𝑇 (7)
Onde 𝐴 é a área transversal do disco e CP e CT os coeficientes adimensionais
de potência e empuxo, respectivamente. Esses foram descritos como a fração de
energia retirada do escoamento e que foram convertidos em potência e empuxo Eq.
(8) e (9) (BRANLARD, 2017).
𝐶𝑃 = 4𝑎(1 − 𝑎)2 (8)
𝐶𝑇 = 4𝑎(1 − 𝑎) (9)
Onde 𝑎 é o fator de indução axial no plano do disco. Esse fator representa o
obstáculo que as pás oferecem ao escoamento, ou seja, o efeito axial que a turbina
oferece à velocidade de corrente livre. Ele pode ser determinado apenas tendo
conhecimento do escoamento resultante do rotor, como pode ser visto na Eq. (10).
𝑎 = 𝑈0 − 𝑉
𝑈0 (10)
Onde V é a velocidade induzida pelo disco. Devido às suas limitações
dimensionais, tal teoria não leva em conta o rotor e seus parâmetros de arranjo na
turbina, como por exemplo, número de pás e elas próprias, e fazem com que o
empuxo, torque, e assim, a potência, cresçam sem precedentes. Por isso, teorias que
consideram os efeitos do fluido na sua estrutura – e que consequentemente a fazem
ter o movimento de rotação – foram desenvolvidas em busca de um melhor
entendimento de sua construção (WOOD, 2011b).
32
2.2.1.2 FroudeDrzewiecki (Teoria do Elemento da Pá)
A Teoria do Elemento da Pá (ou Blade Element Theory BET, em inglês), foi
introduzida por Froude (1878) e Drzewiecki; GauthierVillars (PARYŻ) (1892). Essa
teoria consiste em dividir o volume de controle (pá) em elementos diferenciais e
bidimensionais, e ao aplicar as mesmas leis de conservação da massa, momento
angular e axial, e energia, possibilita analisar os efeitos de sustentação e arrasto no
rotor, com intuito de determinar os valores de corda, torção e velocidade (WOOD,
2011b) – detalhe que essa abordagem passa a tratar as pás do rotor, e não mais em
um disco, como esquematicamente exemplificado na Figura 10.
Ainda, considerase que a velocidade do fluido é conhecida e não varia
radialmente e que as forças que atuam em cada seção são independentes uma das
outras, com valores de arrasto e sustentação bidimensionais (BRANLARD, 2017).
Essa última só é possível pela hipótese do perfil visto como “linha de sustentação”
(BRANLARD, 2017), em que a largura é infinitesimal e tanto a velocidade como as
forças podem ser encontradas.
Figura 10. Representação esquemática da divisão do volume de controle (a); triângulo de
velocidades, a representação das forças e velocidades atuantes em cada elemento da pá (b).
(a) (b)
Fonte: (a) Adaptado de Wood (2011b); (b) Adaptado de Branlard (2017).
O principal resultado de tal teoria é a consideração das forças como agentes
de empuxo e torque, como apresentado nas Eq.(11) e (12):
𝑑𝑇 = 1
2𝜌𝑈0
2𝑐𝐵𝐶𝑛𝑑𝑟 (11)
𝑑𝑄 = 1
2𝜌𝑈0
2𝑐𝐵𝐶𝑡𝑟𝑑𝑟 (12)
33
Onde c é a corda, B o número de pás, Cn o coeficiente de força normal, Ct o
coeficiente de força tangencial e r o raio local do rotor. Esses dois últimos são descritos
nas Eq.(13) e (14):
𝐶𝑛 = 𝐶𝐿 cosΦ + 𝐶𝐷 sinΦ (13)
𝐶𝑡 = 𝐶𝐿 sinΦ − 𝐶𝐷 cosΦ (14)
Onde Φ é o ângulo do escoamento. Pela Figura 10b, determinase esse
ângulo com as Eq.(15) e (16):
Φ = 𝛼 + 𝜃 (15)
tanΦ = (1 − 𝑎)
(1 + 𝑎′)𝜆𝑟 (16)
Em que α é o ângulo de ataque, θ o ângulo de torção da pá, 𝑎′ o fator de
indução tangencial da pá no plano do rotor e λr a razão local de velocidade tangencial
da pá e da corrente livre, segundo Eq.(17). Tal fator de indução tangencial representa
o obstáculo do escoamento relacionado à rotação da esteira. Devese ressaltar que
segundo a teoria de RankineFroude, o valor do fator tangencial é nulo, uma vez que
a rotação do disco é desconsiderada.
𝜆𝑟 = 𝑟𝜔
𝑈0 (17)
2.2.1.3 Prandtl e Betz
A investigação para a correção do número infinito de pás foi iniciada por
Prandtl e Betz (PRANDTL; BETZ, 1927). Foi desenvolvido o fator de correção nas
pontas das pás (ou do inglês TipLoss Factor, F), que procura corrigir os parâmetros
para um número finito como pode ser visto na Figura 11. Da sua definição trivial, a
razão da circulação local pela circulação da corrente livre, faz com que esse efeito se
estenda ao longo do comprimento da pá; sua solução foi investigada por Goldstein;
34
Prandtl (1929) e tem como consequência a geração de um formato de esteira
helicoidal. A formulação inicial é apresentada na Eq. (18) (PRANDTL; BETZ, 1927).
Figura 11. Comparação entre os modelos de escoamento nas pontas das pás, para um número
infinito, adotado por Betz (a), e finito considerado por Prandtl (b)
(a) (b)
Fonte: Adaptado de Branlard (2017).
𝐹 = 2
𝜋cos−1 [exp (−
𝐵(𝑅 − 𝑟)√1 + 𝜆2
2𝑅)] (18)
Onde R é o raio do rotor e λ é a razão de velocidade tangencial da ponta da
pá pela velocidade de corrente livre, segundo Eq. (19).
𝜆 = 𝑅𝜔
𝑈0 (19)
Quando se trata da conversão de energia do escoamento, se considerarmos
o volume de controle tal como a Figura 9 – disco sem restrição, em que a própria
circunferência do rotor é a fronteira , o coeficiente de potência terá seu maior valor se
derivarmos sua função em relação ao fator de indução axial. Esse ponto é conhecido
por Limite de Betz, um de seus desenvolvedores (BETZ, 1920), e é igual a 0,593 para
um valor de 𝑎 em 1/3. Isso quer dizer que, teoricamente, a percentagem máxima de
extração da energia cinética do fluido na condição citada é 59,3%.
35
2.2.1.4 Glauert (Teoria do Momento da Pá)
A partir da combinação das duas últimas teorias, Glauert (1935) formulou sua
teoria de rotor ótimo, a mais generalizada (pela Figura 12, vêse que quase todos os
parâmetros são considerados, como o fator de indução axial e os coeficientes de
potência e torque) e utilizada atualmente em projetos de turbinas, conhecida como
Teoria do Momento da Pá (ou em inglês, Blade Element Momentum BEM). Consiste
em determinar o triângulo de velocidades provenientes da BET, diretamente afetados
pelos fatores de indução da teoria de RankineFroude (Figura13); assim, é possível
determinar o ângulo Φ e os esforços em cada seção da pá. Consequentemente,
determinase o empuxo e o torque segundo as Eq. (11) e (12); assumindo uma
distribuição uniforme da força no diferencial de área, essas são então igualadas às
Eqs. (6) e (7) (BRANLARD, 2017). Deste último passo, encontrase novos valores
para os fatores de indução, e assim continua o loop de convergência para o método
iterativo.
36
Figura 12. Representação do volume de controle considerado.
Fonte: Adaptado de vanKuik; Sørensen; Okulov (2015).
Figura13. Nova representação do triângulo de velocidades, considerado para o BEM.
Fonte: Adaptado de Branlard (2017).
Na tentativa de considerar o número de pás, Glauert (1935) adaptou à Eq. (8)
o fator de perdas nas pontas da pá, e sua nova versão é apresentada na Eq. (20).
𝐶𝑃 = 4𝑎𝐹(1 − 𝑎)2 (20)
O mesmo foi feito com relação ao coeficiente de empuxo. Entretanto, em seus
trabalhos experimentais, Glauert (1935) encontrou desconformidades para altos
valores de fator de indução axial. Essa situação é contraditória à teoria BEM, uma vez
que é válida para pequenas alterações da velocidade de esteira (BRANLARD, 2017).
Dessa forma, propôs correções segundo a Eq. (21).
𝐶𝑇 = 4𝑎𝐹(1 − 𝑎) 𝑝𝑎𝑟𝑎 𝑎 ≤
1
3
𝐶𝑇 = 4𝑎𝐹 (1 −1
4(5 − 3𝑎)𝑎) 𝑝𝑎𝑟𝑎 𝑎 >
1
3
(21)
37
Após algumas considerações, também foi possível determinar os fatores de
indução axial e tangencial no plano do rotor, apresentados nas Eq. (22) e (23),
respectivamente:
𝑎 = 1
4𝐹sin²Φ
𝜎𝐶𝑛 + 1
(22)
𝑎′ = 1
4𝐹 sinΦcosΦ
𝜎𝐶𝑡 − 1
(23)
Onde 𝜎 é a solidez da turbina – apresentada de forma mais detalhada na
seção 2.2.1.9. Além disso, o autor também reformulou o fator de correção nas pontas
das pás, inicialmente apresentado na Eq. (18), para a Eq. (24).
𝐹 = 2
𝜋cos−1 [exp (−
𝐵(𝑅 − 𝑟)
2𝑟 sinΦ)] (24)
Sendo essa teoria a mais utilizada quando se trata de projetos de rotores, seja
para turbinas hidráulicas ou eólicas, alguns autores fizeram uso dela para otimizações,
abordagem parecida com a proposta neste trabalho, e tiveram resultados
interessantes: Riglin et al. (2016) e Schleicher; Riglin; Oztekin (2015) trataram de
turbinas hidrocinéticas quase idênticas ao tipo Kaplan, enquanto que Wang; Yin; Yan
(2019) outra mais parecida com turbinas eólicas.
2.2.1.5 Wilson e Lissaman
Por mais que a corrente livre seja irrotacional, sua interação com uma
máquina girante fará com que ela rotacione, seja no mesmo sentido de rotação da
máquina, comum em hélices, ou no contrário, caso das turbinas (SPERA, 2009). Esse
efeito proporcionado ao escoamento, a rotação da esteira, é significativo à sua
performance.
Isso quer dizer que além de levar em conta a energia angular no plano do
disco e a variação axial da velocidade, o âmbito da esteira também deve ser analisado.
38
O diferencial de empuxo e de torque para esse caso são descritos por Wilson;
Lissaman (1974) e apresentados nas equações (25) e (26).
𝑑𝑇 = 𝜌 (𝜔 +𝛺
2) 𝑟2𝜔𝑑𝐴 (25)
𝑑𝑄 = 𝜌𝑈0𝑟2 𝛺𝑑𝐴 (26)
Onde 𝛺 é a velocidade angular do fluido no plano do disco, cuja consideração
ainda se limita pelos efeitos rotativos na esteira. Já seu coeficiente de potência
desenvolvida é apresentado na Eq. (27).
𝐶𝑃 = 𝑏(1 − 𝑎)2
𝑏 − 𝑎[𝑏 +
(2𝑎 − 𝑏)𝜔
𝛺𝑚𝑎𝑥] (27)
Onde 𝑏 é o fator de indução axial no plano da esteira. Partindo do
entendimento de que a circulação é essencialmente causada pela força de
sustentação da pá, Wilson; Lissaman (1974) a incluiu nas formulações, apresentadas
nas Eq. (28) e (29). Devese destacar que o fator de indução tangencial do plano do
rotor é praticamente igual ao desenvolvido por Glauert (1935) na teoria BEM.
𝑎𝐹(1 − 𝑎𝐹)
(1 − 𝑎)²= 𝜎𝐶𝐿cosΦ
4sin²Φ (28)
𝑎′𝐹
(1 + 𝑎′)=
𝜎𝐶𝐿 4cosΦ
(29)
2.2.1.6 Spera
Ainda ao considerar a rotação da esteira, o fator tangencial é impossível de
ser aproximado sem conhecimento do escoamento por si só, já que a velocidade
rotacional ao longo da pá é constante. Em contrapartida, o fator axial no plano do rotor
pode ser aproximado pela Eq. (30) em função do fator de indução axial na esteira
(SPERA, 1994).
39
𝑎 = 𝑏
2(1 −
(1 − 𝑎)𝑏²
4𝜆²(𝑏 − 𝑎)) (30)
Tavares Dias do Rio Vaz et al. (2013), em seu trabalho, fizeram uso da
consideração dos efeitos de esteira junto da teoria BET com intuito de encontrar um
design de turbina ideal. Ainda, Spera propôs a Eq. (31).
𝐶𝑇 = 4𝑎𝐹(1 − 𝑎) 𝑝𝑎𝑟𝑎 𝑎 ≤ 𝑎𝐶𝐶𝑇 = 4𝐹(𝑎𝐶
2 + (1 − 2𝑎𝐶)𝑎) 𝑝𝑎𝑟𝑎 𝑎 > 𝑎𝐶 (31)
Onde 𝑎𝐶 é o valor crítico de 𝑎. Para o autor, este valor é igual a 0,2, enquanto
Wilson; Lissaman (1974) encontrou 0,37
2.2.1.7 Shen
Shen et al. (2005) salientou que pelo BEM se tratar de uma teoria
bidimensional as forças não se anulam nas pontas, e isso representa uma
inconsistência física, uma vez que elas deveriam tender à zero por conta do fluido se
transportar diretamente do lado pressurizado da pá para o de sucção. Por isso, propôs
uma correção tridimensional às forças envolvidas nos fatores, como é apresentado
nas Eq. (32) e (33).
𝑎𝐹(1 − 𝑎𝐹)
(1 − 𝑎)²=
𝜎𝐶𝑡𝐹1
4𝐹sin²Φ (32)
𝑎′𝐹(1 − 𝑎𝐹)
(1 − 𝑎)(1 + 𝑎′)=
𝜎𝐶𝑛𝐹1 4sinΦcosΦ
(33)
Onde F é o fator de perdas na ponta da pá proposto por Glauert (Eq. (24)) e
F1 o proposto por Shen (Eqs. (34) e (35)). Esse último por sua vez foi modelado de
forma semiempírica, que inclusive também teve como resultado correlações para o
coeficiente de empuxo(36).
𝐹1 = 2
𝜋cos−1 [exp (−𝑔
𝐵(𝑅 − 𝑟)
2𝑟 sin𝛷)] (34)
40
𝑔 = exp [−1
8(𝐵𝜆 − 21)] + 0.1 (35)
𝐶𝑇 = 4𝑎𝐹(1 − 𝑎)𝑝𝑎𝑟𝑎 𝑎 ≤ 𝑎𝐶 =
1
3
𝐶𝑇 = 4[𝑎𝐶2𝐹² + (1 − 2𝑎𝐶𝐹)𝑎𝐹]𝑝𝑎𝑟𝑎 𝑎 > 𝑎𝐶 =
1
3
(36)
2.2.1.8 Demais correções
2.2.1.8.1 Fatores de indução axial e tangencial
Vries (1979) apontou uma inconsistência nas Eq. (28) e (29), quanto ao
resultado não ortogonal entre a velocidade induzida e a velocidade relativa. Contra
isso, corrigiu apenas o fator tangencial (Eq. (37)).
𝑎′𝐹(1 − 𝑎𝐹)
(1 − 𝑎)(1 + 𝑎′)=
𝜎𝐶𝐿 4sinΦ
(37)
Sob uma perspectiva experimental, Buhl (2005) propôs uma correção
segundo a Eq. (38) para o valor do fator de indução axial quando este for maior ou
igual a 0,4 (pelos altos valores de expansão da esteira), enquanto Jonkman (2003)
apresentou a Eq. (39) referente ao fator de indução tangencial. Ambos foram utilizados
e investigados por Lanzafame; Messina (2013).
𝑎 = 18𝐹 − 20 − 3√𝐶𝑇(50 − 36𝐹) + 12𝐹(3𝐹 − 4)
36𝐹 − 50 (38)
𝑎′ = 1
2(√1 +
4
𝜆𝑟2𝑎(1 − 𝑎) − 1) (39)
Ainda, em tese, um maior valor de fator de indução axial no rotor representa
uma maior absorção da energia contida no escoamento. Entretanto, como visto por
Betz (1920), o valor ótimo do fator de indução axial é igual a 1/3, enquanto que para
o tangencial é zero. Por isso, desenvolveuse a Eq. (40) para que o valor do fator de
indução axial em função da razão local de velocidades seja aproximadamente igual a
41
1/3, e da mesma forma, para que a Eq. (41) forneça valores mínimos para o fator de
indução tangencial (HANSEN, 2013; BRANLARD, 2017).
16𝑎3 − 24𝑎2 + 3𝑎(3 − 𝜆𝑟2) − 1 + 𝜆𝑟
2 = 0 (40)
𝑎′ =1 − 3𝑎
4𝑎 − 1 (41)
2.2.1.8.2 Fator de perda nas pontas das pás
Ao longo da história foram definidas novas abordagens teóricas e
semiempíricas para o fator de correção na ponta da pá (como o próprio formulado e
tratado por Shen), como são relatados em Branlard (2017) e CliftonSmith (2009) e
apresentadas a seguir.
Moriarty; Hansen (2005) abordou essa perda não só na ponta das pás, como
também junto ao cubo. Dessa forma, desenvolveram as Eq. (42), (43) e (44).
𝐹𝑇𝑖𝑝 = 2
𝜋cos−1 [exp (−
𝐵(𝑅 − 𝑟)
2𝑟 sinΦ)] (42)
𝐹𝐻𝑢𝑏 = 2
𝜋cos−1 [exp (−
𝐵(𝑟 − 𝑅𝐻𝑢𝑏)
2𝑟 sinΦ)] (43)
𝐹 = 𝐹𝑇𝑖𝑝𝐹𝐻𝑢𝑏 (44)
Onde RHub é o raio do cubo. Já Burton et al. (2011) desconsiderou o fator de
indução tangencial e encontrou a relação apresentada na Eq. (45):
𝐹 = 2
𝜋cos−1 [exp(−
𝐵
2(𝜆
𝜆𝑟− 1)√1 + (
𝜆𝑟1 − 𝑎
)2
)] (45)
2.2.1.8.3 Coeficientes de empuxo e potência
A partir da formulação do BEM, diferentes abordagens para os coeficientes de
potência e empuxo são possíveis e consequentemente resultam em equacionamentos
distintos. Destacase que pela natureza da teoria, os valores de CP são locais ao longo
42
da pá, e ao final, a integral no seu comprimento representa o coeficiente do rotor.
Dentre elas, apresentase a Eq. (46) (BRANLARD, 2017), que desconsidera a rotação
da esteira, enquanto que as Eq. (47) e (48) trazem as propostas de Letcher (2017) e
Jonkman (2003) o fator de perdas nas pontas da pá já é considerado implicitamente.
𝐶𝑃 = 8
𝜆2∫ 𝑎′
𝜆
𝜆ℎ𝑢𝑏
(1 − 𝑎)𝜆𝑟3𝑑𝜆𝑟 (46)
𝑑𝐶𝑃 = 𝐵𝜆2(1 − 𝑎)(1 + 𝑎′)
𝑐
𝑅∗𝑟
𝑅𝐶𝑡
2𝜋 sinΦ cosΦ (47)
𝐶𝑃 = 2
𝜆𝑅∫ 𝜆𝑟
2sin²Φ
𝑅
𝐻
[cosΦ − 𝜆𝑟 sinΦ][sinΦ + 𝜆𝑟 cosΦ] [1 −𝐶𝐷 cosΦ
𝐶𝐿 sinΦ]𝑑𝑟 (48)
Ao mesmo tempo outros pesquisadores também desenvolveram suas
correções ao coeficiente de empuxo, como foi discutido por Pratumnopharat; Leung
(2011). Manwell; McGowan; Rogers (2011) partiram da proposta de Glauert e
desenvolveram a correção apresentada na Eq. (49).
𝐶𝑇 = 4𝑎𝐹(1 − 𝑎) 𝑝𝑎𝑟𝑎 𝑎 ≤ 0.4
𝐶𝑇 = 0,96 + 𝐹(𝑎 − 0,4)[𝐹(𝑎 + 0,4) − 0,286]
0,6427 𝑝𝑎𝑟𝑎 𝑎 > 0.4
(49)
Buhl (2005) desenvolveu sua correção a partir de uma regressão quadrática,
segundo a Eq. (50).
𝐶𝑇 = 4𝑎𝐹(1 − 𝑎) 𝑝𝑎𝑟𝑎 𝑎 ≤ 0.4
𝐶𝑇 =8
9+ (4𝐹 −
40
9) 𝑎 + (
50
9− 4𝐹)𝑎2 𝑝𝑎𝑟𝑎 𝑎 > 0.4
(50)
2.2.1.9 Efeito cascata das pás
Ao revisar as considerações impostas pela BET, vêse que a que diz respeito
sobre uniformidade da força em todas as seções da pá é facilmente violada: ela
considera um espaço infinito de escoamento, enquanto que na realidade, está limitado
43
à distância que separa as pás, a chamada solidez (σ, ou, mais comum, em inglês,
solidity) (WOOD, 2011c). Essa é a razão entre o espaço ocupado pelas pás, e a
circunferência local do rotor. Tal relação é traduzida na equação (51) e representada
na Figura 14.
𝜎 =𝐵𝑐(𝑟)
2𝜋𝑟 (51)
Figura 14. Representação gráfica do parâmetro de solidez.
Fonte: Branlard (2017).
Do ponto de vista da dinâmica dos fluidos, a consideração feita até então da
“linha de sustentação” não prevê alterações no escoamento por conta da espessura
e largura dos perfis aerodinâmicos; entretanto, ao passo que a solidez passa a levá
los em conta, a distorção do fluxo pode ser significativa – como visto na Figura 15.
Esse efeito é chamado de “efeito cascata”.
44
Figura 15. Distorção da linha de corrente do escoamento, ao considerar o efeito cascata.
Fonte: Adaptado de Spera (2009).
Segundo Spera (2009), essa alteração leva à mudança das componentes
tangencial e axial da velocidade, e como consequência, varia o ângulo de ataque aos
perfis segundo as Eqs. (52), (53), (54) e (55). Sendo a Eq. (53) especificamente para
a família de perfis NACA com quatro dígitos.
∆𝛼 = ∆𝛼1 + ∆𝛼2 (52)
∆𝛼1 = 0.11𝐵 (𝑐
𝑅) (𝑡
𝑐)𝑚𝑎𝑥
𝜆
√(1 − 𝑎)2 + (𝑐
𝑅) ²𝜆²
(53)
∆𝛼2 ≈∆Φ
4 (54)
∆Φ = tan−1 [(1 − 𝑎)𝑥
(1 + 2𝑎′)] − tan−1[(1 − 𝑎)𝑥] (55)
Onde Δα1 e Δα2 são as variações do ângulo de ataque devido às alterações
da espessura e da largura, respectivamente, t é a espessura e x a coordenada axial
do perfil. Além disso, empiricamente o comportamento do escoamento se mostra
transformado, ao ponto de atrasar o ângulo de estol do perfil principalmente ao redor
do cubo, região em que o valor de solidez é maior; não se sabe as causas para tal,
mas é comumente associado à alteração do campo de pressão, já que o escoamento
é confinado entre as duas pás. Dixon (2005) investiga amplamente essas
consequências, porém, com base na compressibilidade do fluido e no conhecimento
45
prévio de parâmetros do escoamento (como o ângulo de saída do perfil), fugindo do
escopo do presente trabalho.
Utilizando DFC, Cebrián et al. (2013) fez a análise de arranjo de placas planas,
obtendo correlações para o ângulo de saída do fluido, enquanto Ahmed; Yilbas; Budair
(1998) discutiu os efeitos do escoamento em um perfil NACA 0012, sob um Reynolds
de 3,24E5 num intervalo de ângulo de ataque entre 0 e 24º, para solidez iguais a 0,83
e 0,55. Após a constatação de que os valores de solidez podem ter consequência nos
coeficientes, Fleck (2017) realizou uma série de coleta de dados na simulação de 3
seções transversais da pá otimizada de uma mini turbina eólica, cuja solidez do perfil
aerodinâmico SD7062 era igual a 0,1, 0,2 e 0,3, sob diferentes valores de Reynolds e
de razões de velocidade de corrente livre e de ponta da pá.
Além disso, MicleaBleiziffer et al. (2014) investigou o impacto que correções
dos coeficientes de arrasto e sustentação devido a valores de solidez maiores que 3
causam na performance de uma turbina com pás de perfil hidrodinâmico. Tal efeito de
valores altos de solidez é amplamente investigado em outros estudos sob a
justificativa de que o desempenho da turbina também pode ser afetado. Islam; Islam
(1994) comparou a previsão do funcionamento de moinhos de vento por teorias de
cascata e pela teoria de momento da pá; através da constatação da predominância
desse fenômeno nas operações com baixo valor de razão da velocidade de corrente
livre e de ponta da pá, Singh (2014) estendeu a predição corrigida pelo efeito cascata
para turbinas que possuem alto valor de solidez.
Ainda nesse sentido, quando diz respeito ao número de pás exclusivamente,
Duquette et al. (2003) modificou uma turbina existente, aumentando o número de pás,
com a mesma distribuição de corda, e verificou um decaimento da potência produzida.
Entretanto, Brasil Junior et al. (2019), investigou experimentalmente o coeficiente de
potência para uma hélice contendo duas, três e quatro pás, porém alterando ao
mesmo tempo a distribuição da corda; como resultado, obteve um aumento de ambos
os coeficientes, acompanhado de um decréscimo do tamanho de corda.
2.2.1.10 Cavitação
A cavitação consiste na mudança do estado líquido para vapor, por meio da
queda de pressão localizada além de um valor crítico; esse gradiente está associado
46
a efeitos dinâmicos (seja ele um escoamento ou um campo acústico), enquanto que
esse valor crítico é a pressão de vapor (ARNDT, 1981). O vapor formado fica
aprisionado na forma de bolhas, que podem ser succionadas para regiões de maior
pressão, onde elas se condensam, colapsando as bolhas. Isso gera uma onda de
choque cuja magnitude pode chegar a 400 MPa (DIXON, 2005), causando problemas
estruturais quando ocorrem próximas à superfícies sólidas. No intuito de prever tal
evento, se desenvolveu o índice de cavitação, apresentado na equação (56). Assim,
se este fator for maior que o termo crítico, então a operação está livre de cavitação.
𝜎𝐶𝐴𝑉 =𝑝0 − 𝑝𝑣1
2𝜌𝑈0
2 (56)
Onde p0 é a pressão de referência e pv a pressão de vapor do líquido na
temperatura de trabalho. A consequência direta desse fenômeno, que pode ocorrer
em qualquer máquina de fluxo de reação (um exemplo ilustrativo pode ser visto na
Figura 16), é a erosão superficial. Além de diminuir a vida útil da máquina, causando
falha por fadiga, a performance é afetada. Avellan (2004) detalhou a queda de
eficiência em bombas, turbinas do tipo Kaplan e Francis desta, inclusive, chegando
a uma variação de 3% (KUMAR; SAINI, 2010).
Se tratando de turbinas hidráulicas, a região de sucção das pás é a mais
suscetível, uma vez que há a criação de zonas de baixa pressão (DIXON, 2005).
Nesse sentido, a mínima pressão do fluxo sob o perfil da pá deve ser maior que a
pressão de vapor do fluido (ARNDT, 1981); isso pode ser traduzido pelo coeficiente
de pressão, que varia da mesma forma que os coeficientes de arrasto e sustentação.
47
Figura 16. Exemplo de bolhas decorrentes de cavitação em hélices marítimas.
Fonte: Arndt (1981).
Dentre as formas encontradas para se evitar tal evento, Kumar; Saini (2010)
aponta que o rotor deve ser projetado apropriadamente. Com este intuito, se utilizando
da teoria BEM em turbinas hidrocinéticas, Silva et al. (2017) propôs um design
otimizado de corda que evite a cavitação; os autores adaptaram a equação (56) para
a velocidade relativa do fluido e trataram a pressão de referência do fluido (p0) como
a soma da pressão atmosférica (patm) com a coluna de água local h, em função da
posição radial. Essas adequações levaram à uma velocidade característica para a
cavitação (57), em que ao se comparar com a velocidade relativa (58), a corda pode
ser corrigida ou não (59).
𝑉𝐶𝐴𝑉 = √𝑝𝑎𝑡𝑚 + 𝜌𝑔ℎ − 𝑝𝑣
−1
2𝜌𝐶𝑝𝑚𝑖𝑛
(57)
𝑊 = √[𝑈0(1 − 𝑎)]2 + [𝑟𝜔(1 + 𝑎′)]2 (58)
𝑐𝑢𝑐 =
2𝜋𝑟
𝐵
𝐶𝑇𝐶𝑛(𝑈0𝑊)2
𝑝𝑎𝑟𝑎 𝑊 < 𝑉𝐶𝐴𝑉
𝑐𝑐𝑜 = 𝑐𝑢𝑐 [𝑊
(1 − 𝑓𝑠)𝑉𝐶𝐴𝑉]2
𝑝𝑎𝑟𝑎 𝑊 ≥ 𝑉𝐶𝐴𝑉
(59)
Onde Cpmin o coeficiente de pressão mínimo da seção transversal da pá, cuc a
corda não corrigida e cco a corda corrigida. Para uma dada configuração de turbina ao
48
utilizar a simulação numérica como avaliação dessa abordagem Silva et al. (2017)
identificaram que a pá com corda corrigida produziu 12 vezes menos vapor do que a
não corrigida
2.2.1.11 Difusores
Até o momento, o desempenho das turbinas segundo empuxo e potência são
previstos pela teoria BEM; esta, parte da consideração de um volume de controle sem
restrição, em que seu limite é uma fronteira “invisível” de influência das pontas das
pás. Ao adicionar um difusor em volta do rotor, essa fronteira passa a ser física, o
volume de controle confinado à sua estrutura e a performance pode ser beneficiada.
Difusores são componentes que desaceleram o fluido aumentando sua
pressão (DIXON, 2005). A estrutura se baseia no aumento da área de seção
transversal, e a proporção entre as áreas de entrada e saída do difusor é um
importante fator de projeto: se for feito de maneira brusca, ou seja, em um
comprimento curto, a camada limite da parede tende a se separar diminuindo o
coeficiente de sustentação (HANSEN, 2013) – coeficiente esse, que quanto maior seu
valor, maior será a sucção de fluido, ocasionado pela circulação anelar da esteira
(HANSEN; SØRENSEN; FLAY, 2000); como também não pode ser feito em um longo
comprimento, já que passase a ter maior contato com a parede, e assim, maiores
perdas por atrito (DIXON, 2005).
Essa força gerada pela parede do difusor possibilita que, por exemplo, o
coeficiente de potência supere o limite de Betz (HANSEN, 2013). A razão entre os
coeficientes com e sem o difusor pode ser vista na equação (60) e um exemplo da
influência nas velocidades na Figura 17.
𝐶𝑃𝐷𝐶𝑃
=(1 − 𝑎)
(60)
Onde ε é a razão entre a velocidade de corrente livre e a velocidade no rotor
com difusor e 𝐶𝑃𝐷 é o coeficiente de potência com o difusor.
49
Figura 17. Ilustração dos efeitos do difusor nas velocidades.
Fonte: Adaptado de Tavares Dias do Rio Vaz et at. (2014).
Esse benefício do aumento de potência produzida é tema recorrente de
pesquisas voltadas às turbinas eólicas e hidráulicas. Investigouse o impacto dos
difusores nos parâmetros de funcionamento da turbina já descritos anteriormente,
tanto nas teorias clássicas unidimensionais e BEM (TAVARES DIAS DO RIO VAZ et
al., 2014), quanto na teoria BET (VAZ; WOOD, 2016). Sob um novo triângulo de
velocidades, foram adaptadas as equações (16), (46), (22) e (40), para as (61), (62),
(63) e (64). Essas análises, entretanto, levam em consideração as perdas no difusor
apenas se tratando da performance do rotor; em uma exposição mais a fundo, Vaz;
Wood (2018) considerou essa eficiência nos parâmetros de operação.
tanΦ = 𝑎′𝜆𝑟𝛾𝑎
(61)
𝑎
(1 − 𝑎)=
𝛾²𝜎𝐶𝑛
4𝐹sin²Φ (62)
16𝑎3 − 24𝑎2 + 𝑎 (9 − 3 (𝜆𝑟𝛾)2
) − 1 + (𝜆𝑟𝛾)2
= 0 (63)
𝐶𝑃 = 8
𝜆2∫𝛾𝑎′
𝜆
0
(1 − 𝑎)𝜆𝑟3𝑑𝜆𝑟 (64)
50
Figura 18. Representação do perfil de velocidade na seção transversal da pá, e consequentemente, a razão de aceleração da velocidade.
Fonte: Adaptado de Vaz et al. (2019).
Onde γ é a razão de aumento da velocidade no difusor, como pode ser visto
na Figura 18. Com tal base teórica, do Rio Vaz; Vaz; Silva (2018) otimizou o fator de
indução axial (através da Eq. (63)) em uma turbina hidrocinética ampliada por um
difusor, incluindo uma correção para evitar cavitação; obteve um incremento de 42%
na geração de potência, mas em contrapartida, um aumento de quase 64% da corda.
Já Vaz et al. (2019), fazendo uso da teoria BEM, das Eq. (61), (62) e (64) e de modelos
dinâmicos mecânicos, discutiu o aumento das perdas, da potência gerada e da
velocidade rotativa em turbinas eólicas e hidrocinéticas; os autores concluíram que
mesmo havendo um aumento nas perdas devido ao arrasto a rotação aumentou em
20% e houve um ganho na potência gerada de quase 40%.
Esses coeficientes são dependentes exclusivamente da geometria do difusor
e geralmente obtidos experimentalmente, uma vez que o comportamento do fluido é
difícil de se prever analiticamente. Neste cenário, inúmeras tratativas experimentais
ou de se utilizar DFC para os adquirir são investigadas. Em uma abordagem mista,
Barbosa et al. (2015) desenvolveu uma relação analítica para o perfil de velocidade
do fluido ao longo do comprimento e da altura do difusor, cujos coeficientes variam de
acordo com a geometria e perfil do mesmo. Já Gaden; Bibeau (2010) utilizando DFC,
investigaram o impacto da relação entre áreas de entrada e saída, e o ângulo de
abertura em difusores com formato plano, para turbinas hidrocinéticas; como
resultado, a relação de 1,56 e um ângulo de 20º produziram 3,1 vezes mais potência
que uma turbina sem esse componente. Analogamente, Song; Wang; Yan (2019)
validou as simulações numéricas experimentalmente, de uma turbina hidrocinética
com difusor de perfil hidrodinâmico NACA 0012, cujo aumento de potência foi de, em
média, 45%.
51
2.2.2 Dinâmica dos fluidos computacional
Resolver problemas relacionados à mecânica dos fluidos através do TTR nos
dá informações gerais quanto ao seu desempenho, por conta da sua abordagem
integral. Em contrapartida, os detalhes acerca dos parâmetros de escoamento e o
campo gerado por eles não são obtidos. Para que isso seja possível, o volume de
controle, passa a ser dividido e tratado como um volume diferencial. Assim, as porções
integrais da equação se reduzem ao próprio diferencial da variável. Dessa forma, a
conservação da massa para escoamentos incompressíveis passa a ser representada
pela equação (65), mais conhecida como equação da continuidade.
𝜕𝑢𝑖𝜕𝑥𝑖
= 0 (65)
Onde u é a componente da velocidade e x a coordenada cartesiana,
caracterizadas pelo índice i, o qual varia de um até três. Essa formulação permite
relacionar a massa específica com os campos de velocidade. Esses por sua vez,
podem ser encontrados quando a 2ª lei de Newton é aplicada ao TTR, obtendo assim
o equacionamento para quantidade de movimento. Ao considerar que o fluido em
questão é newtoniano, se obtém as equações de NavierStokes (66) nesse caso
escritas para um escoamento incompressível.
𝜌 [𝜕𝑢𝑖𝜕𝑡
+ 𝑢𝑗𝜕𝑢𝑖𝜕𝑥𝑗
] = 𝜌𝑔𝑖 −𝜕𝑃𝑖𝜕𝑥𝑖
+𝜕
𝜕𝑥𝑗(𝜇𝜕𝑢𝑖𝜕𝑥𝑗
) (66)
Devese destacar que essa formulação relaciona o efeito das forças externas,
sejam de superfície (tensão cisalhante e gradiente de pressão) e/ou de corpo
(gravidade), com o campo de velocidade do escoamento.
No total, as equações (65) e (66) envolvem quatro variáveis a serem
resolvidas de forma parcial (ou seja, dependem entre si) e para cada um dos volumes
diferenciais de controle. Sem hipóteses a serem consideradas, de forma geral isso
inviabiliza a solução analítica. Neste contexto, o uso de ferramentas computacionais
para a solução dessas e outras equações se mostra uma solução viável, e que
caracteriza a dinâmica dos fluidos computacional.
52
2.2.2.1 Modelos de turbulência
Escoamentos turbulentos podem ser vistos como uma superposição de
flutuações aleatórias (irregulares) sobre uma corrente principal. Essas flutuações são
capazes de gerar um elevado efeito de mistura no escoamento, levando a uma alta
difusividade e alta dissipação. Em alguns casos esse efeito pode ser equiparado a um
aumento de milhares de vezes na viscosidade do fluido (SCHLICHTING; KESTIN,
1979). Outras características marcantes de escoamentos turbulentos são a
multiplicidade de escalas e a aleatoriedade no comportamento das variáveis. Dessa
forma, é intrínseco à sua definição a variação temporal da velocidade e da pressão no
escoamento, por exemplo.
Essas características tornam a análise teórica completa do problema
extremamente complexa, e, até o momento impossível (SCHLICHTING; KESTIN,
1979; FERZIGER; PERIC, 2002). Por outro lado, de um ponto de vista estatístico,
esse desenvolvimento instantâneo das variáveis pode ser representado por uma
velocidade média no tempo (𝑢) e pelas flutuações pontuais (𝑢𝑖′), conforme equação
(67). Isso faz com que os problemas passem a serem abordados segundo seu
comportamento médio temporal, ao passo que a aplicação do TTR é feita sob as
variáveis no formato da equação (67).
𝑢𝑖 = 𝑢 + 𝑢𝑖′ (67)
As equações da continuidade (65) e de NavierStokes (66) são reescritas.
Considerando a Eq.(66), e utilizando a hipótese da viscosidade turbulenta, a mesma
assume a forma apresentada na eq. (68), sendo essas chamadas de equações
médias de Reynolds (ou do inglês RANS, ReynoldsAveraged NavierStokes
equations).
𝜌 [𝜕𝑢𝜕𝑡
+ 𝑢𝜕𝑢𝜕𝑥𝑗
] = 𝜌𝑔𝑖 −𝜕𝑃𝑒𝑓𝑓
𝜕𝑥𝑖+
𝜕
𝜕𝑥𝑗((𝜇 + 𝜇𝑡)
𝜕𝑢𝜕𝑥𝑗
) (68)
Onde 𝑃𝑒𝑓𝑓 é a pressão efetiva, a qual também contempla a energia cinética
turbulenta e μt é a chamada viscosidade turbulenta. Esses termos têm por origem a
53
média do produto das flutuações, no entanto, por mais que provenha do teor
convectivo seu efeito viscoso é predominante nas pequenas escalas do escoamento.
Tais propriedades dependem não só do fluido, como também do escoamento, fazendo
com que novas variáveis sejam adicionadas à solução do problema. Neste sentido, os
modelos de turbulência baseados na hipótese de viscosidade turbulenta se
diferenciam na forma com que as modelam e as resolvem.
O modelo kε (LAUNDER; SHARMA, 1974) propõe o fechamento da
turbulência através de duas equações: a primeira, relacionada a energia cinética
turbulenta (k), e a segunda, à sua taxa de dissipação (ε). É útil em prever escoamentos
viscosos, porém pela carência de modelagem, é impreciso quando há gradientes de
pressão adverso e descolamento da camada limite (ARGYROPOULOS; MARKATOS,
2015).
Em contrapartida, o modelo kω proposto inicialmente por Kolmogorov (1941)
e posteriormente aprimorado por Wilcox (1988) tem como principal característica a
modelagem da camada limite e dos gradientes de pressão – e por isso é considerado
mais indicado do que o modelo kε para escoamentos confinados ou em regiões
parietais; entretanto, carece quando há necessidade de modelagem de correntes
livres (ARGYROPOULOS; MARKATOS, 2015). Também é formulado por duas
equações: referentes à energia cinética turbulenta e à sua taxa de dissipação
específica (ω).
O modelo kω SST proposto por Menter (1994), de forma genérica, separa o
escoamento entre a camada limite e a corrente livre: para o primeiro, ou seja, próximo
de superfícies sólidas o modelo SST busca a formulação do modelo kω, enquanto
que na corrente livre transiciona para o modelo kε. Este modelo sofreu várias
correções ao longo dos anos e atualmente é comumente utilizado até mesmo em
simulações de turbinas, como feito por Wang; Yin; Yan (2019), Riglin et al., (2016),
Schleicher; Riglin; Oztekin (2015), Nuernberg; Tao (2018), Abuan; Howell (2019) e
Hoghooghi; Durali; Kashef (2018).
Durante o desenvolvimento desse trabalho o modelo kω SST é utilizado nas
simulações numéricas via Dinâmica dos Fluidos Computacional. No entanto, ressalta
se aqui que não faz parte do escopo do presente trabalho uma revisão detalhada a
respeito de modelagem da turbulência, sendo o leitor interessando direcionado à
pesquisa realizada por Asnaghi; Svennberg; Bensow (2019). Da mesma forma, o
54
trabalho desenvolvido por Javadi; Nilsson (2017) exemplifica a aplicação do modelo
kω SST com diferentes correções, para a simulação de uma turbina Kaplan.
2.2.2.2 Métodos numéricos
A solução das equações que governam o escoamento começa por sua
modelagem matemática, a qual deverá ser capaz de representar os principais
fenômenos físicos do escoamento. A escolha de cada simplificação possível regime
permanente ou transiente, bidimensional ou tridimensional, turbulento ou laminar,
compressível ou incompressível deve ser sempre ponderada entre simplificação
matemática e representação adequada do problema.
Existem diversos métodos numéricos que podem ser utilizados para encontrar
uma solução aproximada para o modelo matemático proposto. No entanto, ao se
considerar problemas envolvendo escoamento de fluidos em domínios com
geometrias complexas o método mais utilizado é o de Volumes Finitos. Este método
busca discretizar o domínio, dividindoo em elementos infinitesimais. Considerando a
aplicação de malhas nãoestruturadas cada uma das variáveis envolvidas no
problema (as quais correspondem às equações diferenciais parciais) é tratada no
centro do volume de controle, ou seja, em seu centroide (ou simplesmente nó)
(FERZIGER; PERIC, 2002).
As variáveis de interesse são então aproximadas por expansão de séries de
Taylor, e as derivadas segundo a definição trivial de derivações, são consideradas em
relação às diferenças finitas entre os pontos, tornando a solução passível de ser feita
por sistemas lineares. Ademais, em cada nó, essa diferença pode ser interpolada com
referência aos valores anteriores (do inglês Backward Difference Scheme, BDS),
futuros (do inglês Forward Difference Scheme, FDS) ou ambos (do inglês Central
Difference Scheme, CDS), a partir das condições iniciais do problema – seja forças
externas e/ou valores de contorno para as variáveis.
O método upwind altera a interpolação segundo o escoamento: se este flui
positivamente, usase o BDS; caso contrário, usase o FDS. Isso faz com que a
solução seja altamente difusiva (FERZIGER; PERIC, 2002) e esse esquema seja
normalmente utilizado apenas para obtenção de um campo inicial. A acurácia da
solução depende da ordem de expansão da série de Taylor, segundo derivadas de 1ª
55
ordem ou 2ª ordem (não é comum a utilização de ordem superior a três em conjunto
com o método dos Volumes Finitos aplicado à geometrias complexas).
Todo este processo é repetido inúmeras vezes até que os parâmetros em
cada nó, com relação à iteração anterior, atinjam um erro menor que a tolerância
especificada. Para maiores detalhes quanto aos algoritmos e esquemas de resolução,
recomendase os dados disponibilizados por OpenFOAM, 2019.
2.3 OTIMIZAÇÃO
O objetivo das técnicas de otimização é encontrar uma solução que melhor
satisfaça problemas cujo intuito seja encontrar seu valor máximo ou mínimo, sujeito à
determinado limite para as variáveis. Para que isso seja possível, é necessário de
antemão modelar tal problema: primeiramente, as variáveis de entrada que o
modelam; posteriormente o que se deseja otimizar, e por fim, a região confinada do
domínio das variáveis que se pretende obter a resolução (NAGY, 2017).
As variáveis de entrada são intrínsecas à cada problema. Podem ser
caracterizadas tanto como constantes, como em variáveis de domínio das funções.
Ou seja, são elas as responsáveis pela variação da imagem da otimização.
Definese o que se pretende otimizar como função objetivo, seja ela uma
maximização ou minimização. Essa pode ser única (SO, do inglês Single Objective),
que demanda menor custo computacional e menor complexidade em resolvêla; ou
múltipla (MO, do inglês MultiObjective), cuja abordagem é peculiar: as funções
costumam competir umas contra as outras, e se espera uma gama de resultados, uma
vez que a modelagem das consequências de cada uma pode alterar o ponto ótimo
(NAGY, 2017).
Tendo estabelecido tanto as variáveis quanto a(s) função(ões) objetivo(s),
delimitase o espaço a ter a solução buscada: é passível de se aplicar limites tanto na
entrada quanto na resposta desejada, por equações ou inequações. Essas definições
podem ser resumidas pelas equações (69) e (70):
min 𝑥∈ℝ
𝑓𝑖(𝑥), (𝑖 = 1,2,3, . . , 𝑀) (69)
56
𝑆𝑢𝑗𝑒𝑖𝑡𝑜 𝑎 ℎ𝑗 = 0, (𝑗 = 1,2,3, … ,𝑁)
𝑔𝑘 ≤ 0, (𝑘 = 1,2,3, … , 𝑃) (70)
A estratégia de otimização é separada em dois campos: a abordagem
determinística e a heurística. A primeira é caracterizada por uma solução exata, em
muitas vezes por meio de técnicas de cálculo para solução da função objetivo ou para
seu ponto crítico. Um exemplo desse tipo de categoria é o método de gradiente
descendente, utilizado em redes neurais para o processo de aprendizagem. Quanto à
segunda, é sinônimo de uma solução aproximada (segundo uma tolerância), porém
abrangente e irrestrito à complexidade dos problemas. Essa é alcançada de forma
aleatória, com métodos de busca no espaço do domínio e avaliação do resultado.
Essa é a categoria dos métodos a serem detalhados nas próximas seções:
PSO (do inglês, Particle Swarm Optimization), ABC (do inglês, Artificial Bee Colony),
FPA (do inglês, Flower Pollination Algorithm) e SA (do inglês, Simulated Annealing).
Além desses, devese destacar outros métodos recém aplicados em pesquisas de
engenharia. Fetanat; Khorasaninejad (2015) aplicou o método baseado em
colonização de formigas (ACO, do inglês Ant Colony Optimization) com intuito de a
partir de uma demanda elétrica, encontrar a configuração de menor custo em uma
geração híbrida entre eólica e solar, variando o número de turbinas, placas e baterias;
Kamjoo et al. (2016) otimizou o custo relacionado a configuração de uma produção
híbrida de energias renováveis, baseado na teoria genética de Darwin.
2.3.1 Particle Swarm Optimization
O método PSO foi desenvolvido por Kennedy; Eberhart (1995), em que sua
inspiração teórica é o movimento de busca de alimento dos pássaros. Partese da
ideia de que cada pássaro em um bando procura individualmente o local de maior
abundância de alimento e comunica os outros para que, se for a melhor fonte, todos
os outros o seguirem.
No contexto computacional, o alimento é sinônimo da solução ótima do
problema. Cada pássaro é tratado como uma partícula, que possui vetores posição e
velocidade; esses são atualizados a cada iteração, sujeitos à componentes cognitivas
(máximo individual) e sociais (máximo global, ou seja, do grupo de partículas). A
57
representação matemática das equações que governam esse método pode ser vista
nas equações (71) e (72) abaixo.
X 𝑖𝑘+1
= X 𝑖𝑘+ V 𝑖
𝑘+1 (71)
V 𝑖𝑘+1
= 𝑤V 𝑖𝑘+ 𝑐1𝑟𝑖1
𝑘 (C 𝑖 − X 𝑖𝑘) + 𝑐2𝑟𝑖2
𝑘 (C 𝑔 − X 𝑖𝑘) (72)
Onde X e V são os vetores posição e velocidade com dimensão do número
de variáveis de entrada, de cada partícula i em determinada iteração k, w o termo de
inércia, c1 e c2 as constantes cognitivas e sociais, r1 e r2 valores aleatórios entre 0 e 1
e C o máximo valor da função objetivo da partícula i e global g. Pela sua baixa
complexidade e baixo custo computacional, tem sido ferramenta de otimização em
problemas de engenharia: por exemplo, Borhanazad et al. (2014) adotou como
entrada a performance de geração solar, eólica e diesel a fim de otimizar o custo e
tamanho de uma micro rede de geração híbrida de energia elétrica; já Djavareshkian;
Esmaeili (2014) simulou numericamente três perfis NACA com corda e espessura
diferentes e utilizou os dados de arrasto, sustentação e coeficiente de pressão a fim
de desenvolver um hidrofólio que ofereça desempenho ótimo.
2.3.2 Artificial Bee Colony
De maneira análoga ao PSO, o método ABC desenvolvido por Karaboga
(2005) é baseado na busca de néctar por enxame de abelhas. A estratégia consiste
em separálas em três grupos: as abelhas operárias (do inglês employed), as
espectadoras (do inglês onlooker) e as observadoras (do inglês scout). As
observadoras vasculham possíveis fontes de néctar ao redor da colméia, enquanto
que as espectadoras enviam as operárias para que as investiguem; ao retornarem,
comunicam às espectadoras sobre a localização e qualidade dessa fonte (através da
linguagem da dança), e essas avaliam a possibilidade de abandonar ou adotar
determinada fonte como a principal.
No cenário computacional, a fonte de alimento é uma possível solução para o
problema, enquanto que a probabilidade de quantidade de néctar e a distância dessa
para a colméia traduzem a sua qualidade. Normalmente, a quantidade de abelhas
operárias é igual ao número de abelhas espectadoras, e a cada iteração a fonte
58
principal pode ser atualizada. Do ponto de vista matemático, essa estratégia pode ser
descrita segundo as equações (73), (74) e (75).
V 𝑖 = X 𝑖 + 𝑟1𝑖(X 𝑖 − X 𝑘) (73)
X 𝑖 = X 𝑚𝑖𝑛 + 𝑟2(X 𝑚𝑎𝑥 − X 𝑚𝑖𝑛) (74)
𝑝𝑖 = 𝐶𝑖
∑ 𝐶𝑛𝑁𝑛=1
(75)
Onde X representa o vetor de uma das soluções do problema i, de dimensão
dependente da quantidade de variáveis, X k um vetor solução aleatório, V é o vetor de
procura de novas fontes de alimento pelas observadoras, r1i um valor aleatório entre
1 e 1, r2 um valor aleatório entre 0 e 1, Xmax e Xmin os valores limites das possíveis
soluções e pi a função probabilidade da qualidade da fonte, em que C é o valor da
função objetivo particular da solução i e de todas as abelhas n, sujeito à uma
população de N abelhas operárias. Por meio desse método, Yildiz (2013) otimizou
estruturalmente um componente automotivo, e comparou estritamente com outros
métodos.
2.3.3 Flower Pollination Algorithm
Naturalmente, de forma semelhante ao método anterior, o FPA tem como
inspiração a natureza: a reprodução das plantas. Tal algoritmo fora desenvolvido por
Yang (2012), e consiste na sobrevivência do “melhor indivíduo”, ou seja, o pólen,
tratado como uma possível solução do problema, é testado segundo a função objetivo
e o quê apresentar menor/maior valor, sobrevive.
O método é dividido em quatro regras, executadas em cada iteração:
i. na primeira, a busca global ou polinização biótica, que analogicamente
na natureza é o transporte do pólen feito por animais e insetos,
traduzido no formato matemático da equação (76). Têmse como meio
de busca o voo de Lévy, ou a distribuição de Lévy;
59
X 𝑖𝑘+1
= X 𝑖𝑘+ 𝐿 (X 𝑖
𝑘− 𝑔∗) (76)
Onde X representa o vetor de solução, de cada pólen i em determinada
iteração k, L a probabilidade de Lévy e g* a melhor solução encontrada.
ii. na segunda, a busca local ou polinização abiótica, e representa o papel
do vento e da água na dispersão regional do pólen;
iii. a terceira, a manutenção de constância da reprodução, ou seja, que
apenas mesmas espécies de flores se reproduzem entre si, garantindo
a qualidade da solução. Tanto essa, quanto a regra anterior são
sintetizadas segundo equação (77);
X 𝑖𝑘+1
= X 𝑖𝑘+ 𝜖 (X 𝑗
𝑘− X 𝑚
𝑘) (77)
Onde ε representa um número aleatório entre 0 e 1, para que uma distribuição
linear seja feita, e j e m pólens diferentes de i. Ou seja, a solução é trocada quando
as soluções j e m forem originadas da mesma população.
iv. e por fim, a troca entre reprodução local e global a partir de uma
distribuição de probabilidade – parecido com o aplicado por ABC.
Como exemplo, um estudo aplicado de tal método em estruturas, Bekdaş;
Nigdeli; Yang (2015) buscou otimizar o peso de treliças, comparando o desempenho
de alguns métodos com o da sua pesquisa.
2.3.4 Simulated Annealing
Por fim, um dos percussores nos métodos metaheurísticos, o método SA fora
desenvolvido por Kirkpatrick; Gelatt; Vecchi (1983) em inspiração ao arranjo dos
átomos no processo de recozimento: os contornos de grãos buscam os espaços cuja
energia (ou entropia) é a menor.
Neste sentido, partindo de uma solução inicial, seu valor segundo a função
objetivo é classificado como sua energia em um estado. Então, a solução é
60
aleatoriamente e localmente variada um determinado número de tentativas por “ciclo
de resfriamento””, e a energia desta nova solução é comparada com a que se
encontra. A decisão de adotar tal solução, ou seja, a seleção de um bom movimento
ou não, é ou em comparação com a menor energia global, ou em base probabilística
do resultado da Equação (78) ser menor que um número aleatório entre 0 e 1, o
chamado “teste de Metropolis”. Tal processo é iterativo segundo a “diminuição” da
temperatura (Equação (79)), cujo valor é então calculado pelas probabilidades
mínimas e máximas de se aceitar uma solução (Equação (80). Ressaltase que as
equações aqui apresentadas estão modificadas pela melhora de performance
apresentado por Hasançebi; Çarbaş; Saka (2010).
P(T) = ψ ∗ exp (−∆𝐸
𝐾𝑇) (78)
∆𝑇 = (𝑇𝑓𝑖𝑛𝑎𝑙
𝑇𝑖𝑛𝑖𝑐𝑖𝑎𝑙)(
1
𝑁−1)
(79)
𝑇 = −1
ln 𝑝 (80)
Onde P(T) é a probabilidade de aceite da solução, ψ o fator de correção, ΔE
é a diferença de energia entre os estados, K é o parâmetro de Boltzman,T a
temperatura, N o número de “ciclos de resfriamento” e p a probabilidade de aceite da
solução. A título de exemplo de tal aplicação, Zhang et al., (2018) aplicou tal algoritmo
para um dimensionamento economicamente viável de um sistema de geração de
energia residencial por fontes renováveis (eólica e/ou solar) acoplado a um
armazenador de energia elétrica (de forma química, por hidrogênio, ou eletroquímica,
por baterias).
61
3 METODOLOGIA
Tendo como objetivo atender uma demanda elétrica, uma turbina hidrocinética
é então desenvolvida na tentativa de melhor aproveitar o recurso hídrico, ou seja, com
um design otimizado. Para que isso seja possível, devese de antemão replicar o
desempenho de um rotor utilizando as equações apresentadas na Seção 2.2 por meio
de algoritmos, fazendo uso da linguagem de programação Python. Essa reprodução
é considerada fiel a partir da validação com não só dados experimentais, como
também trabalhos anteriores com abordagem semelhante.
Com tal objetivo atingido, o primeiro passo para a otimização é com relação à
sua performance hidrodinâmica. Primeiramente, o espaço reduzido entre dois perfis e
seu efeito nos esforços de arrasto e sustentação são investigados por meio de DFC.
Em turbinas com baixa solidez, este fenômeno pode ser desprezado, uma vez que
tais valores se aproximam dos apresentados em perfis únicos (cujo valor de solidez
tende a zero); entretanto, como visto na Seção 2.2.1.9, se este valor crescer, os efeitos
hidrodinâmicos sob o perfil tendem a se alterar. Segundo a Eq.(51), o aumento deste
número pode ter como origem a diminuição do diâmetro do rotor e um maior número
de pás.
Em seguida, é analisado um artifício de mecânica dos fluidos para estimular
a transição da camada limite sob um corpo: protuberâncias na região de
desenvolvimento do escoamento representam um obstáculo para o fluido,
aumentando sua turbulência e o mantendo próximo ao corpo. Dessa forma, é
adicionado um ressalto na região da borda de ataque ao modelo DFC dos perfis único
e em cascata utilizados anteriormente. O perfil hidrodinâmico escolhido foi o NACA
4412, já utilizado em trabalhos semelhantes, enquanto que o software de código
aberto OpenFoam é adotado para realizar as simulações numéricas.
Além das alterações citadas, fazse uso dos resultados e de um modelo de
difusor apresentado na literatura para construção do algoritmo final. As técnicas de
otimização utilizadas e comparadas no quesito custo computacional e precisão do
resultado são as descritas na Seção 2.3, em que os objetivos adotados são a
maximização do coeficiente de potência e minimização do momento de inércia –
assim, o torque inicial para início da operação da turbina é reduzido. Os valores de
entrada – e consequentemente, os limites de projeto , diferente de alguns trabalhos
com abordagem semelhante e aqui citados, são: o diâmetro do rotor, o número de pás
62
e a rotação das pás. As constantes se reduzem à velocidade do rio, densidade da
água e profundidade de disposição da turbina.
Para isso, partese de uma demanda elétrica mensal de 182 kWh por
consumidor residencial, valor esse previsto para 2026 segundo a EPE; ONS (2017)
(Empresa de Pesquisa Energética e Operador Nacional do Sistema elétrico), a fim de
prédimensionar a área necessária, segundo as teorias de rotor apresentadas. A
velocidade dos rios é intrínseca ao local de operação; entretanto, valores encontrados
na literatura para alguns rios brasileiros são: entre 1,7 e 2,5 m/s para o rio Iguaçu e
quase 1,6 m/s para o rio Paraná, ambos situados no estado do Paraná; em torno de
1,5 m/s para o rio Jamari, localizado no estado de Rondônia, e superior a 1,2 m/s para
o rio CuruáUna, no estado do Pará. Assim, é possível prédimensionar os limites do
diâmetro. Quanto ao restante das variáveis de entrada, seus limites são estabelecidos
por boa prática apresentada na literatura. Por outro lado, as constantes citadas são
dependentes do local de aplicação.
Por fim, com intuito de validar o desempenho apresentado pelo design ótimo,
a turbina é então simulada numericamente por meio do DFC.
Figura 19. Ilustração do fluxograma de atividades.
Fonte: Próprio autor.
63
4 RESULTADOS
A forma com que os resultados são divididos é: na primeira seção, há a
validação dos modelos apresentados no Capítulo 2 segundo resultados
experimentais; tendo replicado o funcionamento da turbina, buscase otimizar sua
potência e a inércia, variando seu tamanho, layout e operação com os modelos de
otimização já apresentados; a partir das configurações obtidas, um modelo
tridimensional é construído na intenção de validar a performance através de DFC. Por
fim, realizase as otimizações estruturais: a correção de solidez, a adição do difusor e
o ressalto na borda de ataque do perfil. Tais etapas podem ser representadas segundo
fluxograma na Figura 19.
4.1 VALIDAÇÃO
4.1.1 MODELOS TEÓRICOS
Os resultados experimentais a serem validados são os obtidos por Maniaci;
Kelley; Chiu, 2015, que buscaram obter os parâmetros aerodinâmicos do
funcionamento de uma turbina eólica denominada G1, cujos testes foram realizados
na Universidade Tecnológica de Munique (TUM). Ela possui rotor de três pás e cubo
com raio de 0,550 m e 0,054 m, respectivamente; fora projetada em perfil
aerodinâmico RG14 para maior parte da pá (nos primeiros 13,3% é composto por
formato cilíndrico, seguido por uma suavização até o perfil por mais quase 14%) e
com λ igual à 6,6 (ou velocidade de rotação igual a 850 rpm). Neste estudo, os dados
utilizados apresentam este valor igual a 7.
O processo de validação adotado consiste em avaliar uma variável por vez,
em que a primeira faz uso de todos os dados obtidos em teste para seu cálculo, e ao
passo que os modelos repliquem o comportamento físico, tornar as representações
posteriores independentes de tais dados experimentais – à exceção da distribuição
da corda e de torção da pá, intrínsecos ao problema. A ordem dos coeficientes
validados é: fator de indução axial (a), fator de correção nas pontas (F), fator de
indução tangencial (a’), coeficiente de potência pontual (CP), ângulo de ataque (α),
64
distribuição do coeficiente de arrasto (CD) e sustentação (CL), e a curva de coeficiente
de potência. A pá é dividida em 100 seções transversais para que cada valor seja
calculado localmente e a apresentação do resultado, com exceção do coeficiente de
potência, é dada em função de sua distribuição ao longo do comprimento da pá.
As Equações utilizadas na verificação do fator de indução axial são as (22),
(28) e (32). Para a primeira, resolvese analiticamente; para a segunda e a terceira,
solucionase numericamente através da ferramenta fsolve, presente na biblioteca
linguagem de programação Python, denominada scipy; como requisito, são
necessários estimativas iniciais aqui iguais a 0,33. Notase que todas os três modelos
dependem do fator de perda nas pontas (F); como a validação nesta etapa indispõe
desses dados, assumese tais valores como 1. O resultado pode ser visto na Figura
20.
Figura 20. Resultado da validação dos modelos para fator de indução axial.
Fonte: Próprio autor.
Destacase a boa previsibilidade dos modelos de Glauert (1935) e Wilson e
Lissaman (1974) ao longo de praticamente toda a pá e com os desvios naturais da
simplificação das perdas nas pontas; em contrapartida, vêse que o modelo proposto
65
por Shen et al., (2005), ao menos neste caso, não oferece previsibilidade suficiente
do comportamento físico do escoamento.
Como os resultados em suma foram aparentemente prejudicados, fazse
necessária na próxima etapa de validação dos modelos de perdas nas pontas da pá
– Equações (18), (24), (34), (44) e (45) , a comparação entre tanto o fator de indução
axial proposto por Glauert (1935), como por Wilson e Lissaman (1974). Ainda, adota
se como correção do elevado fator de indução a proposta de Buhl (2005) (Equação
(38)) e da mesma forma que aplicado por Marten et al., 2010. Tais resultados podem
ser vistos nas Figura 21 e Figura 22.
Figura 21. Resultado da validação dos modelos de perdas na ponta da pá, para o fator de
indução axial de Glauert (1935).
Fonte: Próprio autor.
66
Figura 22. Resultado da validação dos modelos de perdas na ponta da pá, para o fator de
indução axial de Wilson e Lissaman (1974).
Fonte: Próprio autor.
É confirmado pelas figuras acima a hipótese apresentada da origem da falha
em prever o fator de indução axial na ponta da pá morar na redução do fator de perdas
à 1. Entretanto, o modelo previsto por Moriarty e Hansen (2005), que leva em conta
as possíveis perdas decorrentes do cubo, mostrase falho para ambos os modelos de
fator de indução axial. O mesmo podese dizer da combinação do fator de perdas na
ponta da pá previsto por Shen et al., (2005) junto ao fator de indução axial de Glauert
(1935), e do previsto por Burton (2011) com o modelo de Wilson e Lissaman (1974).
Fora as combinações citadas, as demais apresentam resultado satisfatório no
sentido de prever o comportamento físico; destacase a combinação do fator de
indução de Wilson e Lissaman (1974) utilizando as perdas na ponta de Shen et al.
(2005), uma vez que esses valores aparentam menor desvio do experimental por
quase toda a pá contra uma pequena diferença dos modelos apresentados por Glauert
(1935) e Prandtl (1927); além disso, é subestimado na ponta – ao contrário dos demais
resultados do fator de indução de Glauert (1935).
Devese ressaltar que a provável justificativa de tal resultado, mora na
complexidade dos modelos: a equação prevista por Wilson e Lissaman (1974), como
67
já apresentado no Capítulo 2, aprimorou a indução axial ao considerar o escoamento
no plano da esteira e a força de sustentação decorrente – força motriz do movimento
rotativo da pá; do mesmo modo, o fator de perda na ponta das pás previsto por Shen
et al., (2005) passa a incluir a correção física de inconsistência das forças atuantes
neste local.
Já quanto ao fator de indução tangencial, por carência de resultados
experimentais diretos, sua validação é feita de forma indireta através do ângulo de
escoamento (Equação (16)): passase então a combinar os modelos de indução axial
e perdas nas pontas da pá previamente analisados com às possíveis equações do
fator de indução tangencial ((23), (33), (37) e (39)). Como o ângulo é dependente de
ambos os fatores, um processo iterativo é feito para que se alcance a convergência
do resultado para uma tolerância predeterminada – no caso, 1E6. Tal erro é
mensurado pela diferença entre o valor presente e o anterior do fator de indução
tangencial, já que este é o parâmetro a ser validado. O resultado pode ser visto na
Figura 23.
Figura 23. Resultado da validação dos modelos de fator de indução tangencial.
Fonte: Próprio autor.
68
A figura acima evidencia a evolução dos modelos teóricos em prever o
desempenho fluidodinâmico das turbinas: notase como os mais recém desenvolvidos
preveem quase que por completo o ângulo de escoamento. Com exceção do proposto
por Glauert (1935), que diverge no que tange a ponta da pá, a combinação do já
validado fator de indução axial com os modelos restantes de fator de indução
tangencial pode ser aplicada na replicação da performance dos parâmetros citados.
Por conveniência de adoção de modelo completo, tomase o desenvolvido por Vries
(1979) para as validações a seguir. Sob os parâmetros já analisados, é possível então
a verificação do coeficiente de potência no ponto experimental. Nesta etapa, são
comparadas as Equações (20), (46), (47) e (48), além da utilizada em estudo
semelhante por Brasil Junior et al., (2019) – Equação (81). O erro passa a ser a
diferença entre os valores passado e presente do ângulo de escoamento, e a integral
fora calculada numericamente pelo Método Trapezoidal, presente na biblioteca numpy
da linguagem de programação Python. Os resultados podem ser vistos na Tabela 1,
cuja diferença percentual é com relação ao experimental.
𝐶𝑃 = 8
𝜆2∫𝑎′𝜆
0
(1 − 𝑎)(1 −𝐶𝐷𝐶𝐿cotΦ)𝜆𝑟
3𝑑𝜆𝑟 (81)
Tabela 1 Validação dos modelos de coeficiente de potência no ponto experimental.
Proposta Resultado Diferença (%)
Experimental 0,4204
Glauert (1935) 0,3582 14,80
Jonkman (2003) 0,1384 67,08
Letcher (2017) 0,2004 52,33
Branlard (2017) 0,5123 21,86
Brasil Junior et al., (2019) 0,4145 1,40 Fonte: Próprio autor.
Uma das dificuldades encontradas em prever aquele que talvez seja o
coeficiente mais importante em respeito ao desempenho de turbinas, o de potência,
consiste na interpretação e replicação das equações. O presente estudo seguiu o
procedimento dos respectivos autores, e tal obstáculo – bem como a ineficácia do
método pode ser a origem da alta disparidade dos resultados das propostas de
Jonkman (2003) e Letcher (2017).
69
Interessante notar que, diferente das propostas anteriores, a previsão de
Glauert (1935) para o coeficiente em questão oferece resultado considerável e não
será descartado em análises futuras. Dessa forma, juntase à proposta de Brasil
Junior et al., (2019), cujo resultado impressiona pela quase exatidão. Segundo o autor,
seu método leva em conta as perdas decorrentes do escoamento em torno da pá, ao
passo que considera seu ângulo de incidência e os coeficientes de arrasto e
sustentação.
Muito além de um único ponto de coeficiente de potência, um modelo de boa
qualidade é capaz de prever essa performance ao longo de inúmeros regimes de
operação, traduzido pela razão entre a velocidade tangencial da pá e a velocidade de
corrente livre (λ). Essa variedade pode acarretar em mudanças drásticas no
escoamento, como o fenômeno de estol. O fator crucial que prevê tal condição
fluidodinâmica, como também o regime pleno, é o ângulo de ataque (Equação (15)).
Sua validação pode ser vista na Figura 24, e fora realizada de forma análoga ao
ângulo de escoamento.
Figura 24. Resultado da validação da previsão do ângulo de ataque.
Fonte: Próprio autor.
70
Uma boa previsão do ângulo de ataque, como apresentado acima, tem como
consequência a precisa determinação dos coeficientes de arrasto e sustentação sob
as diferentes condições, como já citado. Dessa forma, o presente trabalho procurou
replicar o processo utilizado em pesquisa semelhante por Martínez et al., (2005), e os
resultados podem ser vistos nas Figura 25 e Figura 26 .
Figura 25. Resultado da validação de previsão do coeficiente de arrasto.
Fonte: Próprio autor.
71
Figura 26. Resultado da validação de previsão do coeficiente de sustentação.
Fonte: Próprio autor.
Como apresentado em ambas as figuras, de uma forma geral, ou seja, em
maior parte da pá, os coeficientes podem ser previstos adequadamente. O desvio na
ponta da pá do coeficiente de sustentação, bem como em quase a totalidade do
arrasto, em relação ao experimental é na escala milesimal. Por outro lado, devese
ressaltar a carência do modelo em prever o coeficiente de sustentação perto do cubo.
Tal localização coincide não só com ângulos de ataque maiores (vide Figura 24) que
pode ser uma evidencia da dificuldade em prever, ao menos analiticamente, esse
coeficiente em regime de estol; como também é o intervalo no comprimento da pá em
que a seção transversal não é o perfil aerodinâmico em estudo – mas sim, o formato
cilíndrico e a suavização até o perfil.
Assumindo que tanto os fatores intrínsecos ao funcionamento da turbina como
os parâmetros de escoamento são previstos de forma considerável, fazse possível a
última etapa da validação – e talvez mais importante: a curva de coeficiente de
potência. Como ressaltado anteriormente, serão comparados os coeficientes de
potência de Glauert (1935) e Brasil Junior et al., (2019), e também os dados
apresentados por software desenvolvido pela Universidade Tecnológica de Berlim (do
alemão TU Berlin), de código aberto e amplamente utilizado em casos como esse,
72
denominado QBlade (MARTEN et al., 2013). Os critérios de erro e processo
construtivo seguiram a validação do coeficiente de potência pontual, com a única
diferença na tolerância do erro: por questões de convergência, tal valor passou a ser
de 3E2. A curva pode ser vista na Figura 27.
Figura 27. Resultado da validação da curva de potência.
Fonte: Próprio autor.
Prever o coeficiente de potência sob condições externas auxilia diretamente
em projetos de aplicação das turbinas, sejam elas eólicas ou hidrocinéticas. Através
da figura acima, tal reflexão é explícita: na hipótese do projetista se basear no modelo
de Glauert (1935), a máquina pode ser tanto subdimensionada (por exemplo, quando
λ está entre 3 e 5,5, ou 8 e 9), quanto superdimensionada (para λ está entre 5,5 e 8).
Ainda, devese destacar o comportamento deste valor. Enquanto o outro modelo e os
resultados do QBlade seguem o desempenho experimental, no que tange à um
crescimento limitado seguido de uma queda posterior, os valores de Glauert (1935)
não decrescem e ainda tendem a aumentar, como pode ser visto quando λ é igual a
9.
Ademais, é enfatizado como a construção segundo inúmeras propostas da
literatura de parâmetro por parâmetro do funcionamento foi fundamental à previsão: a
73
curva gerada pelo roteiro obteve precisão considerável, superando o resultado de um
software consolidado neste tipo de estudo.
Na tentativa de validar a curva de coeficiente de empuxo através de roteiro
semelhante ao da curva de coeficiente de potência, os resultados não foram
satisfatórios, seja por imprecisão dos modelos aqui apresentados, seja por
dificuldades semelhantes às apresentadas pela implementação do coeficiente de
potência. Assim, ao invés de integrarse os valores locais ao longo da pá, buscouse
implementar um valor médio deles, segundo Eq. (82). São comparadas as Equações
(21), (31), (36), (49) e (50), e o resultado pode ser visto na Figura 28.
𝐶𝑇,𝑚é𝑑𝑖𝑜 = ∫ 𝐶𝑇𝑅
𝐻(𝑟)𝑑𝑟
∫ 𝑑𝑟𝑅
𝐻
(82)
Figura 28. Resultado da validação dos valores de empuxo.
Fonte: Próprio autor.
Em contraste com as correlações de coeficiente de potência, todos os
modelos aplicados replicaram o coeficiente de empuxo de forma satisfatória.
Especialmente quando se trata de valores intermediários de λ, cujos parâmetros de
todos os modelos praticamente coincidiram com os dados experimentais. Destacase
74
que, mais uma vez, os regimes de altos valores de λ tendem a ser delicados quando
o intuito é predizer o coeficiente de empuxo: enquanto Spera (1994) os
superdimensionam, os outros os subdimensionam. Embora em geral o modelo
proposto por Manwell et. al (2011) tenha apresentado o melhor desempenho, fora
adotado neste trabalho o modelo de Buhl (2005), já que além de sua performance ter
sido próxima, este é utilizado em outro software de mesmo intuito, desenvolvido pela
NREL (do inglês, National Renewable Energy Laboratory) denominado AeroDyn
(PRATUMNOPHARAT; LEUNG, 2011).
Um resumo dos modelos teóricos implementados e validados está contido na
Tabela 2.
Tabela 2 Resumo dos modelos validados em replicar funcionamento da turbina.
Variável Modeloteórico Equação
Fator de indução axial (a) Wilson e Lissaman (1974) (28)
Fator de perdas nas pontas da pá (F) Shen et al., (2005) (34) e (35)
Fator de indução tangencial (a’) Vries (1979) (37)
Coeficientes de arrasto e sustentação (CD e CL) Martínez et al., (2005)
Coeficiente de potência (CP) Brasil Junior et al., (2019) (81)
Coeficiente de empuxo (CT) Buhl (2005) (50)
Correção Buhl (2005) (38) Fonte: Próprio autor.
Assim sendo, assumese que o algoritmo de validação implementado replica
o funcionamento da turbina em diferentes características operacionais.
4.1.2 PARÂMETROS DE SIMULAÇÃO DFC
Os coeficientes fluidodinâmicos validados são parte de experimento de
Wadcock (1987), realizado junto à NASA em túnel de vento com o perfil NACA 4412.
O teste fora feito sob as condições de número de Reynolds igual a 1,64E+6 e
intensidade de energia cinética turbulenta igual a 0,1% e comprimento da frequência
de mistura turbulenta de 1 m. O perfil possui corda de 0,9 m e os ângulos de ataque
variaram de 0º à 18º. A validação consiste em simular computacionalmente o mesmo
perfil e condições de teste, a fim de encontrar uma distribuição e quantidade de
75
elementos de malha numérica tal qual os coeficientes de arrasto (CD) e sustentação
(CL) sejam suficientemente replicados. Essa etapa é denominada independência de
malha, e dela fazse possível obter os parâmetros em análise para um domínio
computacional com elementos infinitesimais, através do método GCI (do inglês Grid
Convergence Index). As Equações (83), (84) e (85) abaixo evidenciam o
procedimento, que consiste em extrapolar os dados a partir dos resultados de uma
malha grosseira, outra intermediária e outra refinada (ROACHE, 1998).
𝜑𝐺𝐶𝐼 ≅ 𝜑𝐴𝜏𝐵𝐴
𝑞 − 𝜑𝐵𝜏𝐵𝐴𝑞 − 1
(83)
𝑞 = 1
ln(𝜏𝐵𝐴)|ln |
𝛿𝐶𝐵𝛿𝐵𝐴
| + ln(𝜏𝐵𝐴
𝑞 − 1. 𝑠𝑔𝑛 (𝛿𝐶𝐵
𝛿𝐵𝐴)
𝜏𝐶𝐵𝑞 − 1. 𝑠𝑔𝑛 (𝛿𝐶𝐵
𝛿𝐵𝐴))| (84)
𝛿𝐵𝐴 = 𝜑𝐴 − 𝜑𝐵𝜑𝐴
(85)
Em que φ é o coeficiente em análise, τ a razão entre as malhas, q a ordem de
convergência e δ o erro relativo dos coeficientes em análise. Para isso, as três malhas
Cgrid bidimensionais são construídas por meio do software ICEM, com o aerofólio de
corda 1 m localizado no centro de um raio igual a 50 m. A elaboração da malha tem
por base o valor de y+ igual a 3,5 na parede, o qual resulta em uma distância de 5,0E
5 m para o comprimento normal do primeiro elemento junto à parede do perfil. De
forma genérica, os contornos do domínio computacional são: a entrada e saída do
fluido (Inlet e Outlet, respectivamente), o próprio perfil (aqui chamado de NACA_TOP)
e a periodicidade nas porções inferior e superior. A denominação, a quantidade de
elementos e a qualidade mínima – esse segundo critério de Determinante 3x3x3 – de
cada uma das malhas podem ser vistas na Tabela 3.
Tabela 3 Detalhes das malhas construídas para a validação do perfil isolado.
Denominação Número de elementos Qualidade mínima Razão entre
malhas (τ)
Malha A 292 000 0,7290 1,28
Malha B 375 000 0,8300 1,28
Malha C 480 000 0,8650 Fonte: Próprio autor.
76
O software empregado para as simulações computacionais (processamento
dos dados) foi o OpenFOAM (OPENFOAM, 2019). A partir do valor do número de
Reynolds e do comprimento característico – ou seja, a corda do perfil , e com intuito
de trazer as simulações para o padrão deste trabalho, tomouse por base as
propriedades da água para se estimar a velocidade do escoamento. Assim, sob uma
densidade de 998,2 kg/m3 e uma viscosidade de 1,003E3 kg/m*s, esse valor foi de
1,64788 m/s. Como se trata de um escoamento turbulento, o modelo utilizado foi o k
ω SST. Esse modelo foi escolhido por se tratar de um modelo já consolidado para
esse tipo de simulação, tendo sido aplicado com sucesso por outros estudos
(CEBRIÁN; ORTEGACASANOVA; FERNANDEZFERIA, 2013; SCHLEICHER;
RIGLIN; OZTEKIN, 2015; TIAN et al., 2016; HALDER; RHEE; SAMAD, 2017; JAVADI;
NILSSON, 2017). Um resumo das condições iniciais pode ser visto na Tabela 4.
Figura 29. Imagens representativas das malhas criadas e do zoom em torno do perfil NACA.
Fonte: Próprio autor.
77
Tabela 4 Condições iniciais para as simulações de perfil isolado. Região U p k omega nut
internalField (0,0 0,0 0,0) 0,0 uniform 0,001 uniform 10,0 uniform 1e5
Inlet uniformFixedValue constant (UxUy 0,0)
zeroGradient
turbulentIntensityKineticEnergyInlet
uniform 1,0 intensity 0,001
turbulentMixingLengthFrequencyInlet uniform 10,0
mixingLength 1,0
calculated uniform
0,0
Outlet zeroGradient fixedValue uniform 0,0
turbulentIntensityKineticEnergyInlet
uniform 1,0 intensity 0,001
turbulentMixingLengthFrequencyInlet uniform 10,0
mixingLength 1,0
calculated uniform
0,0
NACA_TOP fixedValue
uniform (0,0 0,0 0,0) zeroGradient
kqRWallFunction uniform 1e10
omegaWallFunction value 10,0
nutUWallFunction uniform
0,0
Periodic/Shadow Cyclic cyclic cyclic cyclic cyclic Fonte: Próprio autor.
Onde U é a velocidade, p a pressão, k a energia cinética turbulenta, omega a
taxa de dissipação específica de turbulência e nut a viscosidade cinemática turbulenta.
A construção do desenvolvimento dos coeficientes passa pela análise dos ângulos de
ataque, que foram: 0º, 5º, 10º, 12º, 15º e 18º. A mudança do ângulo de ataque se dá
pelas componentes de velocidade na entrada, segundo Equação (86).
𝑈𝑥 = 𝑈0 cos 𝛼𝑈𝑦 = 𝑈0 sin 𝛼
(86)
Pela complexidade do escoamento, optouse por um regime transiente. Antes,
porém, contouse com uma simulação inicial para 0° em regime permanente (solver
SIMPLE), resíduo igual a 1E6 e equações resolvidas em 1ª ordem (ou upwind) para
obtenção de um campo inicial; deste campo, alterase para solução de 2ª ordem (ou
linearUpwind) já em regime transiente (solver PISO), para resíduo igual a 1E5. O
passo de tempo fora de 1E3, e os fatores de subrelaxação para o regime transiente
são: 0,1 para pressão e 0,5 para velocidade, energia cinética turbulenta e taxa de
dissipação específica de turbulência. Exclusivamente para 0°, é simulado apenas 10
segundos físicos e posteriormente outro mesmo intervalo de tempo para que se
obtenha os resultados médios; para os demais, como é tomado a simulação de ângulo
78
de ataque anterior como campo inicial, são simulados 70 segundos físicos e da
mesma forma novamente outros 10 segundos para se obter os resultados médios.
São exportados os coeficientes de arrasto e sustentação, e conforme o vetor
direção varia segundo o ângulo de ataque, ambos são determinados segundo
Equação (87). Os resultados da independência de malha podem ser vistos nas Figura
30 e Figura 31, em que conta também com os valores apresentados pelo software de
código aberto desenvolvido pelo Instituto Tecnológico de Massachusetts (do inglês
MIT) denominado XFOIL (DRELA, 2001).
𝐶𝐿 = (−sin𝛼 , cos 𝛼 , 0)𝐶𝐷 = (cos𝛼 , sin 𝛼 , 0)
(87)
Figura 30. Validação do coeficiente de arrasto.
Fonte: Próprio autor.
79
Figura 31. Validação do coeficiente de sustentação.
Fonte: Próprio autor.
É interessante separar a análise de ambas as Figura 30 e Figura 31 em dois
regimes: o plenamente estabelecido e o de estol. A diferenciação desses regimes por
parte dos resultados das simulações numéricas coincide com o experimental, ou seja,
até um ângulo de ataque igual a 12º.
No que diz respeito à primeira região, vêse que os resultados das três malhas
em questão, tanto para o coeficiente de arrasto quanto para o de sustentação, são
próximos entre si até o ângulo de 10º. Dali em diante, naturalmente pela menor
quantidade de elementos, a Malha A começa a destoar de forma significativa dos
resultados das outras Malhas B e C. Estes resultados destas malhas, inclusive, valem
a nota de quão adjacentes são entre si para ambos os coeficientes, traduzido pelo
valor contíguo de GCI. Embora os coeficientes de sustentação numéricos
apresentaram diferença sistêmica com relação aos dados experimentais, destacase
os coeficientes de arrasto alcançados de forma quase exata. Diferentemente do
primeiro, em que o XFOIL forneceu um comportamento preciso, os parâmetros de
arrasto acompanharam os números obtidos pelo XFOIL até o ângulo de 10º, e depois
ainda o superaram quando se chega à 12º.
80
A segunda região, o estol, ocorre conforme o ângulo de ataque aumenta,
fazendo com que o ponto de descolamento da camada limite se desloque em direção
à borda de ataque, e a esteira viscosa cresça. Consequentemente, os efeitos
tridimensionais do escoamento também são amplificados. Por isso, qualquer
simulação numérica bidimensional de aero/hidrofólios neste regime tende a
apresentar números subdimensionados com relação ao experimental, como já
previsto por Matyushenko; Kotov; Garbaruk (2017) e confirmado pelas Figura 30 e
Figura 31. Apesar disso, mesmo a malha mais robusta (Malha A) fora capaz de prever
mesmo que imprecisamente o comportamento das curvas dos coeficientes de
maneira mais razoável que o XFOIL. No que diz respeito às malhas intermediária e
fina (B e C, respectivamente), essas apresentaram mais uma vez parâmetros
próximos entre si, o que novamente pode ser traduzido pela curva contínua do GCI.
Diante do comportamento dessas curvas e do GCI (a de arrasto e a de
sustentação) e dos próprios coeficientes obtidos em geral, podese afirmar que a
independência de malha e a validação dos parâmetros com relação aos dados
experimentais foram atingidas pela Malha B. Dessa forma, as distribuições dos
elementos aplicados nessa malha, bem como os valores de y+ dos elementos, serão
utilizados nos estudos que envolverem a simulação numérica computacional.
4.2 OTIMIZAÇÃO
4.2.1 DESIGN
4.2.1.1 Algoritmo
Assumindo o fato de que a replicação do funcionamento de uma turbina fora
alcançada por meio de implementação dos modelos teóricos em um algoritmo, pode
ser feita então a otimização de seu design com intuito de maximizar a potência e
minimizar sua massa. Esta se dá pela utilização dos métodos heurísticos descritos no
Capítulo 2, de forma multiobjetiva e restritiva a uma região no domínio de entrada.
81
É avaliado então a geração local e isolada de energia elétrica fazendo uso de
turbina hidrocinética, cujas condições ambientais são hipotéticas: temperatura
ambiente de 25 °C, velocidade de corrente igual à 1,64788 m/s e profundidade de 1
m. Tal velocidade fora assumida para a manutenção do número de Reynolds validado
de 1,64E6. Neste sentido, adotouse o perfil NACA 4412 por demonstrar, para este
Reynolds, maior razão entre sustentação e arrasto.
Para tal, fezse uso do XFOIL (DRELA, 2001), em que o resultado da razão é
de aproximadamente 148 e coeficiente de pressão mínimo de 1,656, para um ângulo
de ataque igual a 5,4° (números consideráveis frente ao experimental desenvolvido
por Wadcock, (1987)). A corda e a torção, variáveis intrínsecas ao problema, são
determinadas segundo as Equações (59) e (15), com um fator de segurança de 5%.
Um resumo dessas constantes e outras demais está presente na Tabela 5.
Tabela 5 Constantes de projeto e operação.
Constante Valor
Densidade (kg/m3) 997
Velocidade de corrente livre (m/s) 1,64788
Pressão atmosférica (Pa) 1E5
Pressão de vapor (Pa) 3,17E3
Gravidade (m/s2) 9,81
Profundidade (m) 1,0
Coeficiente de arrasto 0,00723
Coeficiente de sustentação 1,0707
Ângulo de ataque (°) 5,4
Coeficiente de pressão mínimo 1,656
Fator de segurança (%) 5 Fonte: Próprio autor.
As variáveis de layout a serem analisadas sob uma região do domínio
delimitado serão o diâmetro, a velocidade rotativa e o número de pás. Partindo da
demanda mensal por consumidor residencial, previsto por EPE; ONS (2017) para
2026, em 182 kWh, e assumindo uma geração elétrica constante, ou seja, em que a
velocidade de corrente livre seja a mesma ao longo do dia e do mês, a potência
necessária é de 252,78 W.
Assim, através da Equação (6) definese os valores mínimo e máximo de
diâmetro: respectivamente, apenas para tal delimitação, assumese um CP máximo
82
de Betz (0,5793) e um mínimo de 0,2. Dessa forma é possível determinar as rotações
mínima e máxima, conforme a razão de velocidade tangencial na ponta da pá com a
velocidade de corrente livre (λ) iguais a 3 e 10. Intervalo esse frequentemente
apresentado na literatura, tal como a recomendação do número de pás: 3 e 12. Os
valores limites das variáveis podem ser visto na Tabela 6.
Tabela 6 Delimitação da região de domínio.
Variável Valor mínimo Valor máximo
Coeficiente de potência 0,2000 0,5793
Diâmetro (m) 0,4991 0,8493
Número de pás 3 12
Velocidade rotativa (rpm) 189,1884 370,5412 Fonte: Próprio autor.
Devese salientar que por mais que os coeficientes hidrodinâmicos do
escoamento não apresentem diferença significativa com a mudança do número de
Reynolds decorrente da alteração da velocidade rotativa (MARTÍNEZ et al., 2005;
WOOD, 2011c), permanecendo assim praticamente constantes, procurouse
minimizar qualquer efeito semelhante ao aplicar o menor valor de λ ao menor
diâmetro, e viceversa. Caso contrário, aliás, o valor de rotação máxima seria alto do
ponto de vista operacional e físico – 630,628 rpm.
Além do CP, é otimizado também a massa da turbina com intuito de não só
diminuir possíveis custos de produção futuros, como também em oferecer menor
resistência ao movimento ou seja, a inércia da turbina. Este último, inclusive, foi
investigado por Amarante Mesquita et al., (2014) e cuja proposta é apresentada na
Equação (88) e utilizada neste trabalho. Já o cubo, por sua vez, tem seu diâmetro
determinado pela Equação (1). Devese salientar que a massa deste último e qualquer
densidade material foram desconsideradas, uma vez que tal abordagem foge do
escopo do atual trabalho.
𝐽 = 𝐵 ∑𝑚𝑖𝑟𝑖2
𝑁
𝑖=1
(88)
83
Em que i representa a seção transversal em questão, e m a massa entre a
seção transversal presente e anterior. Essa pode ser descrita segundo Equação (89),
desenvolvida por Mark Drela et al., (2006).
𝑚 = 𝜌𝑚𝐾𝐴𝑐𝑡𝑑𝑟 (89)
Onde ρm a densidade do material (considerada unitária), KA é o coeficiente de
área, normalmente igual a aproximadamente 0,6 (MARK DRELA et al., 2006), e t a
máxima espessura do perfil hidrodinâmico, que no caso do NACA 4412 é igual a 0,12
(ou 12%).
A multiobjetividade da função é feita segundo aplicado por Wood, (2011d), em
que se pondera cada uma das variáveis otimizadas – essas, importante salientar,
devem ser normalizadas para que sejam equiparadas. Além disso, devido a infinidade
de combinações, a potência nominal a ser gerada também pode assumir inúmeros
valores e assim, atender ou não a demanda de 252,78 W. Contra isso, um coeficiente
em forma de função parabólica minimiza a função objetivo quando a potência
produzida é a de projeto. Assim sendo, a função objetivo a ser minimizada é
apresentada segundo Equações (90) e (91).
𝑂 = 𝛾 ((1 − )(𝐶𝑃𝑟𝑒𝑓
𝐶𝑃) + (
𝐽
𝐽𝑟𝑒𝑓)) (90)
𝛾 = 9,999(𝑃
𝑃𝑟𝑒𝑓)
2
− 19,998(𝑃
𝑃𝑟𝑒𝑓) + 10 (91)
Onde o subscrito ref se refere aos valores de referência, γ o coeficiente
parabólico da demanda de potência e ε o coeficiente de peso. Inicialmente, o
coeficiente de peso é igual a 0,5, o valor de referência de CP é o de Betz (0,5973),
enquanto que a referência para a inércia de 9,438E5. Este fora determinado tomando
o valor encontrado por Vaz et al., (2018) e dividindo pela densidade do material de
fabricação da pá (SINGH, 2014). A tolerância para o critério de convergência passa a
ser igual a 1E6.
84
A otimização é realizada de forma híbrida e da seguinte maneira: iniciase
com algum dos modelos ABC, FPA ou PSO por 300 iterações; posteriormente, com o
melhor resultado obtido, aplicase o modelo SA em 100 iterações para evitar o mínimo
local e caminhar para o mínimo global (DEKKERS; AARTS, 1991). Todos iniciam
como uma população de 100 possíveis soluções e os modelos ABC e FPA, pelas suas
naturezas, são comumente simulados inúmeras vezes para vasculhar o maior número
de populações possíveis; neste caso, este número de repetições é 20. Essas e outras
constantes de cada um dos modelos podem ser vistas na Tabela 7.
Tabela 7 Constantes dos modelos de otimização.
Constante Valor
População inicial 100
ABC
Máximas tentativas 100
FPA
Probabilidade de troca 0,8
PSO
Termo de inércia 1
Constante de desaceleração 0,99
Constante cognitiva 2
Constante social 2
SA
Probabilidade de aceite inicial 0,8
Probabilidade de aceite final 0,001
Ciclos por esfriamento 50 Fonte: Próprio autor.
Assim, nas Figura 32, Figura 33 e Figura 34 são representadas o
desenvolvimento da solução ótima – no caso, o valor da função objetivo para cada
um dos modelos híbridos. A Tabela 8, por fim, apresenta os resultados obtidos.
85
Figura 32. Resultado da otimização do modelo híbrido ABCSA.
Fonte: Próprio autor.
Figura 33. Resultado da otimização do modelo híbrido FPASA.
Fonte: Próprio autor.
86
Figura 34. Resultado da otimização do modelo híbrido PSOSA.
Fonte: Próprio autor.
Tabela 8 Resultado dos designs ótimos para turbina.
Resultado ABCSA FPASA PSOSA
Valor objetivo (x104) 8,30785 8,30410 8,30408
Coeficiente de potência 0,475509 0,475503 0,475503
Inércia total (x105 kg*m2) 4,175588 4,176472 4,176527
Diâmetro do rotor (m) 0, 550777 0, 550836 0,550841
Diâmetro do cubo (m) 0,13770 0,13772 0,13772
Número de pás 3 3 3
Rotação (rpm) 370,5412 370,5412 370,5412
Λ 6,4846 6,4853 6,4854
Potência (W) 252,722 252,773 252,776
Tempo de simulação (s) 539,80 577,73 2220,14 Fonte: Próprio autor.
Tanto as figuras como a tabela acima apresentam resultados expressivos no
que diz respeito ao intuito do trabalho. Por meio das Figura 32, Figura 33 e Figura 34
podese notar que por mais que o método SA faça com que o mínimo local seja
evitado, a estimativa inicial é crucial para um bom resultado. Isso porque
aparentemente existem inúmero mínimos locais como resultado e que dificultam
87
alcançar o global. No entanto, na diferença do valor objetivo final entre os métodos:
tomando como referência o apresentado por PSOSA , o método ABCSA é defasado
em aproximadamente 0,04540% enquanto que o FPASA é 0,00024%. De um ponto
de vista de aplicação, esses números representam uma disparidade irrelevante.
Algo comum entre todos os designs é a velocidade rotativa e o número de
pás, cujo número é o limite máximo e o mínimo de domínio, respectivamente; ambas
constatações tendem a ser uma orientação sempre que possível. Devese destacar,
que por esse motivo a razão entre as velocidades tangencial e de corrente livre é
naturalmente alterada, e junto disso, o ônus do decréscimo do coeficiente de potência.
Tal comportamento é condizente com a teoria e o comportamento apresentado pela
Figura 10.
No que diz respeito à qualidade dos métodos híbridos, há dois quesitos a
serem julgados: tempo de simulação e qualidade do resultado. E ambos aparentam
ser inversamente proporcionais. Enquanto o método PSOSA teve duração quase
cinco vezes maior que os demais, o valor objetivo foi o menor atingido. Isso se dá
muito provavelmente pela natureza de seu método e que neste caso, devido à grande
quantidade de combinações, passa a ser necessário – em vasculhar boa parte do
domínio.
Por outro lado, por mais que os designs de ABCSA e FPASA não tenham
alcançado o suposto mínimo global atingido pelo método PSOSA em uma ordem de
grandeza irrisória, ambos têm um custo computacional muito menor. Já o primeiro
ainda apresenta resultado interessante: o valor objetivo muito próximo do segundo,
porém traduzido em uma turbina ligeiramente mais eficiente e com menor massa e
diâmetro, mas que por outro lado, não atende completamente a demanda de potência
– mas novamente, em uma escala milesimal. Isso pode ter origem em ambas suas
formulações, cujos métodos de busca aleatória podem necessitar de uma maior
quantidade de tentativas; talvez uma alteração no número de iterações, no número de
simulações ou na população inicial possa beneficiálos, mas tal análise foge do
escopo do trabalho.
Dessa forma, sob as condições ambiente impostas, um projeto que equilibre
a potência e a massa através de um custo computacional mínimo tem como melhor
design o resultado de ABCSA. Merece destaque o coeficiente de potência:
aproximadamente 17,92% menor que o ideal de Betz. As Figura 35, Figura 36 e Figura
88
37 a seguir trazem as distribuições no comprimento da pá de corda, solidez e torção
de tal design.
Figura 35. Distribuição da corda ao longo do comprimento da pá.
Fonte: Próprio autor.
89
Figura 36. Distribuição da torção ao longo do comprimento da pá.
Fonte: Próprio autor.
Figura 37. Distribuição da solidez ao longo do comprimento da pá.
Fonte: Próprio autor.
90
Se no lugar dos coeficientes de sustentação e arrasto obtidos pelo XFOIL
(DRELA, 2001), fossem usados resultados do processo de validação do perfil isolado
da Seção 4.1.2, as distribuições de corda e torção se apresentam segundo as Figura
38 e Figura 39.
Figura 38. Comparação da distribuição da corda ao longo do comprimento da pá entre os
coeficientes obtidos pelo XFOIL (DRELA, 2001)e os validados na Seção 4.1.2.
Fonte: Próprio autor.
91
Figura 39. Comparação da distribuição da torção ao longo do comprimento da pá entre os coeficientes obtidos pelo XFOIL (DRELA, 2001) e os validados na Seção 4.1.2.
Fonte: Próprio autor.
Naturalmente, os coeficientes têm influência direta nos esforços ao longo da
pá, o quê segundo a Figura 38 é traduzido em uma corda maior em cada seção
transversal. Consequentemente, o valor de inércia passa a ser quase 15,30% maior
com relação ao design ótimo segundo os dados do software XFOIL (DRELA, 2001).
Além disso, o coeficiente de potência também sofre decréscimo: mais de 5,30%,
passando a 0,4505. Isso demonstra, de forma análoga, a importância que se têm a
escolha do perfil hidrodinâmico – que acarreta intrinsicamente nos coeficientes de
arrasto, sustentação e de pressão mínimo em projetos de turbinas. A Tabela 9 traz
os resultados do processo de otimização segundo os coeficientes em questão.
92
Tabela 9 Resultado da geometria otimizada para turbina, através do modelo ABCSA segundo os coeficientes validados na Seção 4.1.2.
Resultado XFOIL DFC Diferença (%)
Coeficiente de potência 0,475509 0,448012 5,78
Inércia total (x105 kg*m2) 4,175588 4,692843 12,39
Diâmetro do rotor (m) 0, 550777 0,567617 3,06
Diâmetro do cubo (m) 0,13770 0,14191 3,06
Número de pás 3 3 0
Rotação (rpm) 370,5412 370,5412 0
Λ 6,4846 6,6829 3,06
Potência (W) 252,722 252,891 0,07 Fonte: Próprio autor.
Pela qualidade do design otimizado de forma singular à demanda de potência
especificada, desperta uma possibilidade de algum layout para geração de uma maior
quantidade de energia elétrica. Tal hipótese pode ser analisada levando em conta a
quantidade de pessoas que a turbina atenderia. Tanto o design quanto a demanda
estão presentes na Tabela 10, seguindo mesmo roteiro anterior e apenas sob a
delimitação de diâmetro: caso o valor máximo ou mínimo exceda a profundidade com
decréscimo de 10 ou 20 cm, esse critério passa a ser o limite.
Tabela 10 Resultado da geometria otimizada para turbina, através do modelo ABCSA
segundo diferentes números de pessoas como demanda.
Demanda de potência (W)
1 pessoa 2 pessoas 4 pessoas 252,78 505,56 1011,11
Resultados
Coeficiente de potência 0,475509 0,479950 0,480096
Inércia total (x105 kg*m2) 4,175588 10,931000 18,020000
Diâmetro do rotor (m) 0, 550777 0,775759 0,900000
Diâmetro do cubo (m) 0,13770 0,19395 0,22501
Número de pás 3 3 3
Rotação (rpm) 370,5412 349,6910 349,6910
Λ 6,4846 8,6195 10,0000
Potência (W) 252,77 506,03 681,310 Fonte: Próprio autor.
Interessante notar como mais uma vez os resultados evidenciam a
sensibilidade das condições de operação. É observado como as condições externas
93
permitem atender perfeitamente 2 pessoas, e que isso se dá não só pelo aumento do
diâmetro, algo natural, como também pela variação da rotação. Isso não é visto para
uma demanda de 1011,11 W, muito provavelmente pelo fato de tal potência ser
impossível de ser alcançada segundo condições ambientais e de projeto impostas.
4.2.1.2 Simulação DFC da turbina hidrocinética
A partir das distribuições de corda e de torção do layout otimizado (Figura 35
e Figura 36), uma simulação computacional é conduzida com o objetivo de se
comparar a possível potência gerada por tal turbina, com o valor esperado. Devese
lembrar que não há validação experimental das configurações da simulação, e que tal
procedimento pode ser realizado em trabalhos futuros. Assim, um modelo
tridimensional virtual é feito com auxílio do software Solidworks (DASSAULT
SYSTÈMES, ) e pode ser visto na Figura 40 abaixo.
Figura 40. Representação do rotor completo segundo design otimizado.
Fonte: Próprio autor.
Entretanto, uma simulação computacional do rotor inteiro como na Figura 40
pode tornála altamente custosa do ponto de vista numérico. Contra isso, a construção
da malha, realizada pelo software ICEM, parte da importação de apenas um terço do
rotor. Em outras palavras, é optado por se simular apenas uma das pás – alternativa
esta já realizada por outro autor (SILVA et al., 2017) , cuja operação da turbina inteira
94
é replicada pela ferramenta de periodicidade. No que diz respeito à reprodução do
movimento rotacional, a estratégia utilizada fora a chamada MRF (do inglês Moving
Reference Frame) por demandar uma menor quantidade de elementos na malha sem
prejuízo no resultado final desse tipo de simulação, como constatado por outros
autores em trabalhos semelhantes (DASKIRAN et al., 2017; SILVA et al., 2017;
HOGHOOGHI; DURALI; KASHEF, 2018; ABUAN; HOWELL, 2019). Uma
representação do domínio é apresentada na Figura 41 abaixo. Figura 41. Representação esquemática do domínio de simulação DFC da turbina hidrocinética.
Fonte: Próprio autor.
Foram desenvolvidas três malhas para o processo de independência de
malha – que será precedido pelo processo de independência de passo de tempo ,
seguindo o parâmetro de y+ igual a 3,5 validado na Seção 4.1.2 e próximo ao adotado
em outro trabalho (HOGHOOGHI; DURALI; KASHEF, 2018). O número de elementos
de cada uma delas é apresentado na Tabela 11. No que diz respeito às regiões citadas
na Figura 41, Inlet e Outlet correspondem à entrada e saída do fluido,
respectivamente; FLUID_ROTATING e FLUID_STATIC aos volumes rotacional e
estático (já que o processo de discretização do MRF se baseia na aceleração de
Coriolis de um corpo em movimento rotacional em relação a um referencial estático);
95
BLADE e HUB à pá e ao cubo do rotor; SYM_EXT à superfície externa de simetria do
domínio, e Periodic e Shadow os pares relativos à periodicidade para cada volume.
Tabela 11 Detalhes das malhas construídas para a simulação da turbina.
Denominação Número de elementos
Razão entre malhas (τ)
Malha D 2067506 1,48
Malha E 3065748 1,50
Malha F 4591924 Fonte: Próprio autor.
Para as simulações o software OpenFOAM (OPENFOAM, 2019) é utilizado,
e um resumo das condições iniciais para cada região do domínio pode ser visto na
Tabela 12. Naturalmente, os parâmetros que tangem às condições ambientais são
mantidos segundo descrito na Seção 4.2.1.1. O modelo de turbulência utilizado é o k
ω SST, validado na Seção 4.1.2. Já o volume em movimento (FLUID_ROTATING) é
configurado segundo a velocidade rotacional da Seção 4.2.1.1: 370,5412 rpm, ou
aproximadamente 38,8 rad/s.
Tabela 12 Condições iniciais para as simulações da turbina.
Região U p K omega nut
internalField (0,0 1,64788 0,0) 0,0 uniform 0,001 uniform 10,0 uniform 1e5
Inlet uniformFixedValue
constant (0,0 1,64788 0,0)
zeroGradient
turbulentIntensityKineticEnergyInlet
uniform 1,0 intensity 0,001
turbulentMixingLengthFrequencyInlet uniform 10,0
mixingLength 1,0
calculated uniform
0,0
Outlet zeroGradient fixedValue uniform 0,0
turbulentIntensityKineticEnergyInlet
uniform 1,0 intensity 0,001
turbulentMixingLengthFrequencyInlet uniform 10,0
mixingLength 1,0
calculated uniform
0,0
BLADE fixedValue
uniform (0,0 0,0 0,0) zeroGradient
kqRWallFunction uniform 1e10
omegaWallFunction value 10,0
nutUWallFunction uniform
0,0
HUB fixedValue
uniform (0,0 0,0 0,0) zeroGradient
kqRWallFunction uniform 1e10
omegaWallFunction value 10,0
nutUWallFunction uniform
0,0
96
SYM_EXT Symmetry symmetry Symmetry symmetry symmetry
Periodic/Shadow cyclicAMI cyclicAMI cyclicAMI cyclicAMI cyclicAMI Fonte: Próprio autor.
Onde U é a velocidade, p a pressão, k a energia cinética turbulenta, omega a
taxa de dissipação específica de turbulência e nut a viscosidade cinemática turbulenta.
Devido à complexidade do escoamento – e consequentemente da simulação
, houve dificuldade em configurar os parâmetros de resolução das equações e de
seus métodos iterativos. Por isso, uma alternativa recorrida foi a discretização do
tempo (ou seja, uma simulação transiente com solver PISO) em todo o procedimento
de alcance do resultado final. A medida referencial para esse tipo de simulação é o
tempo de residência, que é o tempo físico mínimo que uma partícula leva para
atravessar o domínio, e fora estimado em 3,5 s. O passo de tempo que melhor capta
as mudanças do escoamento e o desenvolvimento da operação da turbina é
determinado pela independência de passo de tempo. Os valores estabelecidos foram
selecionados de forma semelhante aos trabalhos de Wang; Yin; Yan (2019) e Abuan
e Howell (2019), e são: 1E4 (0,1°/s), 4E4 (0,4°/s) e 8E4 (0,9°/s). Dessa forma, as
simulações iniciais com intuito de se obter um campo para essa análise são
conduzidas com o passo de tempo intermediário.
Assim, iniciase com uma simulação durante um tempo de residência, com
resíduo igual a 1E5, equações resolvidas em 1ª ordem (ou upwind) e fatores de sub
relaxação: 0,1 para pressão e 0,5 para velocidade, energia cinética turbulenta e taxa
de dissipação específica de turbulência. Deste campo, alterase a ordem de solução
para 2ª ordem (ou linearUpwind) para um segundo tempo de residência de simulação.
Finalmente, novamente deste campo, três novas simulações por um novo tempo de
residência cada são realizadas para se obter os resultados médios de cada um dos
passos de tempo em análise. O parâmetro a ser investigado é o torque da pá
(exportado das simulações), e seus resultados médios, bem como o desenvolvimento
desta variável ao longo do tempo para cada passo de tempo podem ser vistos na
Tabela 13 e Figura 42. O torque médio ideal e esperado, usado como referência para
cálculo da diferença relativa percentual da Tabela 13, pode ser deduzido pela
Equação (92) abaixo.
𝑃 = 𝑄 ∗ ω (92)
97
Tabela 13 Resultado das simulações de DFC para independência de passo de tempo.
Passo de tempo (s) Torque médio (N*m) Diferença (%)
Ideal 6,51
1E4 5,40 17,13
4E4 8,27 27,03
8E4 7,53 15,63 Fonte: Próprio autor.
Figura 42. Comportamento temporal do torque para cada um dos passos de tempo.
Fonte: Próprio autor.
A princípio, a Tabela 13 aponta para uma inconsistência dos resultados, já
que a evolução destes não segue o refinamento do passo de tempo. Isso porque o
torque médio apresentado pelo passo de tempo intermediário fora o mais distante do
ideal, e não o mais grosseiro. Entretanto, percebese pela Figura 42 que por mais que
dois tempos de residência já tivessem sido simulados, o tempo físico simulado com
as configurações mais precisas pode ter sido insuficiente para a estabilização de um
padrão para a variação temporal do torque. Neste sentido, simulouse 2 tempos de
residência adicionais com o passo de tempo de 8E4 (0,9°/s), e pode ser visto na
Figura 43.
98
Figura 43. Comportamento temporal do torque para um passo de tempo de 8E4 (0,9°/s) ao longo de três tempos de residência.
Fonte: Próprio autor.
Mesmo com quase 10,5 s de tempo físico simulado, é possível inferir que o
torque não alcançou um regime estatisticamente estabelecido. Por outro lado, vale
ressaltar a melhora apresentada no seu valor médio: de 7,53 N*m em apenas 1 tempo
de residência, para 7,20 N*m considerando todos os 3 tempos de residência, uma
melhora relativa percentual de 5%. Isso mostra que uma simulação de operação mais
longa pode não só alcançar um padrão de comportamento, como também uma média
mais próxima do ideal. Entretanto, através dos recursos disponíveis na Universidade
tal análise é inviável. Toda a investigação contou com um cluster de 36 núcleos que
permitia simular em paralelo, mas que devido à complexidade do problema ainda
apresentava um custo computacional alto: o tempo gasto em média, por tempo de
residência e para cada passo de tempo, do mais grosseiro ao mais refinado, foram 2
semanas, 4 semanas, e 8 semanas, respectivamente. Consequentemente, a análise
de independência de malha também fora prejudicada e não fora possível prosseguir
com ela.
Ainda assim, é possível analisar os efeitos do escoamento enquanto
resultados da simulação DFC. O pósprocessamento é feito através do software VisIt
(ADVANCED SIMULATION AND COMPUTING INITIATIVE (ASCI), 2002), e o
resultado final para o passo de tempo mais refinado é tomado como base. As Figura
44 abaixo trazem, na seção da transversal ao longo da altura da turbina, os campos
de velocidade por componente: vertical, horizontal e axial, respectivamente.
99
Figura 44. Campos de velocidade instantânea por componente (da esquerda para direita): vertical, horizontal e axial.
Fonte: Próprio autor.
Uma vez que a velocidade rotacional é no sentido horário, percebese o efeito
que a pá da turbina causa enquanto obstáculo. Ao mesmo tempo em que a água é
acelerada horizontalmente, ela descende verticalmente e toma o sentido axial
contrário. Neste sentido, a Figura 45 abaixo apresenta as linhas de corrente
impactadas pela operação da turbina, a complexidade e imprevisibilidade do
escoamento ficam evidentes. Isso é um indício claro tanto da dificuldade em prever
seu funcionamento, como também de se obter um resultado de simulação DFC
estatisticamente estabelecido.
100
Figura 45. Linhas de corrente decorrentes da operação da turbina, coloridas pela magnitude da velocidade média.
Fonte: Próprio autor.
A consequência direta desse funcionamento é a formação de esteiras
helicoidais, representado na Figura 46 através das isosuperfíces de velocidade
instantânea.
101
Figura 46. Visualização da esteira helicoidal através da isosuperfíce de velocidade instantânea ao longo domínio.
Fonte: Próprio autor.
Através de uma análise minuciosa, a Figura 47 apresenta uma aproximação
à pá, até que seja possível observar o efeito singular da ponta de pá. E isso fica mais
evidente se tomarmos as linhas de corrente que a contornam, representado na Figura
48, e colorida pela velocidade média. É notável o quanto, à medida que se aproxima
da ponta da pá, a velocidade da água aumenta drasticamente. Não à toa, esse
impacto tem sido um dos principais alvos de desenvolvimento de teoria como correção
ao desempenho da turbina.
102
Figura 47. Aproximação à ponta da pá, e visualização de seu efeito pela isosuperfíce de velocidade instantânea.
Fonte: Próprio autor.
Figura 48. Linhas de corrente contornando a pá, coloridas pela magnitude da velocidade
média.
Fonte: Próprio autor.
Se for tomado um ponto ao longo da altura da pá, a visualização das
componentes da velocidade ressalta a importância de se projetar adequadamente o
103
perfil de pá. Através da Figura 49, é visível o ponto de estagnação do fluido, que
quando bem estimado gera a otimização das forças de arrasto e sustentação.
Figura 49. Campos de velocidade instantânea por componente em um ponto na altura da pá
(da esquerda para direita): vertical, horizontal e axial.
Fonte: Próprio autor.
Em contraste, por se tratar de um modelo de turbina que não conta com a
energia do escoamento em forma de pressão, e sim exclusivamente com a sua
energia cinética, a diferença de pressão ao longo do domínio é ínfima. Mas é
interessante investigar este campo no mesmo ponto da Figura 49, e pode ser visto na
Figura 50. Notase que o intuito de se aplicar um perfil aero/hidrodinâmico é atingido,
já que a diferença líquida de pressão da sua porção inferior e superior é positiva, e
atua como força motriz do movimento rotacional.
Figura 50. Campo de pressão instantânea em um ponto na altura da pá.
Fonte: Próprio autor.
104
4.2.2 SOLIDEZ
4.2.2.1 Simulações DFC dos perfis em série
O ponto em comum entre os trabalhos citados na Seção 2.2.1.9 é que o
parâmetro de solidez tem relevância sob os coeficientes de escoamento. Isso se dá
por conta da interação entre o fluido e as pás algo que, no que tange aos efeitos nos
parâmetros, se despreza nas teorias de previsão do desempenho das turbinas. É
neste âmbito que se busca uma possível correlação entre a solidez e os coeficientes
de arrasto, sustentação e pressão mínima, com intuito de não só mensurar o quanto
efetivamente é gerado em potência elétrica, como também dimensionar
adequadamente os esforços em torno das seções transversais.
Assim, são realizadas as seguintes simulações numéricas bidimensionais por
meio de DFC para que os melhores parâmetros de projeto sejam determinados:
• 10 valores de solidez: entre 0,1 e 1,0, com intervalo de 0,1;
• 6 ângulos de ataque, para cada solidez: 0°, 1,25°, 2,5°, 3,75°, 5° e
10°.
As malhas são construídas utilizando o software ICEM e segundo esquema
da Figura 51, em que se varia a distância entre as porções do perfil e o valor de corda
se mantém igual a 1,0. Buscase uma densidade regular de número de nós na altura
do domínio conforme se varia a distância, com intuito de evitar qualquer discrepância
numérica entre as malhas e a manutenção do valor y+ de 3,5 validado na Seção 4.1.2.
O número de elementos para cada uma das malhas, bem como a qualidade mínima
avaliada segundo critério de Determinante 3x3x3, podem ser vistos na Tabela 14.
105
Figura 51. Representação esquemática da construção do domínio numérico para simulação da solidez.
Fonte: Próprio autor.
Tabela 14 Quantidade de elementos para cada uma das malhas de solidez.
Distância (m) Solidez Número de elementos Qualidade minima
10,00 0,1 1118448 0,8990
5,00 0,2 678318 0,9250
3,33 0,3 540238 0,9340
2,50 0,4 514348 0,9130
2,00 0,5 479828 0,9470
1,67 0,6 445308 0,9460
1,43 0,7 419418 0,9520
1,25 0,8 419418 0,9530
1,11 0,9 406473 0,9350
1,00 1,0 406473 0,9240 Fonte: Próprio autor.
No que diz respeito às regiões citadas na Figura 51, Inlet e Outlet
correspondem à entrada e saída do fluido, respectivamente; NACA_TOP e
NACA_BOTTOM as metades superior e inferior do perfil NACA 4412 e Periodic e
Shadow os pares relativos à periodicidade.
Para as simulações o software OpenFOAM (OPENFOAM, 2019) é utilizado.
As condições iniciais do problema são determinadas ao passo que mais se aproximem
da operação apresentada no procedimento experimental e na validação do perfil
isolado. Dessa forma, um resumo delas pode ser visto na Tabela 15.
106
Tabela 15 Condições iniciais para as simulações de solidez. Região U p k omega nut
internalField (0,0 0,0 0,0) 0,0 uniform 0,001 uniform 10,0 uniform 1e5
Inlet uniformFixedValue constant (UxUy 0,0)
zeroGradient
turbulentIntensityKineticEnergyInlet
uniform 1,0 intensity 0,001
turbulentMixingLengthFrequencyInlet uniform 10,0
mixingLength 1,0
calculated uniform
0,0
Outlet zeroGradient fixedValue uniform 0,0
turbulentIntensityKineticEnergyInlet
uniform 1,0 intensity 0,001
turbulentMixingLengthFrequencyInlet uniform 10,0
mixingLength 1,0
calculated uniform
0,0
NACA_TOP/ NACA_BOTTOM
fixedValue uniform (0,0 0,0 0,0)
zeroGradient kqRWallFunction
uniform 1e10 omegaWallFunction
value 10,0
nutUWallFunction uniform
0,0
Periodic/Shadow cyclic cyclic cyclic cyclic cyclic Fonte: Próprio autor.
Onde U é a velocidade, p a pressão, k a energia cinética turbulenta, omega a
taxa de dissipação específica de turbulência e nut a viscosidade cinemática. A
mudança do ângulo de ataque se dá pelas componentes de velocidade na entrada,
segundo Equação (86).
O procedimento contou com uma simulação inicial para 0°, regime
permanente (solver SIMPLE), resíduo igual a 1E6 e equações resolvidas em 1ª ordem
(ou upwind) para obtenção de um campo inicial; deste campo, alterase para solução
de 2ª ordem (ou linearUpwind) em regime transiente (solver PISO), para resíduo igual
a 1E5. O passo de tempo varia de acordo com a complexidade do escoamento, entre
1E3 e 5E5; os fatores de subrelaxação para o regime transiente são: 0,1 para
pressão e 0,5 para velocidade, energia cinética turbulenta e taxa de dissipação
específica de turbulência. Exclusivamente para 0°, é simulado apenas um tempo de
residência – 10 segundos físicos e posteriormente outro para que se obtenha os
resultados médios; para os demais, como é tomado a simulação de ângulo de ataque
anterior como campo inicial, são simulados dois tempos de residência – e da mesma
forma novamente outro para se obter os resultados médios.
São exportados os coeficientes de arrasto e sustentação tanto para a porção
inferior como superior, em que o valor total é a soma dos respectivos das duas
metades. O vetor direção conforme variação do ângulo de ataque é determinado
segundo Eq. (87). Já o coeficiente de pressão mínimo é obtido em pósprocessamento
107
pelo software VisIt (ADVANCED SIMULATION AND COMPUTING INITIATIVE
(ASCI), 2002) por meio da Equação (93); como nenhuma pressão fora imposta, a
pressão de referência é a padrão do software OpenFOAM (OPENFOAM, 2019): igual
a 0.
𝐶𝑝𝑚𝑖𝑛 =𝑝 − 𝑝01
2𝜌𝑈0
2 (93)
Onde p é a pressão local do campo. Os resultados dos coeficientes de arrasto
e sustentação são apresentados nas Figura 52, Figura 53 e Figura 54, e dos
coeficientes de pressão mínimo na Tabela 16. Devese ressaltar que o valor de solidez
igual 0 provém dos resultados de validação do perfil isolado na Seção 4.1.2 para os
ângulos de 0°, 5° e 10°, enquanto que para os demais foise obtido por meio de
interpolação linear. Já as Figuras 55, Figuras 56 e Figuras 57 apresentam,
respectivamente, os contornos de velocidade média, pressão média e energia cinética
turbulenta, de algumas das simulações de solidez realizadas.
Figura 52. Variação do coeficiente de arrasto segundo ângulo de ataque para os valores de
solidez entre 0 e 1.
Fonte: Próprio autor.
108
Figura 53. Variação do coeficiente de sustentação segundo ângulo de ataque para os valores de solidez entre 0 e 1.
Fonte: Próprio autor.
Figura 54. Variação da razão entre o coeficiente de sustentação e arrasto segundo ângulo de ataque para os valores de solidez entre 0 e 1.
Fonte: Próprio autor.
109
Tabela 16 Resultado dos coeficientes de pressão mínimo segundo ângulo ataque para os valores de solidez entre 0 e 1.
Solidez 0° 1,25° 2,50° 3,75° 5° 10°
0,0 0,8048 0,9921 1,1794 1,3667 1,5540 4,5100
0,1 0,7703 0,8985 1,0267 1,1548 1,2830 3,4430
0,2 0,7438 0,8128 0,8823 0,9732 1,1310 2,6030
0,3 0,7260 0,7858 0,8455 0,9130 1,0110 1,9470
0,4 0,8420 0,7668 0,8189 0,8715 0,9393 1,5430
0,5 0,8530 0,7534 0,7993 0,8449 0,9002 1,3780
0,6 0,9163 0,7499 0,7951 0,8401 0,8849 1,3150
0,7 0,9759 0,7802 0,7794 0,8215 0,8693 1,2030
0,8 0,9695 0,7917 0,7779 0,8166 0,8619 1,1950
0,9 1,1160 0,9735 0,8247 0,8581 0,8730 1,0500
1,0 1,1460 1,0080 0,8496 0,8604 0,8717 1,0230 Fonte: Próprio autor.
110
Figuras 55. Campos de velocidade média para o ângulo de ataque igual à 0°, seguindo a solidez de 0,1 (a) até 1,0 (j).
(a) (b) (c)
(d) (e) (f)
(g) (h) (i)
(j)
Fonte: Próprio autor.
111
Figuras 56. Campos de pressão média para o ângulo de ataque igual à 5°, seguindo a solidez de 0,1 (a) até 1,0 (j).
(a) (b) (c)
(d) (e) (f)
(g) (h) (i)
(j)
Fonte: Próprio autor.
112
Figuras 57. Campos de energia cinética turbulenta para o ângulo de ataque igual à 10°, seguindo a solidez de 0,1 (a) até 1,0 (j).
(a) (b) (c)
(d) (e) (f)
(g) (h) (i)
(j)
Fonte: Próprio autor.
113
É evidente pelas Figura 52, Figura 53 e Figura 54 que os coeficientes de
escoamento são diretamente dependentes dos valores de solidez: em suma,
enquanto os números da sustentação decrescem, os de arrasto aumentam. A origem
disso consiste na proximidade entre os perfis contrair o escoamento entre eles, o que
beneficia o gradiente de pressão favorável e acelera o fluido, como pode ser visto na
Figuras 55. Assim, o escoamento da camada limite apresenta uma maior energia
cinética em resistência ao gradiente de pressão adverso, que em condições de
escoamento livre, pode ocasionar o descolamento dessa mesma camada. Tal
tendência pode ser observada conforme o valor de solidez diminui, em que a zona de
recirculação da borda de fuga devido ao gradiente de pressão adverso do
aero/hidrofólio se expande e a energia cinética turbulenta é intensificada (Figuras 57).
Uma consequência direta de tal efeito é o atraso do estol, por exemplo, sendo que já
fora observado em outros trabalhos (CEBRIÁN; ORTEGACASANOVA;
FERNANDEZFERIA, 2013; LABIGALINI et al., 2019). E finalmente, há também o
ônus do decaimento das pressões resultantes (entre porção superior e inferior do
perfil) envolvidas no aero/hidrofólio, visto nas Figuras 56, e que naturalmente prejudica
os coeficientes aero/hidrodinâmicos.
4.2.2.2 Algoritmo
O balanço dos efeitos citados se reflete em uma diminuição da razão entre a
sustentação e arrasto, com tendência inversamente proporcional ao aumento da
solidez. Isso fica evidente na Figura 58, em que são plotadas a variação de razão
ótima entre sustentação e arrasto e seu ângulo de ataque segundo os valores de
solidez.
114
Figura 58. Variação da razão ótima entre os coeficientes de sustentação e arrasto e seus ângulos de ataque para os valores de solidez entre 0 e 1.
Fonte: Próprio autor.
A partir de tais ângulos de ataque, ótimo para cada valor de solidez, é possível
interpolar funções que traduzam a variação dos coeficientes de arrasto, sustentação
e pressão mínima segundo a solidez. Isso é feito através da ferramenta interp1d,
presente na biblioteca scipy da linguagem de programação Python. Ambas as funções
são incluídas no processo iterativo de obtenção da corda e torção, ao passo que para
cada valor de solidez, a próxima iteração tem os coeficientes citados atualizados.
Para se ter dimensão da diferença que a consideração dos coeficientes em
estudo como função da solidez pode provocar, um perfil de pá é obtido com os
mesmos parâmetros de design ótimo encontrados na Seção 4.2.1.1. A comparação
da distribuição da corda e da torção podem ser vistos nas Figura 59 e Figura 60,
respectivamente.
115
Figura 59. Comparação da distribuição de corda com e sem correção de solidez
Fonte: Próprio autor.
Figura 60. Comparação da distribuição de torção com e sem correção de solidez.
Fonte: Próprio autor.
Percebese que, no que diz respeito à estrutura, a mudança principal ocorre
no começo da pá, onde em razão de uma maior solidez os coeficientes de escoamento
116
são mais afetados; isso leva a uma diminuição dos esforços normais (Equação (13)),
acarretando a um aumento da corda (segundo Equação (59)) e inércia, está para
quase 30% em comparação ao design ótimo obtido por meio dos coeficientes do
software XFOIL (DRELA, 2001); este número cai para 6% se for tomado os
coeficientes validados na Seção 4.1.2. Além disso, os parâmetros de desempenho
também são prejudicados da mesma forma: há uma queda de quase 13% na potência
gerada e de 8,15% no coeficiente de potência, respectivamente segundo uso dos
coeficientes. Em termos práticos, isso representa um déficit de 23,66 kWh em energia
elétrica gerada ao final do mês. Contra isso, um novo processo de otimização deve
ser feito com o propósito de se evitar tal diferença.
A otimização é feita novamente segundo o método validado de acordo com a
Tabela 8, e o processo iterativo é representado na Figura 61 e uma comparação entre
o resultado obtido na Seção 4.2.1.1 sem correção de solidez – com referência aos
valores encontrados pelo XFOIL e o presente corrigido pode ser visto na Tabela 17.
Figura 61. Resultado da otimização do modelo híbrido ABCSA corrigido segundo resultados
de solidez.
Fonte: Próprio autor.
117
Tabela 17 Comparação dos resultados de design ótimos para turbina com e sem correção de solidez.
Resultado XFOIL DFC DFC (corrigido) Diferença (%)
Coeficiente de potência 0,475509 0,448012 0,416217 12,47
Inércia total (x105 kg*m2) 4,175588 4,692843 5,848779 40,07
Diâmetro do rotor (m) 0, 550777 0,567617 0,588675 6,88
Diâmetro do cubo (m) 0,13770 0,14191 0,14717 6,88
Número de pás 3 3 3 0,00
Velocidaderotativa (rpm) 370,5412 370,5412 370,5412 0,00
Λ 6,4846 6,6829 6,9308 6,88
Potência (W) 252,722 252,891 252,699 0,02
Tempo de simulação (s) 539,80 496,83 692,98 122,72 Fonte: Próprio autor.
Notase pelos resultados da Tabela 17 que a diferença nos parâmetros de
design ocasionada pela correção dos coeficientes em função da solidez não é irrisória.
Longe disso, uma vez que é possível inferir, que a queda de quase 13% do coeficiente
de potência corrigido seja compensada pelo aumento dos parâmetros geométricos
com a finalidade de manter a potência nominal: praticamente 7% do diâmetro, e
surpreendentes quase 40% da inércia.
Tal resultado singular se mostra relevante mesmo com os baixos valores de
solidez que a turbina apresenta (máximo de 0,25 Figura 37), e como pôde ser visto
na Figura 54, essa faixa caminha para números parecidos aos do perfil isolado. Em
contrapartida, o tempo de simulação teve aumento expressivo, o que do ponto de vista
computacional poderia pôr em xeque correções nesse sentido; mas como tal
aplicação se concentra na fase de projeto, o tempo gasto acaba sendo mínimo quando
comparado a outros processos. As Figura 62, Figura 63 e Figura 64 a seguir trazem
as distribuições no comprimento da pá de corda, torção e solidez para um design ótimo
corrigido pelos efeitos da solidez.
118
Figura 62. Distribuição da corda ao longo do comprimento da pá, corrigido pelos efeitos da solidez.
Fonte: Próprio autor.
Figura 63. Distribuição da torção ao longo do comprimento da pá, corrigido pelos efeitos da solidez.
Fonte: Próprio autor.
119
Figura 64. Distribuição da solidez ao longo do comprimento da pá, corrigido pelos seus efeitos
no escoamento.
Fonte: Próprio autor.
Vale ressaltar que tal desfecho complementa as conclusões dos autores
acerca dos trabalhos de Islam; Islam (1994) e Singh (2014), onde a influência da
solidez se mostrou proeminente em turbinas com alto valor de solidez e baixo valor
de razão entre as velocidades de corrente livre e de ponta de pá, evidenciando que
mesmo para valores moderados de solidez as alterações geométricas necessárias
para a manutenção da potência de projeto são relativamente grandes. Ressaltase
ainda que esse efeito não deve ser desprezado em projetos de geração de energia
localizada, já que em tais condições o acesso à recursos ambientais e financeiros é
na maioria das vezes restrito, e qualquer precisão deve ser explorada.
120
4.2.3 RESSALTO
4.2.3.1 Simulações DFC dos perfis com ressalto
Protuberâncias que têm como objetivo estimular a transição à turbulência na
camada limite, denominadas neste trabalho como ressaltos no ataque do perfil, são
artifícios comuns e já estudados por outros autores (DIETZ; MAI; GEISSLER, 2008;
PECHLIVANOGLOU, 2013), e têm como principal objetivo o adiamento –
consequentemente novo pico do valor de sustentação – do efeito de estol.
Assim, o estudo é previamente analisado para o perfil isolado, e
posteriormente com eles em cascata. Um ressalto de formato triangular equilátero é
adicionado, com 1 mm de altura e à uma distância de 5% da corda com relação à
borda de ataque. Os elementos de malha inseridos serão apenas para contornar o
ressalto, ou seja, não serão acrescentados elementos na direção ortogonal ao perfil.
Partiuse então dos resultados de independência de malha da Seção 4.1.2, e
utilizouse o mesmo procedimento descrito, as mesmas condições iniciais presentes
na Tabela 4 e a mesma metodologia de cálculo dos coeficientes. A comparação dos
resultados com e sem o ressalto no perfil isolado podem ser visualizados nas Figura
65 e Figura 66 abaixo.
121
Figura 65. Comparação do coeficiente de arrasto do perfil com e sem ressalto.
Fonte: Próprio autor.
Figura 66. Comparação do coeficiente de sustentação do perfil com e sem ressalto.
Fonte: Próprio autor.
122
Independente do regime do escoamento (plenamente estabelecido ou estol),
é claro pelas Figura 65 e Figura 66 a piora que o ressalto acarreta nos coeficientes
fluidodinâmicos. No que diz respeito ao regime plenamente estabelecido, destacase
a mudança no comportamento dos valores. Enquanto o coeficiente de arrasto, que
antes apresentava um padrão explícito de crescimento, passa a ser imprevisível; por
outro lado, o coeficiente de sustentação deixa de uma tendência linear, e passa a uma
nãolinearidade. Enquanto o coeficiente de arrasto teve o valor praticamente
inalterado com relação ao ângulo de ataque de 5º, destoando de qualquer padrão que
se formava, o coeficiente de sustentação tem a maior disparidade deste regime com
relação ao perfil sem ressalto, e antecessor de um inesperado e desproporcional pico
no ângulo de ataque de 12º.
Essa discrepância de valores do ângulo de ataque de 10º, que chega a
21,37% para o coeficiente de sustentação, é justificado pelos resultados das Figuras
67 e Figuras 68. Notase pelas primeiras que tanto antes, como logo depois do
ressalto há um aumento considerável da energia cinética turbulenta com relação ao
resultado da simulação sem o ressalto, que pode ser compreendido como uma
dissipação de energia do fluido. E tal efeito se prolonga até a cauda do perfil, em que
há uma diminuição clara da esteira do escoamento. Naturalmente, o reflexo disso é
visualizado nos campos de velocidade média nas Figuras 68: o escoamento se
desenvolve à uma velocidade menor, o descolamento é antecipado, e a região de
pressão adversa ampliada. Consequentemente, a pressão resultante de sustentação
é diretamente prejudicada.
123
Figuras 67. Campos de energia cinética turbulenta para o ângulo de ataque igual à 10°(a) sem e (b) com o ressalto.
(a)
(b)
Fonte: Próprio autor.
124
Além disso, esse adiantamento do descolamento é um resultado contrário
daquilo que se esperava que o ressalto pudesse acarretar, que é a prorrogação do
estol e uma alteração benéfica nos coeficientes. Isso é evidenciado pelo fato de que,
além do ângulo de estol permanecer o mesmo (12º), o coeficiente de arrasto teve um
aumento nos valores de forma mais acentuada, e o coeficiente de sustentação uma
queda mais brusca.
Ao final, para um ângulo de ataque de 18º, o coeficiente de arrasto e de
sustentação são mais de 87% e quase 30% maior e menor, respectivamente. Isso é
Figuras 68. Campos de velocidade média para o ângulo de ataque igual à 10°(a) sem e (b) com o ressalto.
(a)
(b)
Fonte: Próprio autor.
125
nitidamente visto nas Figuras 69, Figuras 70 e Figuras 71. Além da intensificação dos
efeitos citados para o regime plenamente estabelecido (dissipação de energia do
fluido, antecipação do descolamento e ampliação da região de pressão adversa), as
Figuras 70 trazem uma representação direta dessa diferença dos coeficientes. Ao
provocar um ponto de estagnação, o ressalto aumenta os valores da pressão
envolvida na porção superior do perfil, e reforça a tese de que a pressão resultante de
sustentação sofre ônus direto.
Figuras 69. Campos de energia cinética turbulenta para o ângulo de ataque igual à 18°(a) sem e (b) com o ressalto.
(a)
(b)
Fonte: Próprio autor.
126
Figuras 70. Campos de pressão média para o ângulo de ataque igual à 18°(a) sem e (b) com o ressalto.
(a)
(b)
Fonte: Próprio autor.
127
Figuras 71. Campos de velocidade média para o ângulo de ataque igual à 18°(a) sem e (b) com o ressalto.
(a)
(b)
Fonte: Próprio autor.
128
E ao que tudo indica, tais consequências se estendem às configurações de
solidez. Os resultados segundo tal parâmetro, para um perfil com ressalto, são
apresentados nas Figura 72, Figura 73 e Figura 74, respectivamente os resultados
dos coeficientes de arrasto e sustentação, bem como na Tabela 18, os coeficientes
de pressão mínimo.
Figura 72. Variação do coeficiente de arrasto para o perfil com ressalto, segundo ângulo de
ataque para os valores de solidez entre 0 e 1.
Fonte: Próprio autor.
129
Figura 73. Variação do coeficiente de sustentação para o perfil com ressalto, segundo ângulo de ataque para os valores de solidez entre 0 e 1.
Fonte: Próprio autor.
Figura 74. Variação da razão entre o coeficiente de sustentação e arrasto para o perfil com ressalto, segundo ângulo de ataque para os valores de solidez entre 0 e 1.
Fonte: Próprio autor.
130
Tabela 18 Resultado dos coeficientes de pressão mínimo para o perfil com ressalto, segundo
ângulo ataque para os valores de solidez entre 0 e 1. Solidez 0° 1,25° 2,50° 3,75° 5° 10°
0,0 2,0230 2,2540 2,4630 2,6540 2,8320 3,6110
0,1 1,8320 2,0160 2,5390 2,7790 2,9810 3,8550
0,2 2,0860 2,2830 2,9400 3,1560 3,3880 4,4760
0,3 1,8060 2,0020 2,5850 2,8270 3,0650 4,1430
0,4 1,8170 2,0280 2,2360 2,9130 3,1700 4,2810
0,5 1,8020 2,0320 2,2600 2,9560 3,1950 4,4610
0,6 2,1020 2,3600 2,6060 3,5540 3,9080 5,2890
0,7 2,1230 2,4000 2,6700 3,7160 4,0660 5,6970
0,8 1,9660 2,2420 2,5840 2,9160 3,8830 5,5180
0,9 2,1950 2,7645 3,3340 3,9035 4,4730 6,2170
1,0 1,9870 2,2538 2,5205 2,7873 3,0540 5,5200 Fonte: Próprio autor.
As Figura 72, Figura 73 e Figura 74 reforçam a dependência direta dos
coeficientes aero/hidrodinâmicos com relação à solidez que as Figura 52, Figura 53 e
Figura 54 haviam apresentado. E ainda, mesmo com a presença do ressalto, o
comportamento dos valores em sua maioria se replicou, em que enquanto os números
da sustentação decrescem, os de arrasto aumentam.
Entretanto, notase que esses números apresentaram diferença com relação
aos resultados das Figura 52, Figura 53 e Figura 54, algo semelhante do ocorrido para
o perfil isolado, porém neste caso de forma menos acentuada. Isso se dá pela soma
dos efeitos da solidez e do ressalto: enquanto o ressalto provoca uma queda na
velocidade do escoamento, e consequentemente uma antecipação do descolamento
da camada limite e uma ampliação da região de pressão adversa, o escoamento da
solidez atenua isso pela compressão do fluido, que acarreta no aumento da sua
velocidade e não só induz a sua camada limite ser recolada, como também aumenta
as pressões envolvidas na porção superior do perfil. Tal hipótese pode ser confirmada
ao visualizar as Figuras 75 e Figuras 76, em que se constata que é praticamente
imperceptível a distinção dos campos com e sem ressalto.
Não à toa que, a partir do momento em que o impacto da solidez passa a ser
irrelevante (ou seja, em valores baixos), a implicação do ressalto se sobressai e os
coeficientes tendem aos resultados obtidos pelo perfil isolado. Por fim, é observado
131
uma queda significativa dos valores de coeficiente de pressão mínima, se comparado
com esses valores para os perfis sem ressalto (Tabela 18). A consequência direta é a
criação de zonas de pressão ainda menores, que torna mais suscetível a ocorrência
do fenômeno da cavitação.
Figuras 75. Campos de energia cinética turbulenta para o ângulo de ataque igual à 10° e uma solidez de 1,0 de perfis (a) sem e (b) com o ressalto.
(a)
(b)
Fonte: Próprio autor.
132
4.2.3.2 Algoritmo
Naturalmente, a tendência inversamente proporcional entre a razão ótima dos
coeficientes de sustentação e arrasto e a solidez para os perfis com ressalto se repete,
como pode ser visto na Figura 77. Mas destacase através dela, que à exceção dessa
razão em baixa solidez (entre 0 e 0,1), tanto os valores, como o comportamento deles
Figuras 76. Campos de velocidade média para o ângulo de ataque igual à 10° e uma solidez de 1,0 de perfis (a) sem e (b) com o ressalto.
(a)
(b)
Fonte: Próprio autor.
133
se mantiveram com relação aos perfis sem ressalto. Para que isso seja possível, a
alteração – mesmo que ínfima reside no ângulo de ataque ótimo.
Figura 77. Variação da razão ótima entre os coeficientes de sustentação e arrasto para o perfil
com ressalto e seus ângulos de ataque, para os valores de solidez entre 0 e 1.
Fonte: Próprio autor.
O procedimento de interpolação dos dados e a inserção deles como correção
dos coeficientes de escoamento é o mesmo da Seção 4.2.2.2. Novamente, analisase
o impacto da consideração dos coeficientes do escoamento como função da solidez,
mas alterado pela inserção do ressalto no perfil da pá. Para isso, um perfil de pá é
obtido com os mesmos parâmetros de design ótimo encontrados na Seção 4.2.1.1, e
a comparação da distribuição da corda e da torção podem ser vistos nas Figura 78 e
Figura 79, respectivamente.
134
Figura 78. Comparação da distribuição de corda com e sem ressalto.
Fonte: Próprio autor.
Figura 79. Comparação da distribuição de torção com e ressalto.
Fonte: Próprio autor.
135
A mudança na estrutura da pá é o indício mais claro do impacto que a escolha
dos perfis pode acarretar. Ainda, neste caso, a adoção de um perfil com ressalto
agrava o efeito de solidez da pá, como a diminuição dos esforços normais (Equação
(13)), que acarreta um aumento da corda (segundo Equação (59)) e da inércia. Este
parâmetro, inclusive, é a evidência de ambas as implicações: se tomarmos o design
ótimo obtido por meio dos coeficientes do software XFOIL (DRELA, 2001), o aumento
é de 177,61%, enquanto que para este mesmo layout, mas considerando o efeito de
solidez da Seção 4.2.2.2, esse número é 123,21%. Ademais, há uma queda de quase
23% na geração de potência, que significa um déficit de energia elétrica ao final do
mês em torno de 41,46 kWh. Contra isso, um novo processo de otimização deve ser
feito com o propósito de se minimizar tal diferença.
A otimização é feita novamente segundo o método validado de acordo com a
Tabela 8, e o processo iterativo é representado na Figura 80. Uma comparação entre
os resultados obtidos nas otimizações anteriores pode ser visto na Tabela 19, cuja
diferença percentual é com relação aos valores da Seção 4.2.1.1 considerando os
coeficientes provenientes do software XFOIL (DRELA, 2001).
Figura 80. Resultado da otimização do modelo híbrido ABCSA para um perfil com ressalto, e
também corrigido segundo resultados de solidez.
Fonte: Próprio autor.
136
Tabela 19 Comparação dos resultados de design ótimos para turbina de perfil: com e sem
ressalto, e também com e sem correção de solidez.
Resultado XFOIL DFC DFC (corrigido) DFC
(com ressalto) Diferença (%)
Coeficiente de potência 0,475509 0,448012 0,416217 0,363039 23,65
Inércia total (x105 kg*m2) 4,175588 4,692843 5,848779 16,894000 304,59
Diâmetro do rotor (m) 0, 550777 0,567617 0,588675 0,630421 14,46
Diâmetro do cubo (m) 0,13770 0,14191 0,14717 0,157612 14,46
Número de pás 3 3 3 3 0,00
Velocidaderotativa (rpm) 370,5412 370,5412 370,5412 370,5412 0,00
Λ 6,4846 6,6829 6,9308 7,4223 14,46
Potência (W) 252,722 252,891 252,699 252,782 0,02
Tempo de simulação (s) 539,80 496,83 692,98 834,11 54,52 Fonte: Próprio autor.
A Figura 80 representa previamente como o design otimizado levando em
conta os resultados hidrodinâmicos das simulações com ressalto é inferior em termos
estruturais e de eficiência na conversão de potência, já que o valor objetivo final saiu
de 8,30785 x 104 (Tabela 8) para 16,9286 x 104, um aumento de quase 104%. Através
da Tabela 19, isso é confirmado pela queda de 23,65% do coeficiente de potência,
mas principalmente impulsionado pelo aumento da inércia de 304,59%. Algo que não
se descarta ser a razão desse número é o diâmetro do rotor ser mais de 14% maior,
que pode ocasionar inclusive num custo maior da sua fabricação.
As Figura 81, Figura 82 e Figura 83 a seguir trazem as distribuições no
comprimento da pá de corda, torção e solidez para um design ótimo de um perfil com
ressalto e corrigido também pela solidez.
137
Figura 81. Distribuição da corda ao longo do comprimento da pá, segundo coeficientes de perfil com ressalto e corrigido pelos efeitos da solidez.
Fonte: Próprio autor.
Figura 82. Distribuição da torção ao longo do comprimento da pá, segundo coeficientes de perfil com ressalto e corrigido pelos efeitos da solidez.
Fonte: Próprio autor.
138
Figura 83. Distribuição da solidez ao longo do comprimento da pá, segundo coeficientes de
perfil com ressalto e corrigido pelos efeitos da solidez.
Fonte: Próprio autor.
Se observada a Figura 81, levase a crer que a disparidade da inércia pode
estar correlacionada muito por conta do prejuízo nos coeficientes hidrodinâmicos de
arrasto e sustentação, mas também com o fenômeno da cavitação. A partir de um raio
do rotor em torno de 0,25 m, há um aumento repentino da corda devido às correções
impostas para se evitar a cavitação, e que neste caso estaria mais suscetível de
ocorrer devido aos valores menores de coeficiente de pressão mínimo do perfil com
ressalto (conforme sugerido na Seção 4.2.3.1, e vide a diferença entre as Tabela 16
e Tabela 18). Se essa mudança é desabilitada, a distribuição da corda é alterada, e
pode ser comparada com a corrigida na Figura 84. A inércia passa a ser 13,44%
menor, mas por outro lado, podese dizer que a cavitação causa um déficit na geração
de potência de 3,53%.
139
Figura 84. Comparação da distribuição da corda ao longo do comprimento da pá com e sem correção que evita a cavitação, segundo coeficientes de perfil com ressalto e corrigido pelos
efeitos da solidez.
Fonte: Próprio autor.
Por mais que as respostas que o ressalto proporcionou não corresponderam
às expectativas iniciais de benefício ao escoamento e aos coeficientes
aero/hidrodinâmicos, e como consequência não teve a otimização do layout da turbina
em favor da geração de energia elétrica e de sua estrutura, uma nova interpretação
da sua presença pode ser levantada em questão. Se essa protuberância não fosse
algo intencional no perfil – como foi a tentativa deste trabalho , mas ao invés disso
uma incrustação presente na pá, os prejuízos seriam os mesmos. Ou seja, as
simulações numéricas e computacionais de DFC realizadas levam a crer que a falta
de manutenção das turbinas hidráulicas, ou mesmo as consequências naturais de
suas operações, torna a cavitação mais suscetível de ocorrer e prejudica diretamente
a geração de potência. Logo, a falta de um planejamento de manutenção adequado
pode não só diminuir a vida útil das turbinas hidráulicas, como também debilitar a
produção de energia elétrica.
140
5 CONCLUSÃO
A partir do presente trabalho, em que se buscou desenvolver um projeto
hidrodinâmico de turbina hidrocinética, partindo de uma metodologia baseada na
previsão de seu funcionamento e na otimização de seus parâmetros de layout, é
possível concluir que:
• Não é novidade a dificuldade em se prever o desempenho de uma
turbina. Devese salientar a escassez de dados experimentais
robustos e da sua detalhada descrição, o que dificulta uma
investigação mais profunda e um englobamento das implementações
das equações. Além disso, as teorias ao longo da história carecem de
uniformidade da nomenclatura, das adoções de suposições e do
esclarecimento de suas implementações. Mesmo assim, fora possível
desenvolver uma metodologia de validação inédita na literatura, e que
alcançou um resultado surpreendente: a validação experimental do
coeficiente de potência com 1,40% de erro. Este valor, inclusive,
surpreende quando apresenta desempenho mais expressivo do que o
resultado obtido por um software open source e consolidado no ramo;
• Algoritmos de otimização, como os adotados neste trabalho, fazem
parte de uma área antigamente considerada destoante da engenharia:
a inteligência artificial. A aplicação dela é promissora, e aqui fora um
exemplo disso: a partir das restrições de projeto, dentre as infinitas
possibilidades de construção de uma turbina hidrocinética, os quatro
métodos foram capazes de encontrar uma combinação dos parâmetros
que balancearia não só a sua eficiência, mas também a sua massa. O
coeficiente de potência da turbina otimizada fora igual a 0,475509, ou
seja, 18% a menos que o Limite de Betz;
• O que era previamente esperado, se confirmou: simulações DFC
complexas, como quando se tratam de turbinas, são peculiares e
exigem não somente tempo hábil, mas recursos computacionais
abundantes. Mesmo com um cluster e os seus 36 núcleos de
processamento, o custo computacional de simular os mais de 3 milhões
de elementos foi alto e não fora possível obter resultados consolidados
141
do estudo de independência de passo de tempo. E consequentemente,
impossibilitou a análise de independência de malha. Deste modo,
recomendase como trabalho futuro um enfoque exclusivo em
desenvolver uma simulação como esta;
• É inegável o efeito da solidez em escoamentos. Ao longo dos anos,
algumas tratativas foram abordadas e serviram de inspiração para a
hipótese aqui levantada. A correção dos coeficientes
aero/hidrodinâmicos em função desse efeito levou à um decréscimo de
8,15% da potência gerada pela turbina otimizada. Recomendase como
trabalhos futuros, uma investigação mais aprofundada na proposição,
através da simulação DFC do design de uma turbina corrigida pelos
efeitos da solidez, ou mesmo a prototipagem e validação experimental;
• As protuberâncias, enquanto artifícios de indução de turbulência e
instrumentos de melhora do arrasto e da sustentação, não são
novidades no âmbito de estudo fluidodinâmico. Entretanto, o ressalto
adotado neste trabalho se mostrou ineficaz neste objetivo. Pelo
contrário, se mostrou prejudicial aos coeficientes aero/hidrodinâmicos.
E consequentemente, houve uma queda na potência gerada de quase
23%, e um aumento de 177,61% na inércia da turbina. A partir disso,
recomendase em trabalhos futuros uma investigação em torno da
alteração deste ressalto em: formato, tamanho e posicionamento no
perfil aero/hidrodinâmico;
• Recomendase também prosseguir com a análise de otimização com
adição dos difusores, que fora impossibilitado nesse trabalho por conta
do cronograma.
A metodologia e o trabalho abrem margens para inúmeras novas análises, que
vão desde à mudança do fluido (vento dando lugar à água) até a outros perfis
aero/hidrodinâmicos. Assim sendo, é notável a aplicabilidade que a ferramenta
desenvolvida apresenta. Unir tradicionais engenhos mecânicos, com as ferramentas
de otimização impulsionadas pela inteligência artificial, acrescenta uma nova
abordagem não só à sustentabilidade, mas à engenharia clássica. Isso extrapola a
esfera acadêmica, e lança um olhar sobre como se pode aliar a modernidade ao
habitual e consolidado, em prol de uma melhora no aproveitamento dos recursos
naturais, ou mesmo da qualidade de vida.
142
REFERÊNCIAS
ABUAN, B. E.; HOWELL, R. J. The performance and hydrodynamics in unsteady flow of a horizontal axis tidal turbine. Renewable Energy, v. 133, p. 1338–1351, 1 abr. 2019.
ADVANCED SIMULATION AND COMPUTING INITIATIVE (ASCI). VisIt. [s.l.] Lawrence Livermore National Laboratory, 2002.
AHMED, N.; YILBAS, B. S.; BUDAIR, M. O. Computational Study into the Flow Field Developed around a Cascade of NACA 0012 Airfoils. Computer Methods in Applied Mechanics and Engineering, v. 167, n. 1, p. 17–32, 1 dez. 1998.
AMARANTE MESQUITA, A. L. et al. A Methodology for the Transient Behavior of Horizontal Axis Hydrokinetic Turbines. Energy Conversion and Management, v. 87, p. 1261–1268, 1 nov. 2014.
ANA. Quantidade da água. Disponível em: <https://www.ana.gov.br/panoramadasaguas/quantidadedaagua>. Acesso em: 21 set. 2019.
ANEEL. Capacidade de Geração do Brasil. Disponível em: <https://www2.aneel.gov.br/aplicacoes/capacidadebrasil/capacidadebrasil.cfm>. Acesso em: 4 nov. 2019.
ARAÚJO, M. A. de. Prospecção de Parques Hidrocinéticos: Comparação entre Projetos Preliminares no Rio Iguaçu e Paraná. 2016. Universidade Federal da Integração LatinoAmericana, Foz do Iguaçu, 2016. Disponível em: <https://dspace.unila.edu.br/bitstream/handle/123456789/684/TCC%20%20Marcos%20Aurelio%20de%20Araujo.pdf?sequence=1&isAllowed=y>.
ARDIZZON, G.; CAVAZZINI, G.; PAVESI, G. A new generation of small hydro and pumpedhydro power plants: Advances and future challenges. Renewable and Sustainable Energy Reviews, v. 31, p. 746–761, 1 mar. 2014.
ARGYROPOULOS, C. D.; MARKATOS, N. C. Recent Advances on the Numerical Modelling of Turbulent Flows. Applied Mathematical Modelling, v. 39, n. 2, p. 693–732, 15 jan. 2015.
ARNDT, R. E. A. Cavitation in Fluid Machinery and Hydraulic Structures. Annual Review of Fluid Mechanics, v. 13, n. 1, p. 273–326, 1981.
ASNAGHI, A.; SVENNBERG, U.; BENSOW, R. E. Evaluation of Curvature Correction Methods for Tip Vortex Prediction in SST K−ω Turbulence Model Framework. International Journal of Heat and Fluid Flow, v. 75, p. 135–152, 1 fev. 2019.
AVELLAN, F. Introduction to Cavitation in Hydraulic Machinery. Disponível em: <https://infoscience.epfl.ch/record/254965>. Acesso em: 17 out. 2019.
BEKDAŞ, G.; NIGDELI, S. M.; YANG, X.S. Sizing Optimization of Truss Structures Using Flower Pollination Algorithm. Applied Soft Computing, v. 37, p. 322–331, 1 dez. 2015.
143
BETZ, A. Das Maximum der theoretisch moeglichen Ausnutzung des Windes durch Windmotoren. Zeitschrift fur das gesamte Turbinenwesten, v. 20, 1920. Disponível em: <https://ci.nii.ac.jp/naid/10005373900/>. Acesso em: 10 out. 2019.
BLACK & VEATCH. Phase II UK Tidal Stream Energy Resource Assessment. [s.l.] The Carbon Trust, 2005. . Disponível em: <https://www.carbontrust.com/media/174041/phaseiitidalstreamresourcereport2005.pdf>.
BORHANAZAD, H. et al. Optimization of MicroGrid System Using MOPSO. Renewable Energy, v. 71, p. 295–306, 1 nov. 2014.
BRACKEN, L. J.; BULKELEY, H. A.; MAYNARD, C. M. Microhydro power in the UK: The role of communities in an emerging energy resource. Energy Policy, v. 68, p. 92–101, 1 maio 2014.
BRANLARD, E. Wind Turbine Aerodynamics and VorticityBased Methods: Fundamentals and Recent Applications. [s.l.] Springer International Publishing, 2017.
BRASIL JUNIOR, A. C. P. et al. On the Design of Propeller Hydrokinetic Turbines: The Effect of the Number of Blades. Journal of the Brazilian Society of Mechanical Sciences and Engineering, v. 41, n. 6, p. 253, 16 maio 2019.
BREKKE, H. Hydraulic Turbines. Design, Erection and Operation. [s.l: s.n.]
BUHL, M. L. New Empirical Relationship between Thrust Coefficient and Induction Factor for the Turbulent Windmill State. [s.l.] National Renewable Energy Lab. (NREL), Golden, CO (United States), 1 ago. 2005. . Disponível em: <https://www.osti.gov/biblio/15016819newempiricalrelationshipbetweenthrustcoefficientinductionfactorturbulentwindmillstate>. Acesso em: 12 out. 2019.
BURTON, T. et al. Wind Energy Handbook. [s.l.] John Wiley & Sons, 2011.
CARVALHO, E. Considerada fracasso na época, Rio 92 foi “sucesso” para especialistas. Disponível em: <http://g1.globo.com/natureza/rio20/noticia/2012/05/consideradafracassonaepocario92foisucessoparaespecialistas.html>. Acesso em: 9 out. 2019.
CEBRIÁN, D.; ORTEGACASANOVA, J.; FERNANDEZFERIA, R. Lift and drag characteristics of a cascade of flat plates in a configuration of interest for a tidal current energy converter: Numerical simulations analysis. Journal of Renewable and Sustainable Energy, v. 5, n. 4, p. 043114, 1 jul. 2013.
CLIFTONSMITH, M. J. Wind Turbine Blade Optimisation with Tip Loss Corrections. Wind Engineering, v. 33, n. 5, p. 477–496, 1 out. 2009.
CORRÊA DA SILVA, R.; DE MARCHI NETO, I.; SILVA SEIFERT, S. Electricity Supply Security and the Future Role of Renewable Energy Sources in Brazil. Renewable and Sustainable Energy Reviews, v. 59, p. 328–341, 1 jun. 2016.
144
DASKIRAN, C. et al. Transient Analysis of MicroHydrokinetic Turbines for River Applications. Ocean Engineering, v. 129, p. 291–300, 1 jan. 2017.
DASSAULT SYSTÈMES. SOLIDWORKS. San Diego: BIOVIA, [s.d.]
DEKKERS, A.; AARTS, E. Global Optimization and Simulated Annealing. Mathematical Programming, v. 50, n. 1, p. 367–393, 1 mar. 1991.
DIXON, S. L. Fluid mechanics, thermodynamics of turbomachinery. Amsterdam; Boston: ElsevierButterworthHeinemann, 2005.
DJAVARESHKIAN, M. H.; ESMAEILI, A. Heuristic Optimization of Submerged Hydrofoil Using ANFIS–PSO. Ocean Engineering, v. 92, p. 55–63, 1 dez. 2014.
DO RIO VAZ, D. A. T. D.; VAZ, J. R. P.; SILVA, P. A. S. F. An Approach for the Optimization of DiffuserAugmented Hydrokinetic Blades Free of Cavitation. Energy for Sustainable Development, v. 45, p. 142–149, 1 ago. 2018.
DRELA, M. XFOIL. [s.l.] MIT, 2001.
DRZEWIECKI, S.; GAUTHIERVILLARS (PARYŻ). Méthode pour la détermination des éléments mécaniques des propulseurs hélicoïdaux. Paris: imprimerie GauthierVillars et Fils, 1892.
DUQUETTE, M. M.; SWANSON, J.; VISSER, K. D. Solidity and Blade Number Effects on a Fixed Pitch, 50 W Horizontal Axis Wind Turbine. Wind Engineering, v. 27, n. 4, p. 299–316, 1 ago. 2003.
ELETROBRAS. Principais Linhas de Transmissão Sistema Eletrobras 2018, dez. 2018.
EPE; ONS. Projeção de demanda de energia elétrica para os próximos 10 anos (20172026): Estudos da Demanda. Rio de Janeiro: Ministério de Minas e Energia MME, jan. 2017. . Disponível em: <http://www.epe.gov.br/sitespt/publicacoesdadosabertos/publicacoes/PublicacoesArquivos/publicacao245/topico261/DEA%20001_2017%20%20Proje%C3%A7%C3%B5es%20da%20Demanda%20de%20Energia%20El%C3%A9trica%2020172026_VF[1].pdf>.
FERZIGER, J. H.; PERIC, M. Computational Methods for Fluid Dynamics. 3. ed. Berlin Heidelberg: SpringerVerlag, 2002.
FETANAT, A.; KHORASANINEJAD, E. Size Optimization for Hybrid Photovoltaic–Wind Energy System Using Ant Colony Optimization for Continuous Domains Based Integer Programming. Applied Soft Computing, v. 31, p. 196–209, 1 jun. 2015.
FLECK, G. D. Numerical analysis of the solidity effects over the aerodynamic performance of a small wind turbine. 2017. Universidade Federal do Rio Grande do Sul, 2017. Disponível em: <https://lume.ufrgs.br/handle/10183/173195>. Acesso em: 29 mar. 2020.
145
FROUDE, R. E. On the Part Played in Propulsion by Differences of Fluid Pressure. Transactions of the Institution of Naval Architects, v. 30, p. 16, 1889.
FROUDE, W. On the elementary relation between pitch, slip and propulsive efficiency. London: Institution of Naval Architects, 1878.
GADEN, D. L. F.; BIBEAU, E. L. A Numerical Investigation into the Effect of Diffusers on the Performance of Hydro Kinetic Turbines Using a Validated Momentum Source Turbine Model. Renewable Energy, v. 35, n. 6, p. 1152–1158, 1 jun. 2010.
GLAUERT, H. Airplane Propellers. In: DURAND, W. F. (Ed.). Aerodynamic Theory: A General Review of Progress Under a Grant of the Guggenheim Fund for the Promotion of Aeronautics. Berlin, Heidelberg: Springer Berlin Heidelberg, 1935. p. 169–360.
GOLDSTEIN, S.; PRANDTL, L. On the vortex theory of screw propellers. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, v. 123, n. 792, p. 440–465, 6 abr. 1929.
HALDER, P.; RHEE, S. H.; SAMAD, A. Numerical optimization of Wells turbine for wave energy extraction. International Journal of Naval Architecture and Ocean Engineering, v. 9, n. 1, p. 11–24, 1 jan. 2017.
HANSEN, M. O. L. Aerodynamics of Wind Turbines. [s.l.] Earthscan, 2013.
HANSEN, M. O. L.; SØRENSEN, N. N.; FLAY, R. G. J. Effect of Placing a Diffuser around a Wind Turbine. Wind Energy, v. 3, n. 4, p. 207–213, 2000.
HART, C.; SALING, S. To cut poverty in Asia and the Pacific, ‘Energy Plus’ package a must, says UN report. Disponível em: <https://www.undp.org/content/undp/en/home/presscenter/pressreleases/2012/01/19/tocutpovertyinasiaandthepacificenergypluspackageamustsaysunreport.html>. Acesso em: 9 out. 2019.
HASANÇEBI, O.; ÇARBAŞ, S.; SAKA, M. P. Improving the Performance of Simulated Annealing in Structural Optimization. Structural and Multidisciplinary Optimization, v. 41, n. 2, p. 189–203, 1 mar. 2010.
HEINLOTH, K. (ed.). Renewable Energy. [s.l.] Springer Berlin Heidelberg, 2006. v. 3C
HOGHOOGHI, H.; DURALI, M.; KASHEF, A. A new lowcost swirler for axial micro hydro turbines of low head potential. Renewable Energy, v. 128, p. 375–390, 1 dez. 2018.
IEA. Key World Energy Statistics 2018. [s.l.] International Energy Association, 19 set. 2018. . Disponível em: <https://webstore.iea.org/keyworldenergystatistics2018>. Acesso em: 9 out. 2019.
IRENA. Hydropower: Renewable Energy Technologies: Cost Analysis Series. [s.l.] International Renewable Energy Agency, jun. 2012. . Disponível em:
146
<https://www.irena.org/documentdownloads/publications/re_technologies_cost_analysishydropower.pdf>.
ISLAM, Q. Md.; ISLAM, S. A. k. m. The Aerodynamic Performance of a HorizontalAxis Wind Turbine Calculated by Strip Theory and Cascade Theory. JSME International Journal Series B, v. 37, n. 4, p. 871–877, 1994.
JACOBSON, P. Assessment and Mapping of the Riverine Hydrokinetic Resource in the Continental United States. [s.l: s.n.]. Disponível em: <http://www.osti.gov/servlets/purl/1219876/>. Acesso em: 9 out. 2019.
JAVADI, A.; NILSSON, H. Detailed Numerical Investigation of a Kaplan Turbine with RotorStator Interaction Using TurbulenceResolving Simulations. International Journal of Heat and Fluid Flow, v. 63, p. 1–13, 1 fev. 2017.
JONKMAN, J. M. Modeling of the UAE Wind Turbine for Refinement of FAST_AD. [s.l.] National Renewable Energy Lab., Golden, CO (US), 1 dez. 2003. . Disponível em: <https://www.osti.gov/biblio/15005920modelinguaewindturbinerefinementfastad>. Acesso em: 12 out. 2019.
KAMJOO, A. et al. MultiObjective Design under Uncertainties of Hybrid Renewable Energy System Using NSGAII and Chance Constrained Programming. International Journal of Electrical Power & Energy Systems, v. 74, p. 187–194, 1 jan. 2016.
KARABOGA, D. AN IDEA BASED ON HONEY BEE SWARM FOR NUMERICAL OPTIMIZATION. In: Anais...2005.
KAUNDA, C. S.; KIMAMBO, C. Z.; NIELSEN, T. K. A technical discussion on microhydropower technology and its turbines. Renewable and Sustainable Energy Reviews, v. 35, p. 445–459, 1 jul. 2014.
KENNEDY, J.; EBERHART, R. Particle swarm optimization. In: Proceedings of ICNN’95 International Conference on Neural Networks, Anais... In: PROCEEDINGS OF ICNN’95 INTERNATIONAL CONFERENCE ON NEURAL NETWORKS. nov. 1995.
KIRKPATRICK, S.; GELATT, C. D.; VECCHI, M. P. Optimization by Simulated Annealing. Science, v. 220, n. 4598, p. 671–680, 1983.
KOLMOGOROV, A. N. Equations of turbulent motion in an incompressible fluid. [s.l: s.n.]
KUMAR, P.; SAINI, R. P. Study of cavitation in hydro turbines—A review. Renewable and Sustainable Energy Reviews, v. 14, n. 1, p. 374–383, 1 jan. 2010.
LABIGALINI, L. et al. LIFT AND DRAG CORRELATIONS FOR AN AIRFOIL ARRAY THROUGH NUMERICAL SIMULATION ANALYSIS. In: Proceedings of the 25th International Congress of Mechanical Engineering, Anais... In: 25TH INTERNATIONAL CONGRESS OF MECHANICAL ENGINEERING. ABCM, 2019. Disponível em: <http://abcm.org.br/anaisdeeventos/COB2019/1872>. Acesso em: 30 mar. 2020.
147
LANZAFAME, R.; MESSINA, M. Advanced brake state model and aerodynamic poststall model for horizontal axis wind turbines. Renewable Energy, v. 50, p. 415–420, 1 fev. 2013.
LAUNDER, B. E.; SHARMA, B. I. Application of the EnergyDissipation Model of Turbulence to the Calculation of Flow near a Spinning Disc. Letters Heat Mass Transfer, v. 1, p. 131, dez. 1974.
LEHR, J. H.; KEELEY, J. W.; KINGERY, T. B. Alternative energy and shale gas encyclopedia. [s.l: s.n.]
LETCHER, T. M. Wind energy engineering: a handbook for onshore and offshore wind turbines. [s.l: s.n.]
LIMA MARTINS BARBOSA, D. et al. An Investigation of a Mathematical Model for the Internal Velocity Profile of Conical Diffusers Applied to DAWTs. Anais da Academia Brasileira de Ciencias, v. 87, 28 abr. 2015.
MANIACI, D. C.; KELLEY, C. L.; CHIU, P. Assessment of Scaled Rotors for Wind Tunnel Experiments. [s.l.] Sandia National Lab. (SNLNM), Albuquerque, NM (United States); Sandia National Laboratories, Los Angeles, CA, 1 jul. 2015. . Disponível em: <https://www.osti.gov/biblio/1235645>. Acesso em: 17 mar. 2020.
MANWELL, J. F.; MCGOWAN, J. G.; ROGERS, A. L. Wind energy explained: theory, design and application. Chichester, U.K.: John Wiley & Sons, Ltd., 2011.
MARK DRELA et al. Area and Bending Inertia of Airfoil SectionsMIT OpenCourseWare, , 2006. . Disponível em: <https://ocw.mit.edu/courses/aeronauticsandastronautics/1601unifiedengineeringiiiiiiivfall2005spring2006/systemslabs06/spl10b.pdf>. Acesso em: 23 mar. 2020.
MARTEN, D. et al. Integration of a WT Blade Design Tool in XFoil/XFLR5. In: Anais...1 jan. 2010.
MARTEN, D. et al. QBLADE: An Open Source Tool for Design and Simulation of Horizontal and Vertical Axis Wind Turbines. International Journal of Emerging Technology and Advanced Engineering, v. 3, p. 264–269, 1 mar. 2013.
MARTÍNEZ, J. et al. An Improved BEM Model for the Power Curve Prediction of StallRegulated Wind Turbines. Wind Energy, v. 8, n. 4, p. 385–402, 2005.
MATYUSHENKO, A. A.; KOTOV, E. V.; GARBARUK, A. V. Calculations of Flow around Airfoils Using TwoDimensional RANS: An Analysis of the Reduction in Accuracy. St. Petersburg Polytechnical University Journal: Physics and Mathematics, v. 3, n. 1, p. 15–21, 1 mar. 2017.
Meio ambiente. Disponível em: <https://nacoesunidas.org/acao/meioambiente/>. Acesso em: 9 out. 2019.
MENNY, K. Strömungsmaschinen: Hydraulische und thermische Kraft und Arbeitsmaschinen. 5. ed. [s.l.] Vieweg+Teubner Verlag, 2006.
148
MENTER, F. R. TwoEquation EddyViscosity Turbulence Models for Engineering Applications. AIAA journal., v. 32, n. 8, p. 1598–1605, 1994.
MICLEABLEIZIFFER, M.; UNTAROIU, A.; DELGADO, A. Development of a novel design method for marine propellers by computing the exact lift of arbitrary hydrofoils in cascades. Ocean Engineering, v. 83, p. 87–98, 1 jun. 2014.
MORIARTY, P. J.; HANSEN, A. C. AeroDyn Theory Manual. [s.l.] National Renewable Energy Lab., Golden, CO (US), 1 jan. 2005. . Disponível em: <https://www.osti.gov/biblio/15014831>. Acesso em: 10 out. 2019.
NAGY, D. Design optimization. Disponível em: <https://medium.com/generativedesign/designoptimization2ec2ba3b40f7>. Acesso em: 3 nov. 2019.
NUERNBERG, M.; TAO, L. Three Dimensional Tidal Turbine Array Simulations Using OpenFOAM with Dynamic Mesh. Ocean Engineering, v. 147, p. 629–646, 1 jan. 2018.
OBLAS, N. Design, Manufacture and Prototyping of a Hydrokinetic Turbine Unit for River Application. Theses and Dissertations, 1 jan. 2016. Disponível em: <https://preserve.lehigh.edu/etd/2747>.
OLIVEIRA, V. Sem fornecimento da Venezuela, custo para manter energia em RR chega a R$ 1,6 bilhão em um ano. Disponível em: <https://g1.globo.com/rr/roraima/noticia/2020/03/09/semfornecimentodavenezuelacustoparamanterenergiaemrrchegaar16bilhaoemumano.ghtml>. Acesso em: 14 out. 2020.
ONS. Sistema de Transmissao Horizonte 2024, [s.d.]
OPENFOAM. The Open Source Computational Fluid Dynamics (CFD) Toolbox. Disponível em: <http://www.openfoam.com>. Acesso em: 2 nov. 2019.
O’ROURKE, F.; BOYLE, F.; REYNOLDS, A. Tidal current energy resource assessment in Ireland: Current status and future update. Renewable and Sustainable Energy Reviews, v. 14, n. 9, p. 3206–3212, 1 dez. 2010.
PECHLIVANOGLOU, G. Passive and Active Flow Control Solutions for Wind Turbine Blades. 31 jan. 2013. Disponível em: <https://depositonce.tuberlin.de/handle/11303/3784>. Acesso em: 24 nov. 2019.
PENCHE, C.; ESHA. Guide on How to Develop a Small Hydropower Plant. Brussels: European Small Hydropower Association, 2004.
PRANDTL, L.; BETZ, A. Vier Abhandlungen zur Hydrodynamik und Aerodynamik: (Flüssigkeit mit kleiner Reibung; Tragflügeltheorie, I. und II. Mitteilung; Schraubenpropeller mit geringstem Energieverlust). [s.l.] Selbstverlag des Kaiser Wilhelm Instituts für Strömungsforschung, 1927.
PRATUMNOPHARAT, P.; LEUNG, P. S. Validation of various windmill brake state models used by blade element momentum calculation. Renewable Energy, v. 36, n. 11, p. 3222–3227, 1 nov. 2011.
149
QUADROS, T. O histórico dos principais encontros e acordos climáticos mundiais. Disponível em: <https://www.nexojornal.com.br/grafico/2017/11/17/Ohist%C3%B3ricodosprincipaisencontroseacordosclim%C3%A1ticosmundiais>. Acesso em: 9 out. 2019.
RANKINE, W. J. M. On the Mechanical Principles of the Action of Propellers. [s.l: s.n.]
RIGLIN, J. et al. Experimental and numerical characterization of a fullscale portable hydrokinetic turbine prototype for river applications. Renewable Energy, v. 99, p. 772–783, 1 dez. 2016.
RITCHIE, H. Energy Production and Consumption. Our World in Data, 2021. Disponível em: <https://ourworldindata.org/energyproductionconsumption>.
ROACHE, P. J. Verification and validation in computational science and engineering. Albuquerque, N.M.: Hermosa, 1998.
ROSER, M.; RITCHIE, H.; ORTIZOSPINA, E. World Population Growth. Our World in Data, maio 2019. Disponível em: <https://ourworldindata.org/worldpopulationgrowth>. Acesso em: 17 out. 2019.
SANTOS, I. F. S. dos. Análise técnica e econômica de parques hidrocinéticos com base em previsões numéricas (CFD) e dados experimentais. 2019. Universidade Federal de Itajubá, Itajubá, 2019. Disponível em: <http://repositorio.unifei.edu.br/xmlui/handle/123456789/2001>. Acesso em: 9 out. 2019.
SCHLEICHER, W. C.; RIGLIN, J. D.; OZTEKIN, A. Numerical characterization of a preliminary portable microhydrokinetic turbine rotor design. Renewable Energy, v. 76, p. 234–241, 1 abr. 2015.
SCHLICHTING, H.; KESTIN, J. Boundary layer theory. New York: McGrawHill, 1979.
SHEN, W. Z. et al. Tip Loss Corrections for Wind Turbine Computations. Wind Energy, v. 8, n. 4, p. 457–475, 2005.
SIERKSMA, G. Linear and integer programming: theory and practice. New York: Marcel Dekker, 2002.
SIIROLA, J. J. Speculations on global energy demand and supply going forward. Current Opinion in Chemical Engineering, Energy and environmental engineering / Reaction engineering. v. 5, p. 96–100, 1 ago. 2014.
SILVA, P. A. S. F. et al. Analysis of cavitation for the optimized design of hydrokinetic turbines using BEM. Applied Energy, Clean, Efficient and Affordable Energy for a Sustainable Future. v. 185, p. 1281–1291, 1 jan. 2017.
SINGH, K. Blade element analysis and experimental investigation of high solidity wind turbines. 2014. University of Calgary, 2014. Disponível em: <https://prism.ucalgary.ca/handle/11023/1774>. Acesso em: 29 mar. 2020.
150
SONG, K.; WANG, W.Q.; YAN, Y. Numerical and Experimental Analysis of a DiffuserAugmented MicroHydro Turbine. Ocean Engineering, v. 171, p. 590–602, 1 jan. 2019.
SPERA, D. Wind Turbine Technology. [s.l: s.n.]
SPERA, D. A. (ed.). Wind Turbine Aerodynamics Part A: Basic Principles. In: Wind Turbine Technology: Fundamental Concepts in Wind Turbine Engineering, Second Edition. [s.l.] ASME Press, 2009. p. 281–350.
TAVARES DIAS DO RIO VAZ, D. A. et al. Optimum aerodynamic design for wind turbine blade with a Rankine vortex wake. Renewable Energy, v. 55, p. 296–304, 1 jul. 2013.
TAVARES DIAS DO RIO VAZ, D. A. et al. An Extension of the Blade Element Momentum Method Applied to Diffuser Augmented Wind Turbines. Energy Conversion and Management, v. 87, p. 1116–1123, 1 nov. 2014.
TIAN, W. et al. Numerical Simulations of a Horizontal Axis Water Turbine Designed for Underwater Mooring Platforms. International Journal of Naval Architecture and Ocean Engineering, v. 8, n. 1, p. 73–82, 1 jan. 2016.
VAZ, J. R. P. et al. Drivetrain Resistance and Starting Performance of a Small Wind Turbine. Renewable Energy, v. 117, p. 509–519, 1 mar. 2018.
VAZ, J. R. P. et al. Powertrain Assessment of Wind and Hydrokinetic Turbines with Diffusers. Energy Conversion and Management, v. 195, p. 1012–1021, 1 set. 2019.
VAZ, J. R. P.; WOOD, D. H. Aerodynamic Optimization of the Blades of DiffuserAugmented Wind Turbines. Energy Conversion and Management, v. 123, p. 35–45, 1 set. 2016.
VAZ, J. R. P.; WOOD, D. H. Effect of the Diffuser Efficiency on Wind Turbine Performance. Renewable Energy, v. 126, p. 969–977, 1 out. 2018.
VERMAAK, H. J.; KUSAKANA, K.; KOKO, S. P. Status of microhydrokinetic river technology in rural applications: A review of literature. Renewable and Sustainable Energy Reviews, v. 29, p. 625–633, 1 jan. 2014.
VRIES, O. de. Fluid dynamic aspects of wind energy conversion. In: Anais...1979.
WADCOCK, A. J. Investigation of lowspeed turbulent separated flow around airfoils. [s.l.] NASA, 1 ago. 1987. . Disponível em: <https://ntrs.nasa.gov/search.jsp?R=19880003089>. Acesso em: 21 mar. 2020.
WANG, W.Q.; YIN, R.; YAN, Y. Design and prediction hydrodynamic performance of horizontal axis microhydrokinetic river turbine. Renewable Energy, v. 133, p. 91–102, 1 abr. 2019.
WILCOX, D. C. Reassessment of the scaledetermining equation for advanced turbulence models. AIAA Journal, v. 26, n. 11, p. 1299–1310, 1988.
151
WILLIAMS, A. A.; SIMPSON, R. Pico hydro – Reducing technical risks for rural electrification. Renewable Energy, 2007 World Renewable Energy Conference Pacific Rim Region. v. 34, n. 8, p. 1986–1991, 1 ago. 2009.
WILSON, R. E.; LISSAMAN, P. B. S. Applied aerodynamics of wind power machines. NASA STI/Recon Technical Report N, v. 75, 1 jul. 1974. Disponível em: <http://adsabs.harvard.edu/abs/1974STIN...7522669W>. Acesso em: 9 out. 2019.
WOOD, D. Aerofoils: Lift, Drag, and Circulation. In: WOOD, D. (Ed.). Small Wind Turbines: Analysis, Design, and Application. Green Energy and Technology. London: Springer London, 2011a. p. 57–75.
WOOD, D. Blade Element Theory for Wind Turbines. In: WOOD, D. (Ed.). Small Wind Turbines: Analysis, Design, and Application. Green Energy and Technology. London: Springer London, 2011b. p. 41–55.
WOOD, D. Blade Element Calculations. In: WOOD, D. (Ed.). Small Wind Turbines: Analysis, Design, and Application. Green Energy and Technology. London: Springer London, 2011c. p. 77–99.
WOOD, D. Blade Design, Manufacture, and Testing. In: WOOD, D. (Ed.). Small Wind Turbines: Analysis, Design, and Application. Green Energy and Technology. London: Springer London, 2011d. p. 119–143.
YANG, X.S. Flower Pollination Algorithm for Global Optimization. (J. DurandLose, N. Jonoska, Eds.) In: Unconventional Computation and Natural Computation, Berlin, Heidelberg. Anais... Berlin, Heidelberg: Springer, 2012.
YILDIZ, A. R. A New Hybrid Artificial Bee Colony Algorithm for Robust Optimal Design and Manufacturing. Applied Soft Computing, v. 13, n. 5, p. 2906–2912, 1 maio 2013.
ZHANG, W. et al. Optimization with a Simulated Annealing Algorithm of a Hybrid System for Renewable Energy Including Battery and Hydrogen Storage. Energy, v. 163, p. 191–207, 15 nov. 2018.