Embed Size (px)
Citation preview
SIMULAÇÃO NUMÉRICA DO ESCOAMENTO AO REDOR DE AEROFÓLIOS VIA
MÉTODO DE VÓRTICES ASSOCIADO AO MÉTODO DOS PAINÉIS
Daniel Fonseca de Carvalho e Silva
DISSERTAÇÃO SUBMETIDA AO CORPO DOCENTE DA CORDENAÇÃO DOS
PROGRAMAS DE PÓS-GRADUAÇÃO DE ENGENHARIA DA UNIVERSIDADE
FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS NECESSÁRIOS
PARA A OBTENÇÃO DO GRAU DE MESTRE EM CIÊNCIAS EM ENGENHARIA
MECÂNICA.
Aprovada por:
____________________________________________
Prof. Gustavo César Rachid Bodstein, Ph.D.
____________________________________________
Prof. Albino José Kalab Leiroz, Ph.D.
____________________________________________
Prof. Marcelo José Colaço. D.Sc.
____________________________________________
Prof. Luiz Antonio Alcântara Pereira, D. Sc.
RIO DE JANEIRO, RJ – BRASIL
SETEMBRO DE 2005
ii
SILVA, DANIEL FONSECA DE CARVALHO E
Simulação Numérica do Escoamento ao
redor de Aerofólios Via Método de Vórtices
Associado ao Método dos Painéis [Rio de
Janeiro] 2005
XIII, 175 p. 29,7 cm (COPPE/UFRJ,
M.Sc., Engenharia Mecânica, 2005)
Dissertação - Universidade Federal do
Rio de Janeiro, COPPE
1. Método de Vórtices
2. Método dos Painéis
3. Aerofólios
4. Aerodinâmica
I. COPPE/UFRJ II. Título ( série )
iii
Resumo da Dissertação apresentada à COPPE/UFRJ como parte dos requisitos necessários
para a obtenção do grau de Mestre em Ciências (M.Sc.)
SIMULAÇÃO NUMÉRICA DO ESCOAMENTO AO REDOR DE AEROFÓLIOS VIA
MÉTODO DE VÓRTICES ASSOCIADO AO MÉTODO DOS PAINÉIS
Daniel Fonseca de Carvalho e Silva
Setembro/2005
Orientador: Gustavo César Rachid Bodstein
Programa: Engenharia Mecânica
Este trabalho apresenta um estudo numérico de escoamentos bidimensionais,
incompressíveis, não-permanentes e de alto número de Reynolds ao redor de aerofólios. O
algoritmo numérico empregado associa o Método de Vórtices Discretos e o Método dos
Painéis para realizar simulações que cobrem uma ampla faixa de ângulos de ataque e
algumas variações no número de Reynolds, para dois tipos de aerofólios.
No modelo numérico desenvolvido utiliza-se o Método dos Painéis com
singularidades do tipo vórtices com distribuição linear para calcular a contribuição do
aerofólio sobre o escoamento. Para o Método de Vórtices Discretos, utiliza-se o Método do
Avanço Randômico para simular a difusão e a Lei de Biot-Savart para o cálculo da
velocidade convectiva. O esquema de Adams-Bashforth de segunda ordem é usado para o
avanço temporal. As distribuições de vorticidade na superfície do corpo são calculadas de
tal forma a satisfazer a condição de impenetrabilidade e a conservação de circulação. A
cada passo de tempo são gerados novos vórtices, que deslocam-se por ação dos efeitos de
difusão e convecção, formando a esteira à jusante do corpo.
Os resultados obtidos para diversos ângulos de ataque apresentam razoável
concordância com resultados experimentais, sobretudo na captura do fenômeno de “estol”
sobre aerofólios da família NACA. Estes resultados indicam a grande capacidade do
algoritmo e dos modelos físicos empregados neste trabalho para simulação numérica de
escoamentos complexos de altos números de Reynolds.
iv
Abstract of Dissertation presented to COPPE/UFRJ as a partial fulfillment of the
requirements for the degree of Master of Science (M.Sc.)
NUMERICAL SIMULATION OF THE FLOW AROUND AIRFOILS USING THE
DISCRETE VORTEX METHOD IN ASSOCIATION WITH THE PANEL METHOD
Daniel Fonseca de Carvalho e Silva
September/2005
Advisor: Gustavo César Rachid Bodstein
Department: Mechanical Engineering
This work presents a numerical study of two-dimensional, incompressible,
unsteady and high Reynolds number flows around airfoils. A numerical algorithm, that
uses the Discrete Vortex Method in association with the Panel Method, is developed to
carry out simulations for two airfoil geometries, covering a wide range of values of the
Reynolds number and the angle of attack, including the stall and post-stall regions.
The numerical model employs the Panel Method with a linear piecewise-
continuous vortex distribution along the panels to calculate the body contribution to the
flow. A cloud of discrete vortices with a Lamb core is used to model the vorticity field
within the lagrangian Discrete Vortex Method framework. The motion of these vortices
simulates the convective and diffusive transport of vorticity. The Random Walk Method
simulates the transport of vorticity by diffusion. The second-order Adams-Bashforth time-
marching scheme simulates the transport of vorticity by convection, where the convective
velocity is calculated from the Biot-Savart Law. The strength of the nascent vortices that
model the new vorticity generated on the body surface at each time step is calculated in
order to conserve total circulation and to satisfy the no-penetration boundary condition on
the body.
The results obtained for several angles of attack present good agreement with
experiments for NACA airfoils. These results show the enormous capacity of the
numerical algorithm and the physical models used for the simulation of complex high
Reynolds number flows.
v
AGRADECIMENTOS
Depois de cursar as disciplinas do curso de mestrado em Engenharia mecânica,
segue-se uma difícil etapa de elaboração de uma tese em determinada área de interesse.
Várias pessoas contribuíram das mais diversas formas para que este trabalho pudesse ser
concluído. Para algumas delas registro aqui meus agradecimentos.
Agradeço primeiramente a minha mãe, irmã e avó, as quais me deram apoio
incondicional para realização do curso de mestrado.
Agradeço especialmente a minha namorada Mayna, pela sua compreensão pelos
finais de semana de trabalho árduo e apoio nas horas difíceis.
Gostaria também de agradecer aos meus amigos, pela paciência em ouvir todos os
dias algo sobre Mecânica dos Fluidos e dinâmica de vorticidade, por mais inacreditável
que parecesse a existência dos tais vórtices.
Agradeço ao incentivo e aconselhamentos de meu orientador, Prof. Gustavo, que
sempre manteve um nível de exigências adequado para possibilitar um trabalho de boa
qualidade, além de ser botafoguense, o que contribuiu para o bom andamento do nosso
trabalho.
Agradeço também aos professores do departamento de Engenharia mecânica da
COPPE/UFRJ, os quais exibiram conceitos e fundamentos de grande valia para minha
formação profissional durante o curso de mestrado.
Meu sincero agradecimento para os colegas de turma com os quais cursei a maioria
das disciplinas do mestrado, em especial Alexandre e Rodolfo, com os quais era possível
encarar os desafios sempre com bom humor.
Agradeço a equipe do Laboratório de Mecânica dos Fluidos e Aerodinâmica da
COPPE/UFRJ, sobretudo Angelo e Vanessa, que desenvolveram excelentes trabalhos, os
quais consultei frequentemente e deram suporte a todas minhas dúvidas durante a
elaboração do programa desenvolvido.
vi
Agradeço a CAPES pelo apoio financeiro durante o primeiro ano de mestrado e
agradeço a própria COPPE pela oportunidade de realizar um curso de pos-graduação de
qualidade indiscutível em uma entidade pública.
Finalmente, agradeço à concessionária de serviço público de energia elétrica Furnas
Centrais Elétricas S.A. pelo apoio durante parte da execução desta tese de Mestrado
através de um projeto de P&D ANEEL, aprovado no Ciclo 2002/2003. Especificamente,
agradeço aos Engenheiros de Furnas Daniel Mesquita, Edson Teixeira e Maria do Carmo
Reis Cavalcanti.
vii
ÍNDICE
Capítulo 1 – Introdução ................................................................................................. 1
1.1 Motivação ............................................................................................................... 2
1.2 O Escoamento ao Redor de Geometrias Aerodinâmicas ....................................... 5
1.3 Objetivos ................................................................................................................ 9
Capítulo 2 – Revisão Bibliográfica ............................................................................. 10
2.1 Trabalhos Experimentais ...................................................................................... 10
2.2 Modelos Analíticos............................................................................................... 11
2.3 Método dos Painéis .............................................................................................. 11
2.4 Método de Vórtices .............................................................................................. 13
Capítulo 3 – Modelo Matemático ............................................................................... 19
3.1 Hipóteses Simplificadoras ....................................................................................19
3.2 As Equações de Movimento ................................................................................ 20
3.3 A Equação de Transporte de Vorticidade ............................................................ 22
3.4 A Equação de Poisson para Pressão ..................................................................... 23
3.5 Relação Velocidade–Vorticidade: A Lei de Biot-Savart...................................... 23
3.6 O Vórtice de Lamb ............................................................................................... 26
3.7 Cargas Aerodinâmicas na Forma Adimensional .................................................. 27
Capítulo 4 – Método dos Painéis ................................................................................ 29
4.1 Escoamento Potencial........................................................................................... 29
4.2 Método Numérico................................................................................................. 31
4.3 A Condição de Kutta ............................................................................................ 33
4.4 Escolha de Singularidades .................................................................................... 34
4.5 Transformação do Problema em um Sistema de Equações Lineares ................... 36
4.6 Cálculo do Campo de velocidade e das Cargas Aerodinâmicas........................... 40
viii
Capítulo 5 – Método de Vórtices Discretos ................................................................ 42
5.1 O Algoritmo Básico.............................................................................................. 42
5.2 Geração dos Vórtices............................................................................................ 43
5.3 Convecção dos Vórtices ....................................................................................... 44
5.4 Difusão dos Vórtices ............................................................................................ 46
5.5 Avanço Temporal ................................................................................................. 47
5.6 Atualização das Condições de Contorno.............................................................. 48
5.7 Tratamento dos Vórtices que Penetram no Interior do corpo............................... 49
5.8 Cálculo das Cargas Aerodinâmicas .......................................................................50
Capítulo 6 – Implementação Numérica ...................................................................... 52
6.1 Solução do Sistema de Equações Lineares........................................................... 52
6.2 Geração de Números Randômicos ....................................................................... 53
6.3 Raio do Núcleo do Vórtice de Lamb.................................................................... 53
6.4 Velocidades Auto-Induzidas na Nuvem de Vórtices............................................ 54
6.5 Campo de velocidade nas Proximidades do Corpo .............................................. 56
6.5.1 Utilização de Sub-painéis .............................................................................. 57
6.5.2 Método da Imagem........................................................................................ 57
6.6 Reflexão dos Vórtices........................................................................................... 58
6.7 Parâmetros Numéricos.......................................................................................... 59
6.8 Particularidades do Algoritmo.............................................................................. 61
6.9 Rotinas Numéricas................................................................................................ 62
6.10 Fluxograma......................................................................................................... 65
Capítulo 7 – Resultados e Discussões.......................................................................... 67
7.1 Teste Preliminar: Cilindro Circular ...................................................................... 68
7.2 Determinação dos Parâmetros Numéricos ........................................................... 68
7.3 Aerofólios Estudados ........................................................................................... 71
7.4 Camada Limite sobre um Aerofólio .................................................................... 73
7.5 Aerofólio NACA 0012 ........................................................................................ 77
ix
7.6 Aerofólio NACA 4415 ........................................................................................ 88
7.7 Outros Aerofólios ................................................................................................ 90
Capítulo 8 – Conclusões .............................................................................................. 91
Referências .................................................................................................................... 95
Apêndice A – Singularidades de Ordem Superior no Método dos Painéis .......... 105
Apêndice B – Geometria de Aerofólios NACA de 4 Dígitos .................................. 107
Apêndice C – Resultados para o Aerofólio NACA 0012 com Re = 1,7 x 105 ........ 109
Apêndice D – Resultados para o Aerofólio NACA 0012 com Re = 1,7 x 105 e
Ângulo de Ataque 90o ........................................................................ 142
Apêndice E – Resultados para o Aerofólio NACA 0012 com Re = 1,7 x 106 ........ 145
Apêndice F – Resultado para o Aerofólio NACA 0012 com Re = 1,7 x 105 e
Tempo de Simulação mais Longo .................................................... 150
Apêndice G – Resultados para o Aerofólio NACA 4415 com Re = 1,7 x 105 ....... 153
x
LISTA DE SÍMBOLOS
A - Matriz de influência
Ai - Coeficiente de influência
aij - Elemento da matriz de influência
b - Vetor de constantes
bi - Elemento do vetor de constantes
Bi - Coeficiente de influência
bt - Vetor de constantes atualizado
c - Corda do aerofólio
C - Curva Material
Cd - Coeficiente de arrasto
Cl - Coeficiente de sustentação
Cp - Coeficiente de pressão
d - Distância do vórtice ao corpo
dCd - Cd infinitesimal
dCl - Cl infinitesimal
dl - Elemento infinitesimal linear
dS - Elemento infinitesimal de área
dV - Elemento infinitesimal volumétrico
f - Função do contorno do aerofólio
k - Constante
L – Matriz triangular inferior
M - Número de vórtices discretos
n - Vetor normal à superfície do corpo
N - Número de painéis
Nsub - Número de sub-painéis
NT - Número de passos no tempo
p - Pressão absoluta
P - Número randômico
po - Pressão de referência
ps - Pressão de estagnação
q - Módulo de velocidade
xi
q - Campo vetorial de velocidade
Q - Número randômico
∞Q - Velocidade do escoamento uniforme
qt - Velocidade tangencial
r - Coordenada do sistema polar
r1 , r2 - Variáveis auxiliares
Re - Número de Reynolds
S - Superfície
s - Coordenada tangencial
t - Tempo
t - Vetor tangente
u - Velocidade na direção
U – Matriz triangular inferior
Vol - Volume
v - Velocidade na direção y
X - Coordenada horizontal global de um sistema cartesiano
x - Coordenada horizontal local de um sistema cartesiano
Y - Coordenada vertical global de um sistema cartesiano
y - Coordenada vertical local de um sistema cartesiano
α - ângulo de ataque
β - ângulo de inclinação dos painéis
δ - Distância normal do vórtice ao painel mais próximo
1δ - Espoessura de deslocamento da camada limite
i∆ - Comprimento do painel i
r∆ - Avanço difusivo radial
ip∆ - Variação de pressão sobre o painel i
x∆ - Deslocamento na direção x
y∆ - Deslocamento na direção y
θ∆ - Avanço difusivo tangencial
ε - Distância de geração de vórtices
ijkε - Tensor alternante
xii
φ - Potencial de velocidade
∞φ - Potencial de velocidade do escoamento incidente
Γ - Circulação
γ - Intensidade de vorticidade
λ - Intensidade de fonte
µ - Intensidade de um dipólo
ν - Viscosidade cinemática
π - Constante = 3.14156295...
θ - Coordenada do sistema polar
ρ - Massa específica
σ - Raio do núcleo dos vórtices
ω - Vorticidade
ω - Vetor vorticidade
ξ , ψ - ângulos auxiliares
Subscritos
i, j, k, p, n – índices
a – Referente ao início do painel
b – Referente ao início do painel
∞ - Relativo ao escoamento incidente
Superescritos
* - Variáveis dimensionais
a – Referente ao início do painel
b – Referente ao início do painel
^ - Para velocidades calculadas com vorticidade (circulação) unitária
¯ - Média
Operadores
∇ - Gradiente
xiii
⋅∇ - Divergente
×∇ - Rotacional
DtD - Derivada material
t∂∂ - Derivada parcial
dtd - Derivada ordinária
- Módulo
1
CAPÍTULO 1: INTRODUÇÃO
A Mecânica dos Fluidos é uma ciência extremamente rica em detalhes e,
conseqüentemente, de uma complexidade impressionante. Desde as civilizações mais
antigas busca-se o entendimento do escoamento de fluidos, desde a canalização de água
para o consumo, passando pela utilização dos ventos como força motriz de embarcações e
moinhos, até os complexos escoamentos de fluidos multifásicos na indústria de petróleo
atual.
Esta complexidade ainda fascina cientistas e pesquisadores em busca da completa
capacidade de modelar os fenômenos físicos envolvidos no escoamento de um fluido.
Durante o século XVII, Leonardo da Vinci elaborou desenhos de sistemas hidráulicos, um
deles podendo ser visto na Figura 1.1, mostrando já estar ciente da complexidade
envolvida no movimento de um fluido na simples descarga de uma tubulação.
Devido à crescente importância tecnológica da Mecânica dos Fluidos para
Engenharia, foram desenvolvidas gradativamente ferramentas matemáticas capazes de
representar a física desses fenômenos. As famosas equações de Navier-Stokes possuem
solução analítica fechada apenas para uma pequena classe de problemas, onde diversas
simplificações são permitidas. Mesmo as simplificações de Prandtl, que dão origem às
equações de camada limite, e a utilização da hipótese de efeitos viscosos desprezíveis, que
gera as equações de Euler, ainda não possuem solução analítica geral. Neste contexto é que
Figura 1.1 – Desenho esquemático elaborado por Da Vinci (SILVEIRA NETO, 2004).
2
surgem as soluções numéricas, auxiliadas enormemente pelo concomitante progresso da
computação eletrônica. A transformação de equações diferenciais parciais em sistemas de
equações algébricas e as conseqüentes soluções aproximadas obtidas, possibilitam a
análise e a melhor compreensão dos escoamentos complexos que ocorrem em problemas
reais de Engenharia.
Dentre os métodos numéricos desenvolvidos, um dos mais difundidos, possuindo
inclusive grande aplicação em programas computacionais comerciais, é o Método dos
Volumes Finitos, apresentado em PATANKAR (1980) e em MALISKA (2004). Este
método, juntamente com Método dos Elementos Finitos, pertence a uma classe de métodos
numéricos conhecidos como métodos Eulerianos, que possuem como característica
principal a discretização da região fluida com auxílio de uma malha computacional, onde
as grandezas de interesse são avaliadas.
O presente trabalho utiliza o Método de Vórtices Discretos, descrito
detalhadamente em LEWIS (1991). O Método de Vórtices Discretos é um caso particular
do Método de Partículas (COTTET & KOUMOUTSAKOS, 2000), o qual se caracteriza
por ser um método numérico que emprega uma abordagem Lagrangiana, isto é, que
acompanha o movimento de cada partícula individualmente. Esta metodologia não é
restrita à análise mecânica dos fluidos, porém adequa-se naturalmente à simulação de
escoamentos externos e incompressíveis que ocorrem comumente em projetos de
Engenharia, tais como: na análise aerodinâmica de veículos terrestres e aeronaves, nos
cálculos de esforços aerodinâmicos sobre edificações, na determinação de esforços
hidrodinâmicos sobre estruturas aquáticas ou submarinas, nos escoamentos dos diversos
tipos de turbomáquinas e na análise de interações fluido-estrutura em geral.
1.1 Motivação
A utilização de perfis aerodinâmicos em veículos terrestres, aquáticos e aeronaves
impulsiona fortemente a pesquisa e desenvolvimento de modelos para o cálculo de
escoamentos ao redor de superfícies aerodinâmicas. A crescente velocidade dos veículos e
a maior preocupação com a economia de combustível torna o estudo da aerodinâmica cada
vez mais importante. As aplicações para superfícies de controle em aeronaves requerem
um conhecimento ainda mais detalhado das propriedades do escoamento nestas regiões.
3
O comportamento do escoamento ao redor de perfis aerodinâmicos também é
fundamental para o projeto e a análise de escoamentos em turbomáquinas. Cada roda de
pás fixas ou móveis de um estágio de turbina a vapor ou a gás pode ser analisada através de
uma seção transversal ao longo da circunferência da roda das pás. Esta seção, quando
planificada, pode ser analisada como uma fileira de aerofólios, conforme a Figura 1.2.
Turbinas Eólicas de Eixo Horizontal (HAWT - Horizontal Axis Wind Turbine) é um tipo
de turbomáquina que vem ganhando importância, devido à grande preocupação com custo
e eficiência de sistemas eólicos para geração de energia elétrica. O rotor destas turbinas é
fabricado, em geral, com três pás dispostas radialmente a partir de um eixo central
horizontal, como pode ser visto na Figura 1.3a. As pás possuem um perfil aerodinâmico
específico, como mostrado na Figura 1.3b e, por isso, funcionam como uma asa, onde
forças de sustentação se desenvolvem sobre suas superfícies.
Figura 1.2 – Roda de palhetas de turbomáquinas analisada como fileira de aerofólios.
4
O desenvolvimento de turbinas do tipo HAWT requer o conhecimento detalhado do
escoamento do ar em torno das pás do rotor. É através do estudo desse escoamento que se
torna possível melhorar o desempenho aerodinâmico do rotor, permitindo, em última
análise, aumentar a capacidade da turbina através do aumento de seu tamanho e,
principalmente, de sua eficiência. A melhoria do desempenho do rotor pode ser realizada
de várias formas, podendo-se destacar o desenvolvimento de novos perfis e a otimização
da configuração geométrica das pás do rotor.
A passagem do ar por uma turbina eólica produz um escoamento à jusante do rotor
denominado Esteira, que se caracteriza por possuir uma velocidade menor do que o
escoamento do vento incidente, além de ser tridimensional e turbulento. A Figura 1.4
ilustra a esteira atrás de um rotor de uma turbina do tipo HAWT. Conseqüentemente,
turbinas posicionadas na esteira de outras turbinas produzem menos energia, pois o vento
incidente possui menos energia disponível. Este fenômeno é denominado Sombreamento,
ou Efeito Sombra. Assim, pode-se ressaltar duas grandes motivações para se estudar o
escoamento do ar em torno das pás de uma turbina:
(a) o aumento do desempenho aerodinâmico da turbina através de projetos mais
eficientes das pás do rotor;
(b) a redução do efeito sombra para o posicionamento otimizado de turbinas em
fazendas eólicas.
a) turbina de três pás b) Modelo tridimensional de uma pá
(REPOWER, 2004) ; (SILVA et al., 2004);
Figura 1.3 –Turbina tipo HAWT.
5
Naturalmente, há muitas outras motivações para esse estudo, podendo-se
mencionar a geração de ruído causado pelo escoamento do ar em torno do rotor e o projeto
mecânico da turbina, que é função direta das cargas aerodinâmicas atuando sobre ela.
Desta forma, fica clara a importância de se estudar, de um ponto de vista
fundamental, a aerodinâmica de uma pá de uma turbina do tipo HAWT. Devido à
complexidade do escoamento em torno de uma pá de turbina, que é tridimensional e
transiente, o estudo da aerodinâmica do rotor eólico passa, em primeiro lugar, pelo estudo
do escoamento bidimensional em torno do perfil aerodinâmico (ou aerofólio) que
compõem uma seção (ou elemento) da pá. A determinação de perfis aerodinâmicos
adequados para o escoamento através do rotor eólico e o desenvolvimento de metodologias
que permitam a análise teórica desses perfis são assunto de pesquisa básica e de primordial
importância para o desenvolvimento de projetos aerodinâmicos mais eficientes da turbina
como um todo. Este conhecimento compreende o alicerce no qual estudos tridimensionais
mais elaborados necessitam se apoiar para melhorar o desempenho de sistemas eólicos.
1.2 O Escoamento ao Redor de Geometrias Aerodinâmicas
O movimento de um fluido ao redor da superfície de um corpo sólido, produz
regiões onde os efeitos viscosos e inerciais atuam de forma diferenciada. Quanto maior for
o número de Reynolds, menor é a região onde os efeitos viscosos são importantes e vice-
versa. Esta diferença de comportamento do escoamento, aliada a fenômenos de transição
laminar-turbulenta, gera padrões de escoamento bastante diferenciados para uma mesma
geometria.
Figura 1.4 – Esteira formada por uma turbina de duas pás (DAHL et al.,1999).
6
Na classe de escoamentos externos, as regiões onde os efeitos viscosos são
importantes são denominadas camada-limite e esteira, ilustradas esquematicamente nas
Figuras 1.5 e 1.6. A camada limite é uma região formada nas proximidades do corpo pelo
efeito de aderência da camada de fluido que está em contato com a superfície, conhecida
como condição de aderência. A esteira é uma região formada a jusante do corpo, resultado
da própria perturbação que o corpo apresenta ao escoamento (BATCHELOR, 1967).
A camada limite e a esteira possuem um campo de velocidade característico,
diferente do escoamento externo a essas regiões. Nas regiões onde os efeitos viscosos são
consideráveis, assinaladas nas Figuras 1.3 e 1.4, existem camadas cisalhantes de fluido
com diferentes perfis de velocidade. Para números de Reynolds altos, essas camadas de
fluido são instáveis e tendem a se enrolar por efeito da tensão cisalhante viscosa, resultante
do gradiente de velocidade presente no escoamento, conforme mostra a Figura 1.7. Este
enrolamento característico de camadas cisalhantes é conhecido como instabilidades de
Kelvin-Helmholtz (SILVEIRA NETO, 2004).
As instabilidades formadas, sobretudo em escoamentos ao redor de corpos
rombudos, dão origem a uma esteira com regiões de escoamento turbilhonar alternado e
contra rotativo, a esteira de Von Karman. Uma foto de uma visualização experimental
desta esteira é apresentada na Figura 1.8.
Escoamento incidente
Camada Limite
Esteira
Figura 1.5 – Escoamento para Reynolds alto (BATCHELOR, 1967).
7
Escoamento incidente
Camada Limite
Esteira
Figura 1.6 – Escoamento para Reynolds baixo(BATCHELOR, 1967).
Intabilidades de Kelvin-Helmholtz
uuniforme
uesteira /
camada limite
Figura 1.7 – Formação de instabilidades no escoamento (SILVEIRA NETO, 2004)
Figura 1.8 - Esteira de Von Karman (SILVEIRA NETO, 2004)
8
No caso de geometrias aerodinâmicas, mesmo para um número de Reynolds fixo, o
escoamento pode apresentar padrões absolutamente diferentes, dependendo fortemente da
orientação relativa entre o corpo e o escoamento incidente. Por exemplo, para um aerofólio
como o mostrado na visualização experimental da Figura 1.9a, onde o ângulo de ataque é
baixo ou nulo (ângulo entre a corda do aerofólio e a velocidade do escoamento incidente –
ver capítulo 2), o escoamento comporta-se como se não houvesse efeitos viscosos. Nesse
caso há a formação de uma camada limite e de uma esteira extremamente finas. Já para
esta mesma geometria, quando posicionada em um ângulo de ataque alto, como mostrado
na Figura 1.9b, ocorre a descolamento da camada limite, fenômeno denominado “estol”
(stall) quando ocorre sobre a superfície de um aerofólio (SCHLICHTING & GERSTEN,
2000). Este tipo de escoamento forma instabilidades de diversos tipos e escalas, sendo sua
análise tridimensional um problema complicado de ser resolvido mesmo nos dias atuais.
Figura 1.9 – Escoamento ao redor de um aerofólio: (a) com baixo ângulo de ataque;
(b) com alto ângulo de ataque. (DEWI, 1998)
(a)
(b)
9
1.3 Objetivos
Este trabalho se propõe a elaborar modelos matemáticos e numéricos para estudar o
escoamento ao redor de um perfil aerodinâmico, utilizando-se para este fim o Método de
Vórtices Discretos. Especificamente, objetiva-se realizar um estudo numérico do
escoamento bidimensional, transiente e incompressível ao redor de perfis aerodinâmicos.
Os principais objetivos deste estudo são listados a seguir:
• Testar a utilização de distribuições de singularidades tipo vórtice de distribuição
linear e explorar suas vantagens quanto à continuidade da vorticidade ao longo
dos painéis sobre o corpo;
• Estudar a aplicação do Método de Vórtices associado ao Método dos Painéis
para o escoamento ao redor de aerofólios com um número consideravelmente
maior de vórtices que o empregado por ALCÂNTARA PEREIRA (1999);
• Verificar a capacidade do Método de Vórtices de representar o fenômeno de
“estol” sobre a superfície de aerofólios;
• Estudar a influência da curvatura no escoamento para aerofólios tipo NACA,
comparando os aerofólios NACA 0012 e NACA 4415;
• Aplicar uma forma alternativa de tratamento da conservação de circulação ao
redor do aerofólio.
Estes estudos envolvem o desenvolvimento de novos algoritmos do Método de
Vórtices específicos para utilização em simulações numéricas de escoamentos em torno de
perfis aerodinâmicos, para diferentes orientações relativas do escoamento incidente.
10
CAPÍTULO 2: REVISÃO BIBLIOGRÁFICA
2.1 Trabalhos Experimentais
Para a perfeita compreensão dos fenômenos físicos em Mecânica dos Fluidos,
assim como em quase toda a Engenharia, é de suma importância que medições
experimentais sejam feitas em condições controladas, seja no laboratório ou no campo. As
simulações numéricas buscam retratar os fenômenos físicos tão próximos da realidade o
quanto possível. Desta forma, a comparação com resultados experimentais adequados é
fundamental para a validação de um código computacional.
Especificamente no caso de aerofólios, existe um grande número de estudos
experimentais realizados ao longo dos anos para a determinação das cargas aerodinâmicas
que atuam sobre aerofólios das mais variadas famílias. Em geral, os experimentos
objetivam medir a distribuição do coeficiente de pressão ao longo da corda do aerofólio, e
a dependência que os coeficientes de sustentação, arrasto e momento possuem do ângulo
de ataque, tendo o número de Reynolds como parâmetro. Essas medições são tomadas em
condições de regime permanente. A literatura possui inúmeras publicações disponíveis que
compilam esses dados, sendo as mais utilizadas as citadas abaixo.
ABBOTT & VON DOENHOFF (1959) reúnem informações experimentais sobre
diversos tipos de aerofólios da família NACA, mostrando medições para o coeficiente de
sustentação em função do ângulo de ataque para números de Reynolds da ordem de 106 e
distribuição de pressão para alguns perfis.
Em conjunto com dados experimentais para corpos rombudos, tubulações e
diversas outras aplicações, BLEVINS (1984) apresenta diversos resultados para aerofólios,
incluindo escoamentos com Números de Reynolds mais baixos, da ordem de 105.
HOERNER & BORST (1985) expõem dados detalhados para aerofólios juntamente
com diversas teorias relacionadas ao escoamento de ar e água ao redor de aerofólios, flaps,
superfícies de sustentação em geral e aviões.
11
No caso específico de aerofólios para pás de turbinas eólicas, o extenso trabalho de
BERTAGNOLIO et al. (2001) apresenta resultados experimentais em comparação com
resultados numéricos no programa XFOIL e de um outro programa, ELLIPSYS 2D, do
próprio laboratório Risoe, mostrando boa concordância com respeito aos coeficientes de
sustentação e arrasto.
2.2 Modelos Analíticos
Inicialmente, a modelagem de escoamentos no campo da aerodinâmica era
realizada basicamente através de modelos analíticos, grande parte deles utilizando
elementos da Teoria Potencial. Alguns livros abordam o assunto com extrema clareza a
apresentam aplicações da Teoria Potencial, como por exemplo THWAITES (1960),
KARAMCHETTI (1980) e ANDERSON (1991).
Fundamentos básicos de aerodinâmica, e um desenvolvimento didático da Teoria
de Aerofólios Finos, além de algumas técnicas numéricas e diversas aplicações
aeronáuticas são apresentadas em BERTIN & SMITH (1998).
A apresentação das teorias analíticas em Mecânica dos Fluidos é feita também em
MILNE-THOMSON (1955), no campo da hidrodinâmica e em MILNE-THOMSON
(1958), onde se buscam aplicações aerodinâmicas fazendo uso da Transformação
Conforme. Para aplicações do Método das Perturbações é recomendado o livro de VAN
DYKE (1975).
A teoria de camada limite é apresentada de forma extensa detalhada, na última
edição de SCHLICHTING & GERSTEN (2000). Para detalhes específicos sobre a camada
limite laminar, pode-se consultar ROSENHEAD (1963).
2.3 Método dos Painéis
A utilização das teorias analíticas, tal como a Transformada Conforme, torna-se
extremamente trabalhosa para corpos de geometria arbitrária, devido à alta complexidade
atingida pelas manipulações algébricas. Neste tipo de aplicação é que os métodos
12
numéricos apresentam as suas maiores vantagens, pois podem ser aplicados, salvos
restrições específicas, para qualquer geometria.
Existe uma classe de métodos numéricos com grande utilização em Mecânica dos
Fluidos que se baseia na discretização da superfície do contorno, chamado Método dos
Elementos de Contorno BREBBIA et al.(1984). Um desses métodos é denominado
Método dos Painéis, que se destaca pela sua grande aplicação na indústria aeronáutica para
a análise de escoamentos potenciais ao redor de corpos rombudos e corpos aerodinâmicos.
Vários autores consideram que os primeiros passos para o desenvolvimento do Método dos
Painéis foram dados por MARTENSEN (1959), que descreve uma forma de analisar perfis
aerodinâmicos em escoamentos potenciais para aplicações aeronáuticas.
Diversos trabalhos decorrentes foram publicados, destacando-se o artigo de HESS
& SMITH (1966) que apresenta uma extensão do método para outras aplicações. Uma
revisão histórica e características detalhadas do método podem ser vistas em HESS (1990).
Para apresentações didáticas de descrição dos fundamentos do Método dos Painéis e das
teorias potenciais clássicas aplicadas a escoamentos bidimensionais e incompressíveis ao
redor de aerofólios, pode-se destacar os livros de KATZ & PLOTKIN (2001) e MORAN
(1984).
O Método dos Painéis permite uma ampla variedade de escolha de singularidades e
de suas formas de distribuição sobre a superfície discretizada do corpo. Um estudo
comparativo dessas possibilidades pode ser visto em PEREIRA et al. (2004). Aplicações
de singularidades de maior ordem são analisadas em PEREIRA & BODSTEIN (2004).
SILVA et al. (2004) utilizam o Método dos Painéis para cálculo bidimensional de
perfis de pás de turbinas eólicas, calculando a sustentação para baixos ângulos de ataque e
dimensionando o rotor aerodinamicamente segundo o procedimento descrito em
MANWELL et al. (2002).
Neste trabalho o Método dos Painéis é utilizado para calcular a contribuição do
corpo no escoamento, com intuito de possibilitar a simulação de uma geometria arbitrária.
13
2.4 Método de Vórtices
A literatura com respeito a métodos lagrangianos é consideravelmente extensa,
sobretudo no contexto do Método de Vórtices. Porém existem alguns artigos de revisão
sobre o assunto que fornecem uma ampla visão, não só dos trabalhos mais recentes nesta
área, mas também de sua evolução histórica. Dentre estes artigos pode-se destacar a última
revisão de SARPKAYA (1994), pela sua extensa e detalhada compilação de esquemas
numéricos e aplicações. Outros artigos similares são aqueles elaborados por SAFFMAN &
BAKER (1979), LEONARD (1980) e (1985), McCROSKEY (1982), AREF (1983), AREF
& KAMBE (1988), SARPKAYA (1989) e LEWIS (1999). O recente trabalho de
KAMEMOTO (2004) apresenta um resumo das aplicações do Método de Vórtices em
problemas práticos de Engenharia, mostrando as potencialidades do método.
A primeira tentativa de simular um escoamento através da modelagem do campo de
vorticidade por meio de vórtices discretos foi feita por ROSENHEAD (1931), cujos
cálculos manuais utilizando vórtices potenciais já mostravam resultados que apresentavam
uma grande semelhança com o escoamento estudado. Porém a simulação tornava-se
rapidamente instável, gerando um movimento caótico dos vórtices após alguns passos de
tempo. Mesmo com o uso de computadores para discretizar a vorticidade em um número
maior de vórtices, os resultados prosseguiam divergindo, uma vez que este esquema gera
altas velocidades induzidas.
Posteriormente passou-se a utilizar vórtices com núcleos finitos, onde as
velocidades dentro do núcleo são suavizadas por alguma função especial, permitindo,
assim, uma simulação mais real dos efeitos viscosos do escoamento. Vários esquemas para
limitação do campo de velocidade dos vórtices em simulações bidimensionais têm sido
propostos. SPREITER & SACKS (1951) utilizam o vórtice de Rankine, por ser este uma
primeira aproximação do núcleo de um vórtice real. KUWAHARA & TAKAMI (1973)
usam o vórtice de Lamb de núcleo variável, com o argumento que este descreve a solução
exata para a difusão de vorticidade de um único vórtice em um fluido viscoso. CHORIN &
BERNARD (1973), utilizam o vórtice de Lamb de núcleo fixo.
Em CHORIN (1973), a difusão viscosa é simulada através do Método do Avanço
Randômico, baseado na concepção do movimento Browniano de Einstein. LEWIS (1981)
14
mostra a associação do Método de Vórtices com o Método dos Painéis, facilitando a
extensão do Método de Vórtices para geometrias mais complexas. No mesmo ano,
PORTHOUSE & LEWIS (1981) adicionam o avanço randômico no modelo com o Método
dos Painéis, formando a base para o livro de LEWIS (1991) que trata da geração de
vorticidade na superfície do corpo e seu subseqüente desprendimento na esteira, indicando
algoritmos para simulação e incluindo programas em anexo.
KAMEMOTO (1994) classifica os métodos de vórtices em três classes: a primeira
libera vórtices apenas nos pontos de separação que são supostos conhecidos; a segunda
utiliza uma aproximação de camada limite para calcular os pontos de separação sobre a
geometria considerada e os vórtices são gerados a partir da espessura de deslocamento da
camada limite, nos pontos de separação calculados; a terceira classe modela a vorticidade
presente na camada limite por vórtices, que são gerados em toda a camada limite da
superfície do corpo e são transportados por convecção e difusão para todo o escoamento,
fazendo com que a separação ocorra como um resultado natural da simulação. Esta última
classe é mais poderosa e torna-se cada vez mais viável com a utilização de computadores
mais potentes e esquemas de aceleração do algoritmo convectivo, como o Partícula-Caixa
de BARNES & HUT (1986) e o Caixa-Caixa de GREENGARD & ROKHLIN (1987) ou
CARRIER et al. (1988). Estes esquemas de aceleração da etapa convectiva do algoritmo
reduzem a ordem de grandeza do contador de operações de M2 para MlogM e M,
respectivamente, onde M é o número de vórtices. Ainda no artigo de KAMEMOTO
(1994), é apresentado um esquema baseado na aproximação de camada limite, que sugere
que os elementos de vorticidade possuam seu comprimento paralelo ao escoamento mais
relevante. Desta forma, é proposta a liberação de uma linha (ou folha) de vorticidade com
comprimento finito, que permanece como tal até uma certa distância do corpo e, a partir
desse ponto, é transformada em um vórtice pontual com núcleo. Resultados para um
cilindro retangular e uma esfera são apresentados.
Para a simulação de aerofólios com “espoiler” (spoiler), XU (1998 e 1999) parte do
princípio de que a separação ocorre na quina viva da extremidade do espoiler e utiliza um
esquema com liberação de vórtices apenas nos pontos de separação, bordo de fuga e
extremidades do aerofólio, obtendo resultados razoáveis em comparação a experimentos.
15
Apresentando uma revisão para o tratamento da interação de vórtices com corpos
bidimensionais, mais especificamente aerofólios, MOOK & DONG (1994), discutem
simulações com desprendimento de vórtices apenas no bordo de fuga do aerofólio e
empregam um esquema especial para evitar distribuição irregular de vorticidade.
Ao longo do tempo, novos métodos também surgiram para simular a difusão
viscosa. O Método do Crescimento do Núcleo vem sendo utilizado com sucesso em alguns
algoritmos (KAMEMOTO, 1994). Este esquema gerou alguma polêmica por ter sido
criticado por GREENGARD (1985), que argumentou que o método aproximaria a equação
errada e a solução não tenderia à solução da equação de Navier-Stokes. OGAMI &
AKAMATSU (1991) propuseram um novo modelo, o Método da Velocidade de Difusão,
aplicado a um cilindro circular por OGAMI & AYANO (1995) e CLARKE & TUTTY
(1993). Este método apresenta problemas nos pontos em que a vorticidade tende a zero. O
recente Método de Redistribuição de Vorticidade, de SHANKAR & VAN DOMMELEN
(1996), também propõe um novo tratamento para difusão. Esses últimos três métodos
juntamente com o Avanço Randômico são comparados em TAKEDA & TUTTY (1997),
que chegaram à conclusão de que, para uma faixa de Reynolds de 300 a 9500, o Método de
Redistribuição de Vorticidade fornece a implementação mais robusta.
KOUMOUTSAKOS & LEONARD (1995) utilizam um esquema de aceleração
computacional e uma nova técnica de simulação da difusão denominada Intercâmbio da
Intensidade da Partícula (Particle Strength Exchange – PSE). Os resultados obtidos são
bastante representativos da interação do fluido viscoso com um cilindro circular em
estágios iniciais do escoamento.
COTTET et al. (2000) propõe um modelo de variação espacial do raio do núcleo
dos vórtices, tentando desta forma representar a deformação dos vórtices no escoamento e
apresentando resultados para cilindros circulares.
Outros trabalhos exploram uma variedade de aplicações do Método de Vórtices a
problemas de Engenharia, tais como o “estol” dinâmico de aerofólios e o escoamento em
cascatas de turbinas e compressores (LEWIS, 1991). Aplicações a escoamentos
tridimensionais também têm sido realizadas, podendo-se citar como exemplo os artigos de
OJIMA & KAMEMOTO (1999), que simulam o escoamento ao redor de uma esfera com
16
número de Reynolds de 300 a 1000 utilizando o método do crescimento do núcleo para a
difusão, e de PLOUHAMS et al. (2002) para o mesmo escoamento e a mesma faixa de
Reynolds, utilizando, porém, o método PSE para a difusão.
AKBARI & PRICE (2003) utilizam o Método de Vórtices para a simulação do
“estol” dinâmico sobre o aerofólio NACA 0012. A etapa convectiva é resolvida através do
esquema de vórtice-caixa e a etapa difusiva través de um esquema ADI (Alternating
Direction Implicit) de diferenças finitas de segunda ordem. A geometria do aerofólio é
mapeada através da transformação de Joukowski. Também são apresentados resultados
para a análise do “estol” convencional (estático).
O escoamento ao redor de turbinas eólicas é explorado por OJIMA &
KAMEMOTO (2001), que simulam o escoamento com os perfis NACA 0012 e MEL 012
e capturam em sua simulações o formato helicoidal esperado para a esteira formada a partir
da ponta das pás. Os resultados são de boa qualidade quando comparados aos
experimentos, considerando-se que o escoamento é tridimensional e bastante complexo.
Para modelos simplificados, onde o processo de geração de vorticidade é
implementado apenas nos pontos de separação, a distribuição de pressão sobre o corpo
pode ser calculada através da equação de Bernoulli não-permanente (VEZZA &
GALBRAITH, 1985). MUSTTO (2004) propõe uma metodologia mais geral para o
cálculo da distribuição de pressão sobre o corpo inspirada em LEWIS (1991), que toma
como base a equação de Navier-Stokes. Uma forma alternativa para se calcular a
distribuição de pressão sobre o corpo é resolver numericamente a equação de Poisson para
a pressão em todo o domínio fluido, descrita em NAKANISHI et al. (1993), que fornece
consequentemente o campo de pressão em todo o escoamento. Esta metodologia apresenta
bons resultados para simulações bi e tridimensionais. As forças aerodinâmicas podem
ainda ser calculadas a partir da integração da distribuição de pressão sobre a superfície do
corpo ou através da fórmula de Blasius (MUSTTO, 2004). SHINTANI & AKAMATSU
(1994) também apresentam uma metodologia alternativa para cálculo dos carregamentos
que necessita somente dos campos de vorticidade e velocidade do escoamento a cada passo
de tempo.
17
No âmbito nacional, a tese de mestrado de ALCÂNTARA PEREIRA (1999), trata
do escoamento em um cilindro circular e um aerofólio NACA 0012, a análise limita-se a
apenas um ângulo de ataque (5o) e o número de vórtices é limitado. A tratamento da
conservação de circulação é similar ao sugerido por LEWIS (1991) e pode ser aprimorado,
conforme descrito nos capítulos subseqüentes. Devido a limitações computacionais, ainda
é utilizado o esquema de eliminação de vórtices, o que reduz a resolução do campo de
vorticidade.
Na tese de doutorado de ALCÂNTARA PEREIRA (2004) e publicação decorrente
(ALCÂNTARA PEREIRA et al., 2004), o escoamento ao redor de uma grade linear de
turbomáquina utilizando função de interferência de grade, é analisado com o Método dos
Painéis de LEWIS (1991) e um Método de Vórtices com modelagem de turbulência via
Função Estrutura de Velocidade de segunda ordem, mostrando resultados promissores. A
adaptação do modelo de turbulência da Função Estrutura de Velocidade para esquemas
lagrangianos é descrito por HIRATA (2000) e aplicado ao escoamento ao redor de um
cilindro circular em ALCÂNTARA PEREIRA et al. (2002). Este modelo de turbulência
foi posteriormente explorado detalhadamente na tese de doutorado de MUSTTO (2004)
para o escoamento ao redor de um cilindro circular. A influência do chamado efeito solo
no escoamento sobre um corpo de geometria arbitrária é analisado na tese de doutorado de
RICCI (2002), que utilizou o Método de Imagens juntamente com o Método de Vórtices.
Diversos trabalhos de desenvolvimento e aplicação do Método de Vórtices ao
escoamento ao redor de corpos rombudos e aerodinâmicos têm sido realizados pela equipe
do Laboratório de Mecânica dos Fluidos e Aerodinâmica (LABMFA) da COPPE/UFRJ. O
escoamento não permanente ao redor de um aerofólio, com ponto de separação fixo,
durante sua interação com um vórtice potencial é modelado na tese de FONSECA (1996).
Nas teses de mestrado de MUSTTO (1998) e MALTA (1998), o escoamento ao redor de
um cilindro circular bidimensional é estudado e algumas alterações na implementação do
Método de Vórtices são propostas. Na tese de mestrado de CARREIRO (2002), são
estudados os escoamentos ao redor de cilindros bidimensionais com seção transversal na
forma de uma elipse, um aerofólio e uma placa plana, resolvendo o escoamento via
Transformação Conforme associada ao Método de Vórtices. Neste caso a transformação
conforme fornece a solução potencial exata, porém limita a simulação a geometrias com
transformação conhecida, o que é uma grande barreira quando trata-se de aerofólios.
18
Na tese de Doutorado de GUEDES (2003), são tratadas as geometrias de cilindro
retangular, cilindro circular e um interessante estudo de dois cilindros circulares com
diferentes posicionamentos relativos. Esta tese utiliza o Método dos Painéis com fontes e
vórtices de Lamb, satisfazendo às condições de impenetrabilidade e não-escorregamento
simultaneamente, e o Método de Vórtices para evolução da esteira.
A tese de doutorado de MUSTTO (2004), trata do escoamento ao redor de um
cilindro circular, utilizando o Teorema do Círculo e uma modelagem de turbulência via
Função Estrutura de Velocidade de Segunda Ordem. Novamente a utilização do Teorema
do Círculo limita a análise a geometrias mais simples, além de possuir um custo
computacional maior, por requerer o cálculo das velocidades induzidas pelas imagens dos
vórtices no interior do corpo.
Diante dos avanços já alcançados com o Método de Vórtices, particularmente
quando associado ao Método dos Painéis, esta revisão histórica indica a possibilidade de
um futuro promissor para esta técnica, motivando sua aplicação para a geometria de
aerofólios.
19
CAPÍTULO 3: MODELO MATEMÁTICO
Para a solução de qualquer problema real de Engenharia é necessário modelar os
fenômenos físicos presentes através da adoção de hipóteses simplificadoras e de equações
matemáticas que expressem corretamente a física envolvida. As hipóteses devem ser tais
que possibilitem a formulação de um problema matemático que seja bem posto e que
possua solução única. Deseja-se, em última análise, obter uma solução para o problema
formulado de tal maneira que os aspectos mais relevantes do fenômeno em estudo estejam
representados. Tendo em vista esse panorama, descreve-se a seguir o modelo utilizado e as
hipóteses simplificadoras.
3.1 Hipóteses Simplificadoras
Considere um aerofólio imerso em uma região infinita totalmente preenchida por
um fluido que inicia seu escoamento uniforme impulsivamente, isto é, o fluido assume
velocidade constante em todo o domínio instantaneamente. O escoamento incide sobre o
aerofólio fazendo um ângulo com relação à corda do aerofólio, chamado ângulo de ataque.
A Figura 3.1 mostra um desenho esquemático da situação descrita. Nesta Figura, ∞Q é a
velocidade do escoamento incidente, de módulo Q∞ , e α denota o ângulo de ataque.
Também estão representados os sistemas de coordenadas cartesiano (X,Y) e polar (r,θ ),
ambos com origem no bordo de ataque do aerofólio.
θ
α
∞Q X
Yr
Figura 3.1 – Diagrama esquemático do problema abordado.
Escoamento incidente
Corpo
Fluido
Fluido
20
Para o tratamento matemático do problema, são adotadas ainda as seguintes
hipóteses simplificadoras:
a dimensão do corpo no plano normal ao do escoamento é infinita, ou seja, é
muito maior que sua maior dimensão no plano do escoamento, de forma que o corpo pode
ser considerado como bidimensional;
o escoamento se inicia do repouso de forma impulsiva e objetiva-se limitar a
análise a tempos longos mas não extremamente longos, de forma a que efeitos
tridimensionais não se desenvolvam e o escoamento possa ser considerado bidimensional;
a região fluida é infinita e uniforme, ou seja, não há influência de barreiras ou
outros fenômenos interfaciais;
o aerofólio é o único corpo imerso no fluido;
o fluido é homogêneo, newtoniano, com propriedades constantes e seu
escoamento ocorre de maneira incompressível, considerando que as velocidades
envolvidas são muito menores que a velocidade do som no meio;
a força de campo gravitacional é desprezada, por possuir pouca influência no
escoamento.
Nota-se que não é feita a hipótese de regime permanente, pois a dinâmica do
escoamento provoca oscilações consideráveis do campo de velocidade com o tempo, um
dos fenômenos em estudo no trabalho.
3.2 As Equações de Movimento
As equações que regem o movimento de um escoamento incompressível de um
fluido newtoniano e com propriedades constantes, na forma dimensional, desprezando-se
as forças de corpo, podem ser escritas como:
0* =⋅∇ q [Equação da Continuidade] (3.1)
*2******
*
*
*
* 1 qqqqq∇+∇−=∇⋅+
∂∂
≡ νρ
ptDt
D [Equação de Navier-Stokes] (3.2)
21
Nestas equações, q* é o campo vetorial de velocidade, t* o tempo, p* a pressão estática
absoluta,ν é a viscosidade cinemática do fluido e ρ sua massa específica. O símbolo (*)
denota as varáveis na sua forma dimensional.
As equações (3.1) e (3.2) podem ser adimensionalisadas utilizando-se grandezas
características do escoamento:
• a velocidade característica, Q∞;
• o comprimento característico, escolhido como a corda do aerofólio c, indicada
na Figura 3.2;
• o tempo característico, determinado em função das variáveis anteriores: c / Q∞ .
Desta forma obtém-se o conjunto de variáveis adimensionais que possibilita a
análise de resultados para uma classe maior de problemas:
c
*∇≡∇ ;
cQtt ∞≡
*
; ∞
≡Qqq
*
; 2
*
∞
≡Qpp
ρ.
Assim, as equações da continuidade e de Navier-Stokes na forma adimensional
podem ser reescritas como
0=⋅∇ q , (3.3)
qqqqq 2
Re1
∇+−∇=∇⋅+∂∂
≡ ptDt
D , (3.4)
onde νcQ∞=Re é o número de Reynolds.
c
Figura 3.2 – Corda do Aerofólio
22
O sistema de equações diferenciais parciais (3.3) e (3.4) está sujeito às seguintes
condições de contorno:
q = 0 na superfície do corpo, que reflete as condições de impenetrabilidade e
não escorregamento para um corpo estacionário;
1=q , no infinito.
A condição inicial é de 1=q , em t = 0+, para toda região fluida.
3.3 A Equação de Transporte de Vorticidade
O rotacional do vetor velocidade q é chamado vetor vorticidade qω ×∇≡ . O vetor
vorticidade torna-se útil no estudo em questão devido à possibilidade de interpretar o
movimento de um fluido em um escoamento incompressível através da criação e
subseqüente evolução do campo de vorticidade no tempo. Em fluidos homogêneos, a
vorticidade é produzida somente nos contornos de superfícies. A vorticidade também pode
ser gerada no interior de fluidos não homogêneos (SARPKAYA, 1994).
Tomando-se o rotacional da equação (3.4) e utilizando a equação (3.3), obtém-se a
Equação de Transporte de Vorticidade (KUNDU & COHEN, 2002). Pode-se observar que
esta operação elimina o termo de pressão da equação original, fornecendo
ωqωωqωω 2
Re1
∇+∇⋅=∇⋅+∂∂
≡tDt
D . (3.5)
A interpretação para cada termo desta equação pode ser descrita como
(BATCHELOR, 1967):
t∂
∂ω representa a taxa de variação local de vorticidade;
ωq ∇⋅ representa o transporte convectivo (ou advectivo) de vorticidade;
23
ω2Re1
∇ representa o transporte difusivo de vorticidade devido à viscosidade;
qω ∇⋅ representa o esticamento e a deformação das tubos de vorticidade.
No caso bidimensional, o vetor vorticidade torna-se um escalar, possuindo apenas a
componente ω normal ao plano do escoamento não nula. Logo, a equação (3.5) se torna
uma equação escalar, dada por
ωωtω 2
Re1
∇=∇⋅+∂∂ q . (3.6)
É importante observar que o termo qω ∇⋅ é identicamente nulo no escoamento
bidimensional. Este fato representa uma limitação desta equação para a modelagem de
escoamentos turbulentos, pois o fenômeno da turbulência em fluidos está intimamente
ligado à deformação das linhas de vorticidade, como discutido em SARPKAYA (1994) e
SILVEIRA NETO (2004).
3.4 A Equação de Poisson para Pressão
Ao se tomar o divergente da equação de Navier-Stokes, ou seja, da equação (3.4),
recupera-se novamente o termo de pressão, expresso apenas em termos da velocidade. Esta
equação, denominada Equação de Poisson para a Pressão, pode ser expressa como
( )qq ∇⋅⋅−∇=∇ p2 . (3.7)
A equação (3.7) pode ser utilizada para o cálculo do campo de pressão em todo o
escoamento a partir do campo de velocidade, como por exemplo, em SHINTANI &
AKAMATSU (1994). Este tipo de metodologia possui a vantagem de calcular o campo de
pressão para todo o domínio com base nos campos de velocidade e vorticidade.
3.5 Relação Velocidade–Vorticidade: A Lei de Biot-Savart
Como o Método de Vórtices se baseia na solução da Equação de Transporte de
Vorticidade, equações (3.5 e 3.6), o campo de velocidade deve ser expresso em função do
24
campo de vorticidade. O campo de velocidade incompressível induzido pela vorticidade
concentrada em uma região finita é dado pela seguinte integral de volume (BATCHELOR,
1967),
( ) ( ) ( ) ( )∫∫∫ ′′−
′×′−−=
Vol
dVtt rrr
rωrrrq 3
,41,π
. (3.8)
A equação (3.8) é a Lei de Biot-Savart, desenvolvida inicialmente na teoria de
campos eletromagnéticos. Nesta equação, Vol representa o volume da região rotacional,
isto é, da região onde a vorticidade está concentrada. A integração é feita sobre este
volume através da variável auxiliar r’.
Um filamento de vorticidade é definido como um tubo de vorticidade, onde toda a
vorticidade está concentrada em um núcleo finito. Quando a área da seção transversal do
núcleo tende a zero e a vorticidade tende a infinito, o filamento é denominado filamento
singular de vorticidade, conforme ilustram as Figuras 3.3a e b. No caso de um filamento
singular de vorticidade de circulação constante Γ , a equação (3.8) torna-se
( ) ( ) ( )∫ ′−
′×′−Γ−=
Linha
34 rrrdlrrrq
π, (3.9)
onde dl é o vetor tangente à curva que representa o filamento de vorticidade. Caso o
filamento singular de vorticidade seja uma linha reta, este é chamado vórtice-linha e, no
caso bidimensional, de vórtice pontual, singular ou potencial (vide Figuras 3.3c e 3.3d). A
velocidade na direção θ induzida por um vórtice potencial, de circulação horária, é
( )r
rqπ2Γ
−= (3.10)
onde r=r . Equivalentemente, para um vórtice localizado em um ponto P(xo , yo), a
equação (3.10) pode ser escrita em coordenadas cartesianas como
25
( )( ) ( )222
oo
o
yyxx
yyu
−+−
−Γ=
π, na direção x, (3.11a)
( )( ) ( )222
oo
o
yyxx
xxv
−+−
−Γ−=
π, na direção y. (3.11 b)
Nas equações (3.10) e (3.11a, b) pode-se observar que existe uma singularidade para a
velocidade induzida calculada em um ponto no centro do próprio vórtice. Porém, sabe-se
que estruturas de vórtices reais não possuem esta singularidade. Um exemplo de uma
estrutura real que se assemelha à estrutura de um vórtice é um furacão, mostrado na Figura
3.4. Naturalmente, a velocidade no núcleo do furacão é finita, sugerindo a necessidade de
se utilizar um modelo que represente melhor as velocidades próximas ao centro do vórtice,
como apresentado na próxima seção.
(a) (b)
dl
Figura 3.3 – Filamento de Vorticidade: (a) filamento de vorticidade com núcleo de área finita, (b) filamento singular de vorticidade, (c) vórtice linha, (d) vórtice potencial e suas linhas de corrente.
r r’
(c) (d)
x z
y
26
3.6 O Vórtice de Lamb
Os dois modelos mais utilizados para a representação de vórtices reais são os
modelos de Rankine e Lamb (SARPAKAYA, 1994). Ambos os modelos suavizam as
velocidades induzidas a partir de uma distância ao centro do vórtice. A região circular nas
proximidades do vórtice que suaviza as velocidades induzidas pela região potencial do
vórtice é denominada núcleo. Um vórtice então é definido como a região cuja vorticidade é
não nula e que está circundada por uma região com vorticidade nula.
O modelo de Rankine sugere que o vórtice gire internamente como um corpo rígido
na região do seu núcleo e assuma o comportamento do vórtice potencial fora deste. Este
modelo introduz uma descontinuidade na derivada do campo de velocidade induzida
exatamente na circunferência do núcleo. O modelo de Lamb, utilizado por exemplo em
MUSTTO (1998), MALTA (1998), ALCÂNTARA PEREIRA (1999), GUEDES (2003) e
MUSTTO (2004), possui uma distribuição Gaussiana de vorticidade, que é uma solução
exata da equação de Navier-Stokes para um vórtice viscoso em um domínio infinito. Os
campos de vorticidade e de velocidade induzida para um vórtice de Lamb são dados por
( )
−
Γ−= 2
2
2 exp,σπσ
ω rtr , (3.12)
( )
−−
Γ−= 2
2
exp12
,σπr
rtrq , na direção θ , (3.13)
Figura 3.4 – Furacão se aproximando da costa: semelhança com asestruturas de vórtices (SILVEIRA NETO, 2004)
27
onde ( ) ( )22oo yyxxr −+−= , Re4t=σ é o raio do núcleo, (x , y) é a posição de um
ponto qualquer no escoamento, (xo , yo) é a posição do vórtice e Γ sua intensidade.
Na Figura 3.5 pode-se observar que o comportamento do campo de velocidade de
um vórtice de Lamb é contínuo e finito, e com derivadas contínuas e finitas, para todo r. O
campo de velocidade escrito em coordenadas cartesianas é descrito pelas seguintes
equações
( )
−−
−Γ= 2
2
2 exp12 σπ
rr
yyu o , (3.14)
( )
−−
−Γ−= 2
2
2 exp12 σπ
rr
xxv o . (3.15)
3.7 Cargas Aerodinâmicas na Forma Adimensional
De modo a garantir a generalidade dos resultados para problemas similares, as
equações para os carregamentos que atuam sobre a superfície do aerofólio também devem
ser expressas na forma adimensional.
A pressão ao longo da superfície do corpo imerso no escoamento é calculada
através do coeficiente de pressão, definido como
θq
Figura 3.5 – Velocidade tangencial induzida por um vórtice singular e um vórtice de Lamb.
28
2
21
∞
∗
≡Q
pC p
ρ, (3.17)
Para o cálculo efetivo dos carregamentos, basta integrar a equação (3.17) ao longo
da superfície do corpo, projetando as forças nas direções desejadas, ou seja:
na direção do escoamento, tem-se o coeficiente de arrasto:
dSCdCCS
pCorpo
dd ψcos∫∫ == ; (3.18)
na direção perpendicular ao escoamento, tem-se o coeficiente de sustentação:
dSCdCCS
pCorpo
ll ∫∫ == ψsen . (3.19)
Nas equações (3.18) e (3.19), dS é o elemento de área do corpo e βαψ += é o
ângulo entre a linha normal ao contorno do corpo e a direção normal ao escoamento
incidente. Este ângulo é calculado pela soma da inclinação do elemento de superfície e da
direção do escoamento incidente em relação à horizontal, conforme mostra a Figura 3.6.
A parcela correspondente ao atrito viscoso, originária da integração da tensão
viscosa sobre a superfície do corpo, não é avaliada devido a sua pequena contribuição no
coeficiente de sustentação, carregamento de maior interesse na análise proposta neste
trabalho.
∞Q
dCl
dCd
CpdS
ψ
α
β
Figura 3.6 – Carregamentos sobre um elemento de superfície
29
CAPÍTULO 4: MÉTODO DOS PAINÉIS
Os avanços encontrados na literatura científica na área de escoamentos potenciais
são inúmeros, podendo seu cálculo ser realizado por diversas abordagens diferentes.
Mesmo que o objetivo seja a simulação de um escoamento viscoso e rotacional, a solução
potencial pode auxiliar bastante a análise, inclusive compondo o algoritmo principal da
solução do problema, como será apresentado no capítulo seguinte.
Dentre os diversos métodos utilizados para a análise de escoamento potenciais,
existem técnicas analíticas, tais como o Método de Perturbações (VAN DYKE, 1975) e a
Teoria de Variáveis Complexas (CHURCHILL & BROWN, 1984). Dentro da Teoria de
Variáveis Complexas, pode-se mencionar ainda as técnicas de Transformação Conforme,
utilizada em CARREIRO (2002), e o Teorema do Círculo, utilizado em MUSTTO (2004).
Uma das técnicas numéricas mais utilizadas para o cálculo do escoamento potencial
é o Método de Elementos de Contorno (BREBIA et al., 1984). Como uma variação deste
método, pode-se destacar o Método dos Painéis. O Método dos Painéis é um método
numérico que possui a vantagem sobre as ténicas analíticas de ser aplicável a geometrias
arbitrárias, fornecendo resultados comprovadamente representativos e empregado até hoje
em diversos programas computacionais para projetos aeronáuticos. A aplicação do Método
dos Painéis juntamente com o Método de Vórtices também é uma estratégia
freqüentemente adotada, como por exemplo em LEWIS (1991), GUEDES (2003) e
ALCÂNTARA PEREIRA (1999).
Este capítulo destina-se a descrever a metodologia empregada para utilização do
Método dos Painéis no trabalho desenvolvido.
4.1 Escoamento Potencial
O escoamento irrotacional é aquele em que o vetor vorticidade é nulo. Pode-se
mostrar matematicamente (BATCHELOR, 1967) que todo escoamento irrotacional possui
um campo de velocidade associado que pode ser expresso como o gradiente de uma função
30
escalar φ, denominada Potencial de Velocidade. Assim, se ω = 0, então q = ∇φ, e
se q = ∇φ, então ω = 0. Por esse motivo, o Escoamento Irrotacional também é denominado
Escoamento Potencial.
O escoamento potencial incompressível é governado pela Equação de Laplace para
o potencial de velocidade, ∇2φ = 0, que é obtida ao se substituir q = ∇φ na equação da
continuidade para escoamento incompressível, ∇⋅q = 0. Considerando, então, o
escoamento potencial, bidimensional e incompressível ao redor de um corpo rígido em
repouso, pode-se formular o seguinte problema de valor de contorno para o potencial de
velocidade φ , escrito na forma adimensional,
( ) 02 =∇ rφ , no fluido, (4.1)
( ) 0=∂∂
=⋅∇n
nr φφ , na superfície do corpo, (4.2a)
( ) 1,lim =∇∞→
trr
φ , no infinito. (4.2b)
O campo de velocidade pode ser obtido através da expressão q = ∇φ, ou, escrito em
coordenadas cartesianas,
Xu
∂∂
=φ e
Yv
∂∂
=φ . (4.3a, b)
O vetor n é o vetor unitário normal à superfície do corpo, apontando para dentro,
enquanto que o vetor r(X,Y) representa uma posição qualquer no sistema de coordenadas
fixo no corpo em um instante t, conforme mostrado na Figura 4.1.
Figura 4.1 – Descrição do Sistema de Coordenadas fixo no corpo.
X
Y
n
Corpo Esteira
Q∞= 1 α
Figura 4.1 – Descrição do Sistema de Coordenadas fixo no corpo.
31
Utilizando a identidade de Green, pode-se mostrar (KATZ e PLOTKIN, 1991) que
a solução para a equação (4.1) pode ser escrita como
( ) ( ) ∫∫ ∞
+
++∂∂
−=CorpoEsteiraCorpo
dSdS φπλ
πγφ rr
nr ln
2ln
2. (4.4)
Nesta equação é possível notar que o argumento da primeira integral corresponde
ao potencial de um vórtice (ou dipolo) de intensidade γ , assim como o argumento da
segunda corresponde ao potencial de uma fonte (ou sumidouro) de intensidade λ . Já o
último termo, φ∞, corresponde ao potencial do escoamento externo (escoamento uniforme,
no caso estudado). Vale lembrar que a integração do termo de fontes é realizada somente
na superfície do corpo. Portanto, o gradiente do potencial das fontes é contínuo e não
suporta carregamento.
4.2 Método Numérico
Devido às características matemáticas das funções que descrevem os potenciais das
singularidades empregadas, com velocidades induzidas que variam com r de acordo com a
função 1/r, a influência das singularidades na solução global do escoamento tende a zero
para pontos muito distantes, satisfazendo automaticamente a condição de contorno (4.2b).
Resta, portanto, aplicar a condição de contorno (4.2a) na equação (4.4).
A tarefa de satisfazer à condição de contorno expressa pela equação (4.2a)
analiticamente em toda a fronteira de um corpo com geometria arbitrária é complicada. A
utilização do Método dos Painéis torna essa tarefa bastante simples. No Método dos
Painéis (KATZ & PLOTKIN, 1991), a superfície do corpo é discretizada em pequenos
segmentos, retos ou curvos, denominados painéis, conforme mostra a Figura 4.2. Os
pontos inicial e final de cada painel são chamados nós, ou pontos nodais, e o ponto central
é denominado ponto de controle. Sobre os painéis são distribuídas singularidades, tais
como fontes, dipolos ou vórtices. Desta forma, as condições de contorno são impostas
apenas nos pontos de controle dos painéis.
32
No modelo utilizado neste trabalho, o aerofólio é discretizado em N painéis retos ao
longo da superfície do corpo, sobre os quais são distribuídos vórtices. As intensidades
desses vórtices podem apresentar distribuição constante ou linear ao longo do painel. Para
o caso de distribuição constante de vorticidade ao longo do painel, a condição de contorno
(4.2a) pode ser escrita como
nn ⋅−∇=⋅
+ ∞==
∑∑ φλγN
jjj
N
jjj BA
11
, na superfície do corpo. (4.5)
No caso de distribuições lineares, a intensidade das singularidades depende do seu
valor nos extremos do painel, e a equação (4.5) assume a forma geral
( ) ( ) nn ⋅−∇=⋅
′++′+ ∞=
+=
+ ∑∑ φλλγγN
jjjjj
N
jjjjj BBAA
11
11 . (4.6)
Na equação acima, Aj, A’j, Bj e B’j são denominados coeficientes de influência, que
dependem apenas da geometria do corpo, e YX ααφ sencos +=∞ é o potencial de
velocidade do escoamento uniforme que incide sobre o corpo.
As equações (4.5 e 4.6) são numericamente equivalentes à condição de contorno de
Neumann, equação (4.2a), e são usadas em cada ponto de controle para determinar a
intensidade das singularidades escolhidas. Este procedimento origina um sistema de N
Pontos de Controle
5 4
Painel
7
6
3 2
1
98
1011
Figura 4.2 - Aerofólio discretizado, mostrando os pontos de controle e a numeração dos pontos nodais.
33
equações algébricas lineares e N incógnitas, no caso da distribuição de vorticidade
constante, e um sistema de N equações e N+1 incógnitas, no caso da distribuição de
vorticidade linear. Neste caso, há um valor da intensidade da singularidade para cada
extremidade dos painéis, conforme representado na Figura 4.3. Uma vantagem
significativa da distribuição linear é a continuidade da função nos pontos nodais (exceto no
bordo de fuga – Figura 4.3). Esta modelagem representa melhor a distribuição real de
vorticidade ao longo da superfície do corpo. Singularidades de ordem superior também
podem ser utilizadas (PEREIRA & BODSTEIN, 2004), requerendo maior trabalho
analítico para o cálculo dos coeficientes de influência.
4.3 A Condição de Kutta
Para escoamentos permanentes ao redor de aerofólios, a condição de Kutta deve ser
imposta para que haja igualdade das velocidades na superfície superior e inferior do bordo
de fuga. Alternativamente, a condição de Kutta pode ser considerada como uma condição
de continuidade de pressão no bordo de fuga. Expressa matematicamente através da
vorticidade no bordo de fuga, a condição de Kutta é escrita, em geral (KATZ & PLOTKIN,
1991), como
01 =+ Nγγ , (4.7a)
ou
011 =+ +Nγγ , (4.7b)
1+Nγ
1−Nγ Nγ1−Nγ
Nγ
a) constante b) linear
Figura 4.3 – Distribuições de singularidade.
34
sendo a condição (4.7a) para o caso de distribuição de vorticidade constante sobre os
painéis, e a condição (4.7b) para o caso de distribuição de vorticidade linear.
A equação (4.7b) é adicionada ao sistema de equações (4.6), fechando o problema
para o caso de distribuição linear de vorticidade. No caso da distribuição de vorticidade
constante, a equação (4.7a) deve substituir uma das equações (4.5) ou ser adicionada ao
sistema (4.5), formando um novo sistema de equações algébricas lineares que é
sobredeterminado. Este sistema deve ser resolvido empregando-se, por exemplo, o Método
dos Mínimos Quadrados (KREYSZIG, 1999).
4.4 Escolha de Singularidades
A escolha das singularidades adequadas para a correta modelagem do problema em
estudo deve ser feita de acordo com cada situação física considerada. Embora a equação
(4.5) ou (4.6), seja respeitada, há ainda um certo grau de arbitrariedade nesta escolha, e
faz-se necessário observar se o conjunto de equações utilizado é capaz de representar o
fenômeno físico em questão. Para o caso de escoamento externo sobre aerofólios, a
utilização de fontes unicamente só pode ser aplicada para aerofólios simétricos com ângulo
de ataque nulo, ou seja, situações sem forças de sustentação presentes. Para situações em
que há o desenvolvimento de uma força de sustentação, é necessário empregar
singularidades que possam representar uma assimetria no escoamento, tais como vórtices
ou dipólos. Para as simulações de escoamentos ao redor de aerofólios utilizando o Método
de Vórtices acoplado ao Método dos Painéis a serem realizadas neste trabalho, são
empregadas apenas singularidades do tipo vórtice com distribuição linear de vorticidade.
Este tipo de singularidade apresenta resultados satisfatórios para os objetivos deste
trabalho, conforme apontado por PEREIRA et al. (2004).
A distribuição de vorticidade linear pode ser decomposta na soma de uma
distribuição constante e de uma distribuição linear pura ( 0=γ em x = x1), como mostra a
Figura 4.4a. Assim, as equações para o potencial de velocidade e para as componentes de
velocidade nas direções X e Y induzidas pela distribuição de vorticidade linear em um
ponto qualquer no sistema de coordenadas do painel (Figura 4.4b) são dadas,
respectivamente, por
35
( ) ( )
( ) ( )
−−−
−−+−+
−−
+
+−−−=
1
221
2
2
222
2
121
2
12
1
22211
222ln
2
ln2
θθπ
γγ
θθπ
γφ
yxxyxxxxy
rr
xyxx
rryxxxx
ab
a
, (4.8)
( ) ( ) ( )
+−
−−
+−=1
212
1212 ln
22 rr
yxxx
u abap θθ
πγγ
θθπ
γ, (4.9a)
( ) ( ) ( )
+−−−
−−
+=1
21212
121
2 ln2
ln2 r
rxyxx
xxrr
v abap θθ
πγγ
πγ
, (4.9b)
onde
( ) ( )21
211 yyxxr −+−= e ( ) ( )2
22
22 yyxxr −+−= , (4.10a e b)
−−
= −
1
111 tan
xxyy
θ e
−−
= −
2
212 tan
xxyy
θ . (4.11a, b)
As equações (4.8) e (4.9a,b) estão expressas no sistema de coordenadas local do
painel. Portanto, é necessário transformar as coordenadas do sistema global (X,Y) para o
sistema de coordenadas local (x,y). O posicionamento dos dois sistemas pode ser
x
( )xγy
x x1 x2 x1 x2
2θ1θ
r2 r1
aγ
(x,y)
Figura 4.4 - Distribuição linear de vorticidade ao longo do painel. (a) Função matemática, (b) Variáveis auxiliares no sistema de coordenadas do painel.
(a) (b)
bγ
36
observado na Figura 4.5, onde é indicado o ângulo iβ para a inclinação do painel i, através
do qual é possível fazer a rotação dos eixos, utilizando a operação matricial
−−
−=
1
1
cossensencos
yYxX
yx
ii
ii
ββββ
. (4.12)
Através da Figura 4.5 também é possível determinar os vetores normal e tangente à
superfície do corpo no painel i, que são úteis para a especificação das condições de
contorno. Este painel possui vetor normal na direção de seu eixo y local e vetor tangente ao
longo de seu eixo x local. Logo,
=
i
ii β
βcossen
n e
−
=i
ii β
βsen
cost . (4.13)
As expressões para velocidades induzidas para outras singularidades e ordens
superiores estão apresentadas no Apêndice A.
4.5 Transformação do Problema em um Sistema de Equações Lineares
Em algumas regiões da superfície do aerofólio, tais como nas vizinhanças do bordo
de ataque e do bordo de fuga, o escoamento apresenta variações muito acentuadas dos
campos de potencial de velocidade, vetor velocidade e pressão. De modo a calcular com
melhor precisão essas variações, é conveniente implementar um refinamento da
discretização do corpo nessas regiões. Neste trabalho, uma discretização especial é usada
para aumentar o número de painéis nas regiões dos bordos de ataque e de fuga, como
X
Y
x
y
iβ
Figura 4.5. Sistema de coordenadas global e local.
37
mostrado na Figura 4.6. Assim, mesmo para um número moderado de painéis, a condição
de contorno de não penetração é satisfeita em mais pontos destas regiões críticas. A
discretização dos painéis sobre o corpo é, então, realizada utilizando-se as expressões
( )ξcos12
−=CX e ( )XfY = , (4.14)
onde o ângulo ξ está indicado na Figura 4.6, e f é a função que descreve o contorno do
aerofólio.
Para avaliar numericamente os coeficientes da equação (4.6) é necessário calcular
as velocidades induzidas por cada painel sobre cada ponto de controle. Essas velocidades
dependem da intensidade da vorticidade nas extremidades de cada painel quando a
distribuição de vorticidade é linear. Decompondo a velocidade induzida como a parte que
depende de jγ , denotada por ( )a, e a parte que depende de 1+jγ , denotada por ( )b, é
possível colocar em evidência os termos que dependem de aγ e bγ nas equações (4.9a,b) e
escrever
X
Y
Figura 4.6 – Discretização especial do contorno, concentrando painéis nas extremidades (KATZ & PLOTKIN, 1991).
ξ
C
38
( ) ( )( )
−−−
−=
1
2122
12
ln2 r
ryxx
xxu aa θθ
πγ
, (4.15a)
( ) ( )
+−
−=
1
212
12
ln2 r
ryx
xxu bb θθ
πγ
, (4.15b)
( ) ( ) ( )
−+−+−
−=
1
22122
12
ln2 r
rxxyx
xxv aa θθ
πγ
, (4.16a)
( ) ( ) ( )
+−−−
−=
1
21212
12
ln2 r
rxyxx
xxv bb θθ
πγ
. (4.16b)
Assim sendo, as componentes da velocidade total induzida por um painel j sobre o
ponto de controle de um painel i é
1ˆˆ
ˆˆ
+
+
=
j
ijb
b
jij
a
a
ijp
p
vu
vu
vu
γγ , (4.17)
onde a simbologia (^) é utilizada para designar a velocidade induzida por um painel com
vorticidade unitária, ou seja, com γa = γb = 1. Este procedimento facilita o cálculo dos
coeficientes de influência, como se mostra a seguir.
A velocidade induzida no primeiro ponto de controle, por exemplo, é
1111,1
31312
21211
1111
ˆˆ
ˆˆ
ˆˆ
ˆˆ
ˆˆ
ˆˆ
ˆˆ
ˆˆ
+
−
+
+
+
++
+
+
+
+
=
NN
b
b
NN
a
a
Nb
b
a
a
b
b
a
a
b
b
a
a
p
p
vu
vu
vu
vu
vu
vu
vu
vu
vu
γγ
γγγ
(4.18)
A equação (4.18) pode, portanto, ser generalizada na forma
12 1,
11
ˆˆ
ˆˆ
ˆˆ
ˆˆ
+= −
+
+
+
=
∑ N
iNb
bN
jj
ija
a
jib
b
ia
a
ip
p
vu
vu
vu
vu
vu
γγγ . (4.19)
39
Pela condição de contorno (4.6), esta velocidade induzida, quando adicionada à
velocidade induzida pelo escoamento incidente e projetada na direção normal ao painel,
deve ser nula para cada ponto de controle, gerando, assim, N equações tal que
iiip
p
vu
nn
−=
αα
sencos
, i = (1, 2, ... , N). (4.20)
Para montar o sistema de equações (4.20) na forma matricial, basta calcular as
velocidades induzidas com γj = 1. Portanto, para cada coeficiente da matriz, tem-se
( )iiijp
pi
ijp
pij v
uvu
a ββ cossenˆˆ
ˆˆ
=
= n , i = (1, 2, ... , N) , j = (1, 2, ... , N+1). (4.21)
Conseqüentemente, o vetor de constantes se torna
( )iiiib ββαα
αα
cossensencos
sencos
=
= n , i = (1, 2, ... , N). (4.22)
Utilizando ainda, a condição de Kutta, equação (4.7b), chega-se a um sistema de
equações algébricas lineares N+1 × N+1, onde as intensidades dos elementos de
vorticidade jγ são as incógnitas. Deste modo,
=
+
+
+
+
01001
2
1
1
2
1
1,,2,1,
1,2,22221
1,1,11211
N
N
NNNNNNN
NN
NN
b
bb
aaaa
aaaaaaaa
γγ
γγ
. (4.23)
Este sistema pode ser resolvido por técnicas tradicionais de solução de sistemas
lineares. A matriz A não necessita de tratamento especial para a sua inversão por ser
diagonal dominante.
40
4.6 Cálculo do Campo de Velocidade e das Cargas Aerodinâmicas
Uma vez determinadas as intensidades γj, é possível calcular o campo de velocidade
através da equação (4.19), sendo agora i um ponto qualquer do domínio.
Uma alternativa para se determinar a velocidade tangencial, qt , próxima ao painel é
considerar uma curva material C fechada próxima à superfície que envolve um painel de
comprimento j∆ (Figura 4.7), e calcular a circulação ao longo desta curva pela equação
( )∫ ∫==Γpainel Curva
qdsdxxC
γ , (4.24)
onde s é a variável ao longo da curva C. Como a velocidade normal muito próxima ao
painel é nula e a velocidade no interior do corpo também é nula, então, para uma
distribuição linear de vorticidade no painel, tem-se
( )2
1++= jj
jtqγγ
. (4.25)
Uma vez determinadas as velocidades tangenciais nos pontos de controle dos N
painéis, é possível utilizar a equação de Bernoulli para determinar o coeficiente de pressão,
ou seja,
superfície do corpo
Figura 4.7 – Circuito fechado próximo ao painel
jγ
1+jγ
j∆
Curva C
41
( )itip qC 21−= . (4.26)
Uma comparação entre os resultados obtidos com o Método de Painéis que utiliza
uma distribuição linear de vorticidade e uma solução analítica obtida via transformação
conforme para um aerofólio tipo Van de Vooren são mostrados na Figura 4.8. A Figura
apresenta a distribuição do coeficiente de pressão ao longo da corda do aerofólio. Com
base nesses resultados fica patente a capacidade do método de painéis de calcular
adequadamente a solução potencial para o escoamento bidimensional incompressível ao
redor de um aerofólio.
Figura 4.8 – Resultados numéricos para o coeficiente de pressão ao
longo de um aerofólio tipo Van de Vooren obtidos utilizando-se vórtices
lineares (com condição de Neumann, 120 painéis e α =10°)
comparados com a solução analítica obtida via transformação conforme
PEREIRA et al. (2004).
0 0.2 0.4 0.6 0.8 1-1
0
1
2
3
4
5
x/c
-Cp
Vórtice Linear (c.c. Neumann)Sol. Analítica Vande Vooren 15%
-Cp
x/c
42
CAPÍTULO 5: MÉTODO DE VÓRTICES DISCRETOS
Este capítulo detalha os modelos matemáticos e numéricos que compõem o Método
de Vórtices desenvolvido para simular o escoamento bidimensional de um fluido
newtoniano que escoa de forma incompressível e não permanente ao redor de um corpo
aerodinâmico.
5.1 O Algoritmo Básico
O Método de Vórtices Discretos consiste, fundamentalmente, numa ferramenta
numérica extremamente adequada para simular a dinâmica da vorticidade em escoamentos
externos. O método baseia-se na modelagem da vorticidade presente no escoamento como
uma nuvem de vórtices discretos, cujo campo de velocidade induzido pela nuvem é
resultado da dinâmica do escoamento como um todo. A evolução temporal do escoamento
se dá devido ao transporte convectivo e difusivo de vorticidade que é regido pela Equação
de Transporte de Vorticidade, apresentada no capítulo 3 e repetida aqui
ωqωωqω 2
Re1
∇+∇⋅=∇⋅+∂∂
t. (5.1)
A vorticidade gerada na superfície do corpo é calculada nesse trabalho a partir do
Método dos Painéis com vórtices lineares. A vorticidade de cada painel é transformada em
um vórtice discreto para formar a nuvem, que passa a se mover conforme o campo de
velocidade local do escoamento.
Como visto anteriormente, em um escoamento idealizado como bidimensional, a
equação (5.1) pode ser simplificada, pois o primeiro termo do lado direito da equação
(termo de esticamento e deformação dos tubos de vorticidade) é nulo e a vorticidade,
qω ×∇= , torna-se um escalar. Logo, em escoamentos bidimensionais, ω possui apenas
uma componente não nula, a componente normal ao plano do escoamento. Sob essas
simplificações, a equação (5.1) se torna
43
ωωtω 2
Re1
∇=∇⋅+∂∂ q . (5.2)
Embora os efeitos convectivos e difusivos ocorram simultaneamente, pode-se
supor, do ponto de vista numérico, que eles ocorram sucessivamente dentro de uma mesmo
intervalo de tempo t∆ . Com base nessa idéia, CHORIN (1973) propôs que o operador
convectivo-difusivo da equação (5.2) seja decomposto em dois operadores, um puramente
convectivo e outro puramente difusivo, governados respectivamente, pelas equações
0=∇⋅+∂∂ ω
tω q , [Convecção], (5.3)
ωtω 2
Re1
∇=∂∂ , [Difusão]. (5.4)
Sendo assim, a simulação é realizada de tal forma que a etapa convectiva é
calculada através da equação (5.3) ignorando-se os efeitos viscosos e, dentro do mesmo
intervalo de tempo, resolve-se a etapa difusiva através da equação (5.4) ignorando-se os
efeitos convectivos. Esta aproximação é tão melhor quanto menor for o incremento t∆ .
Um diagrama esquemático simplificado do algoritmo básico do método é
apresentado na Figura 5.1.
5.2 Geração dos Vórtices
Para simular os efeitos de geração de vorticidade na fronteira do fluido com o corpo
são gerados vórtices a partir das distribuições de vorticidade calculadas em cada painel.
Esses vórtices possuem intensidade igual à vorticidade média no painel, que, no caso de
Solução Potencial
Convecção de vórtices
Difusão de vórtices
Avanço no tempo
Geração de vórtices
Figura 5.1 – Algoritmo esquemático do Método de Vórtices
Geometria e cond. inicial
44
distribuição linear, resulta na média aritmética das intensidades nas extremidades. Logo,
para um vórtice k, sua intensidade Γk é dada por
kkk
Painelk dxx ∆
+==Γ +∫ 2
)( 1γγγ , k = 1, 2, ... , N. (5.5)
Os vórtices gerados são posicionados a uma distância ε normal ao ponto de
controle de cada painel, conforme a Figura 5.2 ilustra, sendo estes vórtices livres para se
mover por convecção e difusão, ao final da etapa de geração.
Considerando a camada limite laminar e bidimensional, a distância ε é escolhida
de tal forma que a geração de vorticidade coincida com a espessura de deslocamento da
camada limite, conforme sugerido em ROSENHEAD (1988).
5.3 Convecção dos Vórtices
A equação (5.3), que governa o escoamento incompressível, bidimensional e
desprezando os efeitos viscosos descreve a evolução de um escalar conservativo regido
apenas pela convecção no fluido. A conseqüência disto é que, no movimento de um fluido
de massa específica uniforme, escoando incompressivelmente, desprezando os efeitos
viscosos, a circulação ao redor de uma curva material fechada é constante (Teorema de
Kelvin) e as linhas de vorticidade movimentam-se com o fluido (Teorema de Helmholtz),
como discutido em BATCHELOR (1967). Portanto, o movimento de um vórtice potencial,
que ocupa a posição de vetor posição rk no instante t, é determinado pela taxa de variação
do vetor posição, ou seja,
εε
ε
superfície do corpo
painéis pontos de controle
vórtices gerados
Figura 5.2 – Geração de vorticidade na superfície
45
( )tt kk ,rq
r=
∂∂
. (5.6)
Esta expressão lagrangiana, adotada tanto para escoamentos viscosos como
desprezando os efeitos viscosos, é utilizada para calcular a nova posição de cada vórtice da
nuvem devido ao transporte convectivo de vorticidade, em cada passo de tempo.
Para o cálculo da velocidade em determinado ponto do domínio é necessário
adicionar os três campos de velocidade atuantes:
campo de velocidade devido ao escoamento uniforme:
=
∞
∞
αα
sencos
vu
; (5.7)
campo de velocidade devido à presença do corpo (painéis):
12 1,
11
ˆˆ
ˆˆ
ˆˆ
ˆˆ
+= −
+
+
+
=
∑ N
iNb
bN
jj
ija
a
jib
b
ia
a
ip
p
vu
vu
vu
vu
vu
γγγ ; (5.8)
campo de velocidade devido à nuvem de vórtices livres:
ik
M
kk
iv
v
vu
vu
∑=
Γ=
1 ˆˆ
. (5.9)
Na equação (5.9) as velocidades u e v são obtidas pelas equações (3.14) e (3.15),
fazendo 1=Γ , e M é o número de vórtices presentes na esteira.
Portanto a velocidade total convectiva sobre um vórtice é calculada por
+
+
=
∞
∞
vu
vu
vu
vu
v
v
p
p . (5.10)
46
5.4 Difusão dos Vórtices
A equação (5.4), que representa o transporte difusivo de vorticidade, possui solução
analítica para o campo de vorticidade de um vórtice inicialmente potencial. Esta solução
possui a seguinte forma adimensional
−Γ
= tr
et
tr 4Re2
4Re),(
πω , (5.11)
onde 22 yxr += , é a distância de um ponto qualquer do escoamento ao centro do
vórtice.
Para altos números de Reynolds, pode-se mostrar que a solução dada pela equação
(5.11) pode ser reproduzida pelo movimento aleatório de partículas de vorticidade, de
maneira similar ao movimento Browniano molecular. Com base na teoria de movimento
Browniano desenvolvida por Einstein, CHORIN (1973) propôs o assim chamado Método
do Avanço Randômico, que simula o transporte difusivo de vorticidade para qualquer
escoamento viscoso governado pela equação (5.4) através do movimento aleatório de uma
nuvem de vórtices. Quando o número de vórtices tende a infinito, este movimento produz
uma solução que tende à solução exata da equação (5.4). LEWIS (1991) propôs uma
metodologia de cálculo para o Método do Avanço Randômico um pouco diferente da
metodologia de CHORIN (1973). A metodologia implementada nesse trabalho é a de
LEWIS (1991), a qual é descrita abaixo.
Interpretando a função ω(r,t) como a probabilidade de um vórtice se encontrar em
uma dada região em um dado instante de tempo, o deslocamento difusivo dos vórtices ∆r
pode ser calculado invertendo-se a equação (5.11), o que fornece as equações em
coordenadas polares
∆
=∆P
tr 1lnRe4 e Qπθ 2=∆ , (5.12a, b)
47
onde P e Q são números randômicos entre 0 e 1 retirados de uma distribuição de
probabilidade estatisticamente uniforme. Em coordenadas cartesianas, tem-se
( )( ) iiD
D
rr
yx
=
θ∆∆θ∆∆
∆∆
sencos
. (5.13)
5.5 Avanço Temporal
Resolvidas as etapas convectivas e difusivas, é preciso avançar no tempo
atualizando a posição de cada vórtice da nuvem. Como o método do avanço randômico
fornece diretamente o deslocamento para cada vórtice, é necessário apenas integrar o
campo de velocidade convectivo, equação (5.6). Dois diferentes esquemas de avanço são
utilizados em geral: o avanço de Euler de primeira ordem e o avanço de Adams-Bashforth
de segunda ordem. As equações para esses esquemas são dadas, respectivamente, por
ttvtu
yx
iiC
C ∆∆∆
=
)()(
, (5.14)
tttvtvttutu
yx
iiC
C ∆∆∆
∆∆
−+−+
=
)(5,0)(5,1)(5,0)(5,1
. (5.15)
Como o esquema de Adams-Bashforth (RUGGIERO & LOPES, 1997) necessita do
armazenamento do campo de velocidade do passo anterior, o primeiro passo sempre é feito
com o esquema de Euler. Esquemas de predição-correção, tais como Runge-Kutta, não
foram adotados, pois o cálculo de mais um campo de velocidade da etapa convectiva em
um mesmo instante de tempo torna-se excessivamente dispendioso computacionalmente,
como citado em GUEDES (2003).
Calculados os deslocamentos totais, a posição da nuvem de vórtices pode ser
atualizada através de
+
+
=
++
D
D
C
C
ii yx
yx
tytx
ttyttx
∆∆
∆∆
∆∆
)()(
)()(
. (5.16)
48
5.6 Atualização das Condições de Contorno
Com a inserção de vórtices na esteira, suas velocidades devem passar a ser
consideradas para a imposição das condições de contorno na superfície do corpo. Como os
vórtices gerados possuem intensidade conhecida, sua contribuição é adicionada ao vetor
independente da equação matricial
iik
M
k v
vkii
ip
p
vu
vu
nnn ∑=
Γ−
−=
1 ˆˆ
sencos
αα
, i = (1, 2, ... , N). (5.17)
Desta forma, a cada passo de tempo, o vetor independente deve ser atualizado
mediante a nova configuração da esteira. Vale ressaltar que a matriz de influência
permanece inalterada para toda a simulação, sendo necessário alterar apenas o vetor
independente e resolver o novo sistema matricial.
Outra alteração importante no sistema de equações (4.23) é a remoção da condição
de Kutta (4.7b), que é válida para regime permanente. O escoamento aqui simulado ocorre
em regime transiente, porque se inicia do repouso. De modo a obedecer ao teorema de
Kelvin, a circulação inicial de todo o escoamento deve ser conservada, e, como o
escoamento parte do repouso, pode-se escrever
∑ ∑∫= =
+ =Γ+∆+
=Γ+N
j
M
kkj
jjEsteira
Corpo
ds1 1
1 02γγ
γ . (5.18)
Como no primeiro passo de tempo ainda não há vórtices na esteira, a equação se
torna
∑=
+ =∆+N
jj
jj
1
1 02γγ
(5.19)
e, como os vórtices gerados são obtidos a partir da vorticidade do corpo, cada etapa de
geração de vórtices também somará zero. Portanto esta equação substitui a equação (4.7b)
49
e é valida para todo o tempo, devendo ser adicionada ao sistema. Para facilitar sua
implementação computacional, a equação (5.19) pode ser escrita na forma
0222 1
2
11
1 =∆
+∆+∆
+∆
+=
−∑ NN
i
N
i
ii γγγ . (5.20)
O novo sistema, a ser resolvido para cada passo de tempo, possui a forma
( ) ( )
=
∆∆+∆∆+∆∆ +−
+
+
+
02222
2
1
1
2
1
1211
1,,2,1,
1,2,22221
1,1,11211
tN
t
t
N
N
NNN
NNNNNN
NN
NN
b
bb
aaaa
aaaaaaaa
γγ
γγ
, (5.21)
sendo que vetor bt na equação (5.21) deve ser atualizado a cada passo de tempo.
5.7 Tratamento de Vórtices que Penetram no Interior do Corpo
Ocasionalmente, devido ao avanço randômico, as linhas de corrente com curvas
muito acentuadas próximas à superfície do corpo ou aos erros numéricos inerentes ao
método, alguns vórtices penetram na superfície do corpo. Há, então, a necessidade de
eliminá-los ou refleti-los de volta para o escoamento.
O esquema de eliminação, proposto por CHORIN (1973), embora seja uma forma
aproximada de modelar a destruição de vorticidade na fronteira, reduz a resolução da
esteira e a circulação eliminada deve ser compensada no passo seguinte, por conta do
Teorema de Kelvin (LEWIS, 1991).
O esquema de reflexão garante que todos os vórtices gerados permaneçam no
escoamento, mantendo a resolução da esteira, mas acarretanto em um tempo
computacional maior. Uma aproximação satisfatória para o processo de reflexão consiste
em espelhar o vórtice para fora do corpo, rebatendo sua distância normal ao painel mais
próximo (Figura 5.3), conforme sugerido em LEWIS (1991), por exemplo.
50
5.8 Cálculo das Cargas Aerodinâmicas
A equação de Navier-Stokes, equação (3.2), quando avaliada na superfície do
corpo, na direção de sua superfície (s), assume a forma
sp
tqs
∂∂
−=∂
∂21 . (5.22)
Através da análise do escoamento potencial, a velocidade em cada ponto da
superfície é dada por jjjjtq γγγ =+= + 2)()( 1 . Utilizando esta informação e
discretizando a equação (5.22), pode-se escrever
tp ii
i ∆∆
=∆γ2 . (5.23)
A equação (5.23) expressa o diferencial de pressão em cada ponto de controle.
Integrando esta expressão, tem-se
∑=
∆+=i
nni ppp
10 . (5.24)
painéis vórtice no interior
do corpo
vórtice refletido
Figura 5.3 – Reflexão de vórtices
d
d
51
A pressão de referência p0 deve ser encontrada fazendo-se inicialmente p0 = 0 e
determinando o diferencial de pressão máximo no aerofólio, ps. Sendo assim, basta
adicionar (1 – ps) aos valores de pi previamente calculados para determinar a pressão
adimensional ao longo do corpo (LEWIS, 1991). Como a adimensionalização para pressão
foi feita com base na massa específica do fluido e na velocidade do escoamento incidente,
a pressão calculada corresponde ao coeficiente de pressão, definido pela equação (3.17).
Logo,
)1( sipi ppC −+= . (5.25)
Integrando Cp ao longo da superfície, através da regra do retângulo, é possível
calcular os coeficientes de arrasto e sustentação para cada instante de tempo. Assim,
na direção do escoamento, tem-se o coeficiente de arrasto:
( )∑=
+∆=N
iiipid CC
1sen βα ; (5.26)
na direção perpendicular ao escoamento, tem-se o coeficiente de sustentação:
( )∑=
+∆=N
iiipid CC
1cos βα . (5.27)
As cargas aerodinâmicas expressas pelas equações (5.26) e (5.27) são calculadas a
cada passo de tempo, podendo ser integradas no tempo para se determinar seus valores
médios. Esses valores são aqueles que podem ser comparados aos valores experimentais.
52
CAPÍTULO 6: IMPLEMENTAÇÃO NUMÉRICA
Os fundamentos do Método de Vórtices Discretos e do Método dos Painéis foram
descritos nos capítulos anteriores. A seguir, alguns aspectos adicionais necessários para
implementação numérica do método são apresentados. O método foi implementado
utilizando a linguagem FORTRAN, devido a sua robustez e grande aplicação na
comunidade científica para simulações numéricas de elevada complexidade. Todas as
variáveis e parâmetros reais foram declarados com dupla precisão, para minimizar erros de
arredondamento e truncamento.
6.1 Solução do Sistema de Equações Lineares
A solução do sistema de equações algébricas lineares muitas vezes é um ponto
crucial no método numérico, principalmente para métodos eulerianos, onde cada elemento
da malha representa uma equação algébrica linear a ser resolvida. No caso do Método de
Vórtices Discretos, por ser a abordagem lagrangiana, pode-se destacar como uma etapa
crucial do algoritmo o cálculo do campo de velocidade sobre cada vórtice. Porém, a
utilização do Método dos Painéis como parte do algoritmo do Método de Vórtices requer a
solução de um sistema linear a cada passo de tempo.
Conforme apresentado na seção 5.5, os coeficientes de influência dependem apenas
da geometria do corpo, que é fixa no problema estudado, sendo, portanto, calculados uma
única vez. A cada passo de tempo, basta atualizar o vetor de constantes e resolver
novamente o sistema de equações algébricas lineares. Essa característica particular do
modelo viabiliza adotar um esquema que inverta a matriz de influência uma única vez e a
cada passo de tempo seja realizada a operação
γbA =−1 . (6.1)
A inversão da matriz é realizada pela rotina DLINRG da biblioteca IMSL que faz
parte do pacote do compilador Compaq Visual Fortan Professional Edition 6.6C. Essa
53
rotina realiza a fatoração LU da matriz A, calcula a inversa das matrizes triangulares L-1 e
U-1 e a matriz inversa é obtida por
111 −−− = ULA . (6.2)
6.2 Geração de Números Randômicos
Os números randômicos obtidos de uma distribuição de probabilidade constante
requeridos pelo Método do Avanço Randômico (seção 5.3) foram calculados através da
rotina DRAND, uma rotina do compilador FORTRAN para geração de números
randômicos de dupla precisão entre 0 e 1.
6.3 Raio do Núcleo do Vórtice de Lamb
Conforme descrito na seção 3.6, utiliza-se neste trabalho o modelo de Lamb para
limitação do campo de velocidade de um vórtice. No entanto, o cálculo das velocidades
induzidas pelos vórtices utilizando equação (3.13), provoca grande esforço computacional,
devido ao cálculo do termo exponencial. Assim, a equação (3.13) do vórtice de Lamb, é
utilizada somente quando necessário, como descrito a seguir.
A velocidade induzida por um vórtice, segundo o modelo de Lamb, é
( )
−−
Γ= 2
2
exp12
,σπr
rtrq , na direção θ , (6.3)
onde ( ) ( )222oo YYXXr −+−= e Re4t=σ . O valor de r para o qual a velocidade
tangencial é máxima pode ser encontrado derivando-se a Eq. (5.27) e igualando a zero, o
que permite escrever
σ12091,1max =r . (6.4)
54
Com base na equação (6.4), adota-se como definição para o raio do núcleo do
vórtice o valor max0 2r=σ . Com esta definição, a velocidade avaliada em 0σ=r pela Eq.
(6.3) fornece uma diferença em relação ao campo de velocidade de um vórtice pontual (Eq.
2.10) de apenas 0,6%, como proposto por MUSTTO (1998). Pode-se, então, calcular a
velocidade induzida por um vórtice da seguinte forma:
• ( )r
trqπ2
, Γ= se 0σ≥r ; (6.5a)
• ( )
−−
Γ= 2
0
2
02572,5exp12
,σπr
rtrq se 0σ<r . (6.5b)
Para encontrar o valor do raio do núcleo, leva-se em consideração que a região
próxima ao corpo é a região onde a viscosidade possui maior influência, devendo ser
considerada em toda região entre o vórtice e o corpo. Para tal, o vórtice é gerado tal que
raio de seu núcleo toque a superfície do corpo no ponto de controle do painel
correspondente, conforme a Figura 6.1, desta forma,
εσ =0 . (6.6)
6.4 Velocidades Auto-Induzidas na Nuvem de Vórtices
Durante o cálculo da etapa convectiva do método, é necessário calcular a
velocidade induzida em cada vórtice da esteira por todos os demais componentes da
nuvem. Esta etapa é a mais lenta do algoritmo, pois requer o cálculo da interação mútua de
0σ
painéis pontos de controle
vórtices gerados
Figura 6.1 – Raio do núcleo dos vórtices
0σ
0σ
55
todos os vórtices da nuvem. Este cálculo que provoca um número de operações da ordem
de M2 a cada passo de tempo, onde M é o número de vórtices na esteira.
Os algoritmos de expansão em Multipolos propostos por BARNES & HUT (1986),
GREENGARD & ROKHLIN (1987) ou CARRIER et al. (1988) aceleraram a etapa
convectiva, mas não são utilizados neste trabalho. Desta forma, os erros inerentes apenas
ao Método de Vórtices podem ser isolados nas simulações.
Uma alternativa simples para reduzir o tempo computacional, proposta por LEWIS
(1991) e utilizada por ALCÂNTARA PEREIRA (1999), MUSTTO (1998, 2004) e
GUEDES (2003), pode ser implementada observando-se as expressões para velocidade
induzida de um vórtice j sobre um vórtice i
( )( ) ( )222
jiji
jijij
YYXX
YYU
−+−
−Γ=
π na direção X, (6.7a)
( )( ) ( )222
jiji
jijij
YYXX
XXV
−+−
−Γ−=
π na direção Y. (6.7b)
Analisando estas expressões é possível notar que, a menos da intensidade do
vórtice, a velocidade induzida por um vórtice i sobre um vórtice j é o simétrico do valor da
velocidade induzida por um vórtice j sobre um vórtice i, ou seja
ijiji UU ˆΓ−= e ijiji VV ˆΓ−= , (6.8a,b)
onde novamente o símbolo (^) denota o cálculo de velocidades com intensidade unitária.
Desta forma, basta calcular a velocidade induzida do vórtice j sobre o vórtice i que sua
recíproca já está automaticamente calculada. Caso os vórtices estejam afastados um do
outro por uma distância menor que o seu núcleo, basta aplicar as equações corrigidas da
seção (6.3). Esta simplificação reduz significativamente o tempo computacional das
simulações. Uma outra alternativa consiste em promover a coalescência dos vórtices
distantes do corpo (LEWIS,1991), porém este critério não foi utilizado neste trabalho, por
introduzir parâmetros numéricos adicionais e alguma perda de resolução da esteira.
56
6.5 Campo de velocidade nas Proximidades do Corpo
Quando as condições de contorno são atualizadas, isto é, o vetor b é recalculado
para uma nova configuração da esteira, alguns vórtices podem estar muito próximos do
ponto de controle de algum painel do corpo, o que implica na indução de uma velocidade
muito alta sobre este ponto de controle. Da mesma forma, durante o cálculo da etapa
convectiva, a velocidade induzida pelo corpo sobre vórtices muito próximos a ele podem
ser elevadas, sobretudo nas proximidades dos pontos nodais, onde há uma descontinuidade
na derivada da vorticidade (vide capítulo 4). Essas velocidades elevadas representam um
erro numérico, decorrente da discretização do corpo e da discretização da vorticidade em
sua superfície. Nesta seção serão apresentadas algumas formas de contornar estas
dificuldades numéricas.
Similarmente ao sugerido por LEWIS (1991), para o cálculo da velocidade
induzida do vórtice sobre o ponto de controle do painel, é recomendado utilizar o esquema
de sub-painéis, quando o vórtice estiver a uma distância ao ponto de controle mais
próximo(rjk) menor que o comprimento do painel k, k∆ .
Já para velocidade induzida do corpo sobre o vórtice, os vórtices muito próximos
ao corpo devem receber tratamento especial, quando o vórtice estiver a uma distância ao
ponto de controle mais próximo (rjk) menor que 40% do comprimento do painel k, k∆4,0 .
Neste caso, não faz sentido utilizar sub-panéis para velocidade induzida do corpo sobre o
vórtice, pois a distribuição de vorticidade sobre o painel já é contínua e a velocidade
induzida é o resultado de uma integração ao longo do painel (ver capítulo 4).A Tabela 6.1
resume os parâmetros envolvidos no cálculo especial de velocidades induzidas nas
proximidades do corpo.
Tabela 6.1 – Tratamento especial para vórtices próximos ao corpo
Tratamento Especial Etapa de Cálculo Critério
Sub-painéis Velocidade Induzida pelos
vórtices sobre o corpo kjkr ∆<
Método da Imagem Velocidade Induzida pelo
corpo sobre os vórtices kjkr ∆< 4,0
57
6.5.1 Utilização de Sub-Painéis
Nas simulações aqui realizadas é adotado o modelo sugerido por LEWIS (1991),
onde o painel j em análise é dividido em Nsub sub-painéis iguais, conforme mostra a Figura
6.2. A velocidade induzida sobre este painel é calculada através da média aritmética das
velocidades induzidas pelos vórtices próximos ao painel j sobre os sub-painéis, ou seja
∑=
Γ=
subN
p pj
jsubkjv
v
VU
NVU
1 ˆˆ1 . (6.9)
6.5.2 Método da Imagem
Neste caso, a velocidade induzida pelo painel cujo ponto de controle encontra-se
mais próximo do vórtice é substituída pela velocidade induzida pela imagem do vórtice no
interior do corpo, com intensidade de sinal contrário à intensidade do vórtice original.
Desta forma, é possível obter-se um controle muito melhor dos valores da velocidade
induzida pelo painel, pois a imagem força o vórtice a se mover paralelamente ao painel.
Uma representação esquemática desta situação pode ser vista na Figura 6.3.
( Xj , Yj )jΓ
rjk
12
p
Nsub
js∆
Figura 6.2 – Utilização de Sub-painéis, LEWIS (1991)
sub
j
Ns∆
58
6.6 Reflexão dos Vórtices
Quando é necessário refletir um vórtice para o interior ou exterior do corpo, é
necessário calcular sua nova posição considerando a inclinação do seu painel mais
próximo. No esquema representado na Figura 6.4, o vórtice a ser refletido está no interior
do corpo e conseqüentemente será refletido para fora dele. A partir de correlações
trigonométricas e da Figura 6.4 é possível deduzir a seguinte expressão
±
=
j
j
v
v
r
r
YX
YX
ββ
δcossen
2 , (6.10)
onde o sinal positivo é utilizado no caso de reflexão para o exterior do corpo e o sinal
positivo quando a reflexão for para o interior do corpo. O valor de ( ) ivc YY βδ cos−=
representa a distância normal do vórtice até o painel mais próximo. O valor de Yc deve ser
calculado a partir da reta que define o painel j, ou seja
jvjc kXkY ′+= (6.11)
onde jj
jjj XX
YYk
−
−=
+
+
1
1 e jjjj XkYk −=′ (6.12a,b)
painéis Imagem
do vórtice
vórtice original
Figura 6.3 – Método da Imagem
rjk
d
d
Ponto de controle jΓ
jΓ−
59
6.7 Parâmetros Numéricos
Embora haja um certo grau de arbitrariedade na escolha dos valores dos parâmetros
numéricos, existem alguns critérios objetivos baseados na Física que governa o
escoamento que podem ser usados para se estimar os valores adequados necessários para a
obtenção de simulações eficientes. Os critérios adotados neste trabalho são discutidos a
seguir.
Em primeiro lugar, a escolha do número de painéis (N) deve ser feita
cuidadosamente, pois este número também representa o número de vórtices nascentes que
são gerados na camada limite do corpo e liberados para a esteira a cada passo de tempo.
Quanto maior N, maior será o número M de vórtices presentes na esteira, pois M = NTxN,
onde NT é o número de passos de tempo da simulação. Conseqüentemente, quanto maior
M, melhor é a resolução do campo de vorticidade do escoamento. Teoricamente, quanto
maior M, melhor é a simulação da dinâmica da vorticidade na camada limite e na esteira do
corpo. Porém, à medida que M aumenta, maior se torna o tempo de processamento de cada
simulação, pois o contador de operações de ponto flutuante por passo de tempo é
proporcional a M2, devido à etapa convectiva do algoritmo. Por outro lado, se N é pequeno
a simulação será muito rápida, mas os resultados não terão uma boa exatidão devido à
baixa resolução do campo de vorticidade do escoamento.
( Xj , Yj )
( Xj+1 , Yj+1 )
( Xr , Yr )
( Xv , Yv )
jβ
δ
δ
( Xv , Yc )
painel j
Figura 6.4 – Reflexão de um vórtice para fora do corpo
Interior do corpo
60
Testes preliminares para um escoamento potencial comparado à solução analítica
para o aerofólio de Van de Vooren indicam valores de N entre 200 e 300 painéis para a
obtenção de resultados com ótimo grau de exatidão. Para 300 painéis, por exemplo, chega-
se a obter resultados para o coeficiente de sustentação com erro da ordem de 0,005%
(PEREIRA et al. ,2004).
Testes preliminares para as simulações com o Método de Vórtices sugerem que os
mesmos valores de N fornecem resultados com bom grau de exatidão com um tempo de
processamento (tempo de CPU) razoável (até 28 horas de simulação), em uma plataforma
com processador Pentium IV. A partir desses testes, adotou-se o valor N = 300 para todas
as simulações realizadas neste trabalho.
A distância de geração, suposta constante ao longo do corpo, é calculada com base
na espessura de deslocamento da camada limite, onde ∞≈≈ cQνδε 1 . Na sua forma
adimensional, a distância de geração é estimada por
21
Re −= εε k , (6.13)
onde εk é uma constante, da ordem de 2, ajustada com base em testes numéricos e
νcQ∞=Re é o número de Reynolds global do escoamento.
O incremento de tempo, ∆t, é escolhido de tal forma que o deslocamento de um
ponto material sob ação do escoamento uniforme não exceda aproximadamente o
comprimento de um painel, ou seja ∆≈∆tQ . Para um painel de comprimento médio
Nc2≈∆ , pode-se escrever uma estimativa para o incremento de tempo, em variáveis
adimensionais, como
Nk
t t2=∆ , (6.14)
onde kt é uma constante de ordem 1, mas geralmente menor que a unidade quando o
esquema de Euler é utilizado e maior que a unidade quando o esquema de Adams-
61
Bashforth é utilizado para o avanço dos vórtices no tempo. Valores adequados para
produzir boas simulações exigem alguma experimentação numérica, mas pode-se adiantar
que a faixa 0,01 ≤ ∆t ≤ 0,05 é adequada para a produção de bons resultados finais.
6.8 Particularidades do Algoritmo
Algumas particularidades da implementação numérica são muito importantes para
o desempenho geral do algoritmo utilizado e merecem alguns comentários.
O tratamento especial de vórtices muito próximos do corpo é de vital importância
para qualidade da simulação, sobretudo devido ao fato de que a geometria em estudo é a de
um aerofólio, que possui uma canto vivo no bordo de fuga e a uma região de curvatura
acentuada no bordo de ataque. Os testes realizados comprovaram a necessidade de
implementação das estratégias da seção 6.5, sem as quais a simulação torna-se instável e
fornece resultados sem significado físico algum.
Uma outra particularidade do algoritmo aqui utilizado é a condição de conservação
de circulação, representada pela equação 5.20, que pode, a princípio, parecer contradizer o
teorema de Kutta-Jukowsky (KUNDU & COHEN, 2002), uma vez que a soma da
vorticidade criada sobre o aerofólio é nula, mesmo que este possua força de sustentação.
Esta aparente contradição pode ser explicada pelo fato de que o somatório da vorticidade
sobre o corpo representa a circulação ao redor do corpo no caso potencial, quando a
camada-limite é suposta ter uma espessura desprezível. Na metodologia aqui empregada, a
vorticidade gerada na parede movimenta-se por convecção e difusão, formando uma fina
camada de vorticidade ao longo da superfície do corpo, de espessura não nula,
similarmente a um escoamento real. A circulação ao redor do corpo, incluindo a
vorticidade nessa camada ao seu redor é que deve ser considerada como circulação ao
redor do corpo e esta é não nula, mostrando que não há nenhuma contradição no modelo
empregado.
Pode-se ainda destacar uma última particularidade do algoritmo, a qual se refere à
estimativa dos valores dos parâmetros numéricos inerentes ao Método de Vórtices aqui
utilizado. Experiência anterior com esta metodologia, assim como os testes preliminares já
realizados, demonstram que a escolha dos valores dos parâmetros numéricos é crucial para
62
a obtenção de uma boa simulação. A utilização dos critérios descritos acima para a
estimativa dos valores dos parâmetros numéricos sempre proporciona simulações com
resultados coerentes com a Física do problema. Simulações para a obtenção de resultados
finais requerem um refinamento desses valores visando aumentar a resolução do campo de
vorticidade do escoamento.
6.9 Rotinas Numéricas
Um grande número de rotinas numéricas programadas em linguagem FORTRAN
foi desenvolvido para a simulação de escoamentos externos, bidimensionais e
incompressíveis ao redor de aerofólios utilizando o algoritmo do Método de Vórtices
Discretos, enquanto que algumas outras, conforme citado ao longo deste capítulo, que já
estavam disponíveis na biblioteca IMSL do compilador FORTRAN utilizado, foram
incorporadas ao algoritmo aqui desenvolvido. Todas as rotinas desenvolvidas no código
numérico final estão descritas a seguir, em ordem alfabética. Essas rotinas estão definidas
dentro do programa principal, METVORT13.
CONVEC
Rotina para cálculo das velocidades conevctivas induzidas nos vórtices pela nuvem
de vórtices, pelo escoamento incidente e pelo corpo.
CORRENTE
Rotina para cáculo de velocidades em uma malha definida para traçar
posteriormente as linhas de corrente ao final da simulação.
DESLOC
Rotina para deslocar os vórtice segundo o esquema de avanço temporal escolhido
(Adams-Bashforth ou Euler) e adicionar o deslocamento difusivo.
DIFUS
Rotina para cáculo dos deslocamentos difusivos segundo o Método do Avanço
andômico.
63
FILME
Rotina para salvar periodicamente posições intermediárias da nuvem de vórtices e
docoeficiente de pressão para animações.
GEOMETRIA
Rotina para solicitação, via terminal, da escolha da geometria a ser analisada.
GRAVAR
Rotina para salvar os vetores, quando necessário.
GRIDNACA4
Rotina para cáculo dos pontos nodais, pontos de controle, inclinação, vetores
nomrais e tangentes dos painéis. Esta rotina pode ser substituída ou adaptada com
cerata facilidade para análise de outras geometrias.
INFOMETVORT
Rotina para apresentação de uma tela de identificação do programa antes do início
do cálculo. Esta rotina não possui parâmetros de entrada e saída.
INTPRESSAO
Rotina para cálculo dos carregamentos aerodinâmicos via intergração de pressão ao
longo do corpo.
MATSISTEMA
Rotina para construção e inversão da matriz dos coeficintes de influência.
NACA4
Rotina que calcula a coordenada Y, a partir de uma abscissa X, com base na função
que descreve o contorno do aerofólio.
NOMEARQ
Rotina que nomeia os arquivos a serem salvos do programa.
64
NOVOB
Rotina para atualização do vetor independente da equação matricial. Esta rotina
calcula o campo de velocidades induzido pela nuvem de vórtices e pelo escoamento
incidente sobre os pontos de controle do corpo.
REFLEX
Rotina para reflexão dos vórtices que penetram no corpo.
RESULTADOS
Rotina que salva os resultados finais do programa.
SIMDADOS
Rotina para leitura dos parâmetros numéricos da simulação, obtidos via um arquivo
de dados previamente construído eplo usuário.
SOLSIST
Rotina para solução do sistema linear. Esta rotina resume-se à multiplicação da
inversa da matriz pelo vetor independente da equação matricial, atualizada a cada
passo de tempo.
SUBPAN
Rotina para divisão do painel em sub-painéis, acionada pela rotina NOVOB, quando
o vórtice se aproxima do corpo segundo os critérios da seção 6.5.
VORTL
Rotina para o calcula das velocidades induzidas por um painel com distribuição de
vorticidade linear sobre um ponto qualquer do domínio. Utilizam principalmente as
equações 4.15a, 4.15b, 4.16a e 4.16b.
VORTLAMB
Rotina para cálculo das velocidades induzidas por um vórtice pontual, incluindo o
modelo de núcleo viscoso de Lamb (Equações 3.14 e 3.15).
65
VORTNASC
Rotina para detreminação da posição em que serão gerados novos vórtices. Com
estas posições são constantes ao longo de toda simulação, esta rotina é acionada um
única vez.
ROTINAS IMSL
- DLINRG: Rotina para inversão da matriz atraévs de fatoração LU, utilizando
dupla precisão.
- DMURRV: Rotina para multiplicação da matriz invertida pelo vetor
independente (solução do sistema).
- DRAND( ): Rotina para geração de números randômicos entre 0 e 1.
6.10 Fluxograma
Um fluxograma indicando o algoritmo do programa elaborado está representado na
Figura 6.5, fornecendo uma visão global do procedimento de cálculo com o qual foram
obtidos os resultados descritos a seguir.
66
Início
Entrada de Dados Escolha da geometria
Divisão da Geometria em painéis e cálculo de suas propriedades
Cálculo posição dos vórtices nascentes e posicionamento dos
primeiros vórtices na esteira
Cálculo e inversão da matriz de influência (A)
Solução do sistema de equações para t=0 ( γ = A-1 b)
Solução Potencial
Atualização do vetor de constantes (bt)
Solução do novo sistema de equações – nova distribuição de
vorticidade sobre o corpo
Cálculo da velocidade convectiva
Cálculo do deslocamento difusivo
Avanço no tempo – Cálculo da nova configuração da esteira (Adams-Bashforth ou Euler)
Reflexão dos vórtices que penetraram no corpo
Nascimento de novos vórtices
Último passo de tempo? N
S
Gravação dos Resultados finais
Fim
Figura 6.5 – Fluxograma do programa elaborado em linguagem FORTRAN
Cálculo das Cargas Aerodinâmicas
ttt ∆+←
67
CAPÍTULO 7: RESULTADOS E DISCUSSÕES
Uma vez estabelecidos todos os fundamentos teóricos e considerações numéricas
nos capítulos anteriores, foram realizados vários testes com o modelo numérico
desenvolvido, visando à validação do código computacional e determinação dos valores
dos parâmetros adequados para simulação do escoamento incompressível e bidimensional
ao redor das geometrias de interesse.
Os principais parâmetros de comparação entre os resultados obtidos e outros
resultados numéricos ou experimentais disponíveis na literatura são os coeficientes de
arrasto e sustentação. Embora esses coeficientes, calculados para cada passo de tempo,
sofram oscilações no tempo resultantes do processo de desprendimento de vorticidade da
camada limite, seus valores médios temporais podem ser avaliados para possibilitar uma
comparação quantitativa. O formato e a evolução temporal da esteira do corpo também
podem ser utilizados como importantes indicações da qualidade das simulações, além da
distribuição de pressão sobre a superfície do corpo.
Inicialmente, o programa desenvolvido foi testado para cilindros circulares, pois
esta é uma geometria extensamente explorada na literatura e de comportamento
característico bastante conhecido para esteira e carregamentos.
O comportamento da camada limite sobre o aerofólio também é investigado neste
capítulo, pela sua grande influência no escoamento.
Posteriormente foram realizados testes completos com aerofólios para diferentes
ângulos de ataque, tendo sido escolhidos os aerofólios NACA 0012 e NACA 4415. O
primeiro, pela abundância de resultados da literatura e para comparação com o trabalho de
AKBARI & PRICE (2003), que também utilizam o Método de Vórtices; o segundo para
investigar a capacidade do modelo numérico desenvolvido de simular escoamentos ao
redor de aerofólios com curvatura não nula, além deste aerofólio também possuir maior
espessura. O aerofólio NACA 4415 foi escolhido particularmente por ser utilizado em
projetos de pás de turbinas eólicas.
68
7.1 Teste Preliminar: Cilindro Circular
Com intuito de simular uma geometria mais simples e com resultados bem
conhecidos, utiliza-se inicialmente um cilindro circular com número de Reynolds igual a
105. Os resultados obtidos foram semelhantes aos obtidos por GUEDES (2003), mostrando
uma esteira com estruturas de vorticidade contra-rotativas do tipo de Von Karman. A
Figura 7.1 mostra a esteira do cilindro ao final da simulação. A representação da esteira
corresponde à posição adquirida pelos vórtices como resultado natural da simulação,
representando a região rotacional do escoamento. O número de Strouhal, que pode ser
interpretado como a frequência adimensional de emissão de vórtices na esteira, foi da
ordem de 0,2, conforme esperado (SCHLICHTING & GESTERN, 2000). O coeficiente de
arrasto médio obtido estava em torno de 1,89, acima do resultado experimental, da ordem
de 1,2 (BLEVINS,1984), conforme esperado para simulações bidimensionais. Este teste
preliminar foi realizado com intuito de validar o modelo de forma geral, verificando apenas
se as simulações possuem significado físico apreciável. Para uma análise detalhada sobre a
aplicação do Método de Vórtices associado ao Método dos Painéis para cilindros circulares
recomenda-se consultar o trabalho de GUEDES (2003).
7.2 Determinação dos Parâmetros Numéricos
Os testes preliminares foram necessários para ajustar os parâmetros numéricos
inerentes ao Método de Vórtices. Dentre as possibilidades, foram escolhidos os parâmetros
que possibilitariam a melhor representação física do escoamento, sem acarretar em tempos
computacionais muito longos.
Figura 7.1 – Esteira de um cilindro circular ao final da simulação [Esquema de Euler ; Re = 105 ; ε = 0,01; ∆t = 0,1 ; N = 200; 400 passos de tempo]
69
O número de painéis para discretização de aerofólios e consequente número de
vórtices gerados a cada passo de tempo foi fixado em 300, devido a melhor representação
do campo de velocidade e pressão nas proximidades do bordo de ataque e bordo de fuga no
escoamento potencial ao redor de um aerofólio. A Figura 7.2, semelhante à Figura 4.8,
mostra a distribuição do coeficiente de pressão nas regiões do bordo de ataque e bordo de
fuga de um aerofólio de Van de Vooren de espessura 15% e ângulo de ataque de 10o
calculada com 300 painéis e comparada à solução analítica. Como pode ser visto na figura,
há uma boa concordância entre os resultados numéricos e a solução analítica quando N =
300. Por este motivo, utilizou-se 300 painéis em todas as simulações para aerofólios neste
trabalho. Para detalhes sobre a geometria do aerofólio de Van de Vooren, pode-se
consultar KATZ & PLOTKIN (1991).
0 0.01 0.024.5
4.6
4.7
4.8
4.9
5
X
-Cp
LE, N = 50
0.98 0.99 1-1
-0.9
-0.8
-0.7
-0.6
-0.5
-0.4
-0.3
X
-Cp
TE, N = 50
0 0.01 0.024.5
4.6
4.7
4.8
4.9
5
X
LE, N = 80
0.98 0.99 1-1
-0.9
-0.8
-0.7
-0.6
-0.5
-0.4
-0.3
X
TE, N = 80
0 0.01 0.024.5
4.6
4.7
4.8
4.9
5
X
LE, N = 100
0.98 0.99 1-1
-0.9
-0.8
-0.7
-0.6
-0.5
-0.4
-0.3
X
TE, N = 100
0 0.01 0.024.5
4.6
4.7
4.8
4.9
5
X
LE, N = 300
0.98 0.99 1-1
-0.9
-0.8
-0.7
-0.6
-0.5
-0.4
-0.3
X
TE, N = 300
Resultados NumericosSoluçao Analitica
Figura 7.2 – Distribuição de pressão nas vizinhanças do bordo de ataque (LE) e bordo de fuga (TE) do aerofólio com 50, 80, 100 e 300 painéis.
Cp
Cp
70
Todas as simulações deste ponto em diante foram realizadas com o esquema de
Adams-Bashforth. Este esquema foi escolhido pois seu erro é da ordem de ∆t2, ao contrário
de esquemas de primeira ordem, como o esquema de Euler, que possui erro da ordem de
∆t. Já o esquema Runge-Kutta foi preterido devido ao elevado tempo computacional
necessário para cálculo do campo de velocidade duas vezes para cada passo de tempo
(GUEDES, 2003). O incremento de tempo foi calculado utilizando kt = 3,75, de forma a
evitar que fossem necessários tempos muito longos para uma simulação completa. Sendo
assim, através da equação (6.14), calcula-se
∆t = 0,025. (7.1)
Como quanto maior o número de Reynolds, maior o número de graus de liberdade
do escoamento (SILVEIRA NETO, 2004), ou seja, mais rica em vorticidade são a camada
limite e a esteira. Assim, as simulações foram realizadas com Reynolds da ordem de 105,
pois a modelagem empregada não é apropriada para escoamentos totalmente turbulentos de
números de Reynolds muito altos com a resolução numérica empregada. Para valores do
número de Reynolds muito altos é necessário de se aumentar substancialmente o número
de vórtices presentes no escoamento ou utilizar uma modelagem de turbulência para as
pequenas escalas. Segundo KUNDU & COHEN (2002), para Reynolds maiores que 106,
ocorre a transição da camada limite sobre o aerofólio de laminar para turbulenta.
Utilizando kε = 2,0 e Re = 1,7x105 na equação (6.13), chega-se então ao seguinte valor para
a distância de geração de vórtices:
ε = 0,005. (7.2)
Este valor é o mesmo empregado para o raio do núcleo dos vórtices, conforme
descrito na seção 6.3. Quando a geração de vórtices é feita a uma distância mais próxima
ao corpo, as simulações se mostram instáveis e não convergem ao longo do tempo,
fornecendo resultados desprovidos de significado físico.
Quanto ao número de sub-painéis, foram utilizados sempre 5 sub-painéis, conforme
sugerido em LEWIS (1991), uma vez que foi comprovado em testes preliminares que um
71
número maior que este não representa nenhum ganho adicional de atenuação das
velocidades induzidas sobre os pontos de controle.
O número de passos no tempo é idealmente escolhido para que a simulação atinja
um regime periódico permanente, com oscilações periódicas dos carregamentos e esteira.
Porém este tipo de comportamento só é atingido totalmente para tempos longos para altos
ângulos de ataque, ocasionando enorme esforço computacional. Assim, a maioria das
simulações foi limitada a 400 passos de tempo, o que equivale, com o ∆t escolhido, a um
tempo adimensional total de 10 e 120.000 vórtices presentes na esteira ao final da
simulação. Este tempo total é suficiente para a avaliação dos coeficientes de força médios
(no tempo). A Tabela 7.1 resume os dados apresentados nesta seção que serão utilizados
em todas as simulações deste ponto em diante, exceto quando claramente identificado.
Tabela 7.1 – Parametros Numéricos Empregados
Parâmetro Símbolo Valor Empregado
Número de Reynolds Re 1,7 x 105
Incremento de tempo ∆t 0,025
Número de painéis N 300
Distância de geração ε 0,005
Raio do núcleo dos vórtices σ0 0,005
Número de Passos no tempo NT 400
7.3 Aerofólios Estudados
Devido à grande aplicação de aerofólios em diversos equipamentos da indústria em
geral, este tipo de corpo esbelto foi escolhido para as simulações deste trabalho. As
geometrias dos aerofólios empregados, o NACA 0012 e o NACA 4415, estão apresentadas
na Figura 7.3. O aerofólio NACA 0012 é simétrico com 12% de espessura, enquanto que o
aerofólio NACA 4415 é assimétrico, com 15% de espessura e 4% de curvatura máxima
localizada a 0,4c do bordo de ataque. A descrição detalhada das funções matemáticas que
definem o contorno dos aerofólios encontra-se no Apêndice B. A escolha de um aerofólio
simétrico e outro assimétrico permite avaliar melhor o desempenho do modelo numérico
desenvolvido neste trabalho para simular corpos aerodinâmicos em geral.
72
Ao contrário do cilindro circular, o escoamento ao redor de um aerofólio depende
fortemente do ângulo de ataque, podendo apresentar padrões de escoamento totalmente
diferenciados para pequenas variações deste ângulo. Um fenômeno de particular interesse é
o “estol”, denominação utilizada quando ocorre a separação massiva da camada limite
sobre a superfície do aerofólio em uma região a montante do bordo de fuga, originando
escoamento reverso e desprendimento de vórtices. Como conseqüência, ocorre uma queda
no coeficiente de sustentação e um aumento significativo do coeficiente de arrasto. A
análise da região de “estol” é extremamente importante, pois neste ponto ocorre um valor
máximo do coeficiente de sustentação na região de baixo arrasto. Embora o coeficiente de
sustentação volte a crescer depois da região de “estol”, o coeficiente de arrasto cresce
significativamente, não sendo mais vantajoso trabalhar com aerofólios nesta situação.
Para cada um dos aerofólios estudados, simulações foram realizadas para diversos
ângulos de ataque, de forma a possibilitar a identificação da região de “estol” e demais
fenômenos físicos envolvidos. Cada simulação produz os seguintes resultados diretos do
programa: os coeficientes de arrasto e sustentação ao longo do tempo, os coeficientes de
pressão ao longo da superfície do aerofólio, as disposições da esteira para diversos
instantes de tempo e as linhas de corrente ao final de cada simulação. O modelo fornece
naturalmente o campo de velocidade em cada ponto do escoamento e em cada instante de
tempo, caso seja de interesse.
Devido à discretização do campo de vorticidade, os carregamentos sobre os
aerofólios sofrem flutuações no tempo significativas, sendo, portanto, avaliados por médias
temporais, para facilitar sua análise e comparação com outros resultados numéricos e
experimentais.
NACA 0012
NACA 4415
Figura 7.3 – Geometria dos aerofólios empregados nas simulações
73
7.4 Camada Limite sobre um Aerofólio
Visando investigar o comportamento da camada limite sobre um aerofólio, foram
realizadas simulações utilizando o aerofólio NACA 0012, para ângulos de ataque de 0o e
10o e os parâmetros listados na seção anterior. O número de passos no tempo, foi reduzido,
pois não deseja-se capturar o caráter oscilatório da esteira, porém somente a formação da
camada limite. Portanto, com base em experimentação numérica foram utilizados 80
passos de tempo, correspondendo a um tempo adimensional t = 2,0.
A Figura 7.4, obtida a partir da simulação com ângulo de ataque nulo apresenta a
velocidade na direção X nas proximidades do aerofólio ao final da simulação. Como
esperado, a velocidade cresce a partir de um valor próximo de zero, até um valor
aproximadamente constante, correspondente à velocidade fora da camada limite. Nesta
mesma figura observa-se que foram escolhidas 10 regiões diferentes do aerofólio, 5 na
parte superior e 5 na parte inferior. Em cada uma destas regiões a velocidade é calculada
em um ponto de controle e em pontos com a mesma coordenada X que o ponto de controle.
Como este aerofólio é simétrico e a simulação foi feita com ângulo de ataque nulo, é
esperado um comportamento similar das regiões superior e inferior do aerofólio, para uma
mesma distância ao bordo de ataque.As Figuras 7.5a até 7.5e, apresentam o perfil de
velocidade mais detalhadamente para cada coordenada X.
74
Figura 7.4 – Perfil de velocidade próxima ao aerofólio
(a) (b) (c) (e) (d)
0 0.5 10
0.005
0.01
0.015
0.02
0.025
u
Y
0 0.5 10
0.005
0.01
0.015
0.02
0.025
u
Y
0 0.5 10
0.005
0.01
0.015
0.02
0.025
u
Y
0 0.5 10
0.005
0.01
0.015
0.02
0.025
u
Y
0 0.5 10
0.005
0.01
0.015
0.02
0.025
u
Y
Superficie inferiorSuperficie superiorPerfil de Blasius
Figura 7.5 – Detalhe dos perfis de velocidade próximos ao aerofólio
(a) (b) (c)
(d) (e)
Perfil de Blasius Superfície superior Superfície inferior *
Y Y Y
Y Y
u u u
u u
75
Os perfis de velocidade são apresentados juntamente com o perfil de velocidade de
camada limite laminar de Blasius para uma placa plana infinita e sem espessura
(SCHLICHTING & GERSTEN, 2000). Nas Figuras 7.5a e 7.5b, regiões mais próximas do
bordo de ataque, pode-se observar uma razoável concordância dos resultados, devido à fina
geometria do aerofólio. O valor da velocidade calculada fora da camada limite é maior que
o valor teórico de Blasius devido à espessura finita do aerofólio, que se contrapõe ao
resultado de Blasius, válido para uma placa plana sem espessura.
As Figuras 7.5c até 7.5e apresentam as demais regiões onde a camada limite foi
calculada. Nestas regiões, a comparação com o perfil de Blasius já não é mais válida, pois
a geometria do aerofólio provoca variações significativas no perfil de velocidades. Nas
últimas regiões (Figura 7.5e), é possível identificar a presença de escoamento reverso,
característica que indica o início do processo de separação da camada limite laminar em
simulações bidimensionais, (SCHLICHTING & GERSTEN, 2000). Este efeito é esperado
nas proximidades do bordo de fuga, devido à corda finita do aerofólio. Os demais perfis
indicam que a camada limite apresenta oscilações em seu perfil resultantes do mecanismo
oscilatório de desprendimento de vorticidade para a formação da esteira, tendo esses perfis
gradientes significativos de velocidade na direção perpendicular ao escoamento, mesmo
em regiões distantes do bordo de fuga do aerofólio. Os perfis de velocidade mostrados
sugerem a ocorrência de separação prematura do escoamento a montante do bordo de fuga,
mesmo em pequena intensidade. Este efeito é comum em simulações que utilizam o
Método de Vórtices, como observado por SPARLAT & LEONARD (1981) e LEWIS
(1991). As simulações realizadas tentam minimizar este efeito com a utilização de um
grande número de vórtices nascentes (N = 300), visando melhorar a discretização do
campo de vorticidade. Porém, ainda assim, observa-se uma tendência de aparecimento de
estruturas de vorticidade maiores que as visualizadas em experimentos.
As Figuras 7.6 e 7.7 apresentam o mesmo conjunto de resultados, porém para um
ângulo de ataque de 10o. Nestas figuras é possível notar como o perfil de velocidade na
região inferior do aerofólio permanece com comportamento similar ao perfil de camada
limite laminar de Blasius, enquanto que a parte superior apresenta perfis de velocidade
com grandes regiões de escoamento reverso, sugerindo separação intensa do escoamento.
As regiões inflexionais na camada limite provocam o surgimento de estruturas de
vorticidade, similares a instabilidades de Kelvin-Helmholtz (SILVEIRA NETO, 2004).
76
Figura 7.6 – Perfil de velocidade próxima ao aerofólio
(a) (b) (c) (e) (d)
-0.5 0 0.5 10
0.005
0.01
0.015
0.02
0.025
u
Y
-0.5 0 0.5 10
0.005
0.01
0.015
0.02
0.025
u
Y
-0.5 0 0.5 10
0.005
0.01
0.015
0.02
0.025
u
Y
-0.5 0 0.5 10
0.005
0.01
0.015
0.02
0.025
u
Y
-0.5 0 0.5 10
0.005
0.01
0.015
0.02
0.025
u
Y
Superficie InferiorSuperficie SuperiorPerfil de Blasius
(a) (b) (c)
(d) (e)
Figura 7.7 – Detalhe dos perfis de velocidade próximos ao aerofólio
Perfil de Blasius Superfície superior Superfície inferior *
77
7.5 Aerofólio NACA 0012
Conforme já mencionado, as simulações para o aerofólio NACA 0012 são de
grande valia para validação do código computacional e observação dos fenômenos físicos
em análise devido à abundância de resultados disponíveis na literatura para este aerofólio.
Utilizando os parâmetros numéricos apresentados na seção 7.2, foram realizadas
simulações para o perfil NACA 0012 com ângulo de ataque variando de 0 a 75o. As
simulações foram realizadas com intuito de observar o “estol” sobre o aerofólio e, por este
motivo, no intervalo de 0 a 20o onde experimentalmente é observado o “estol”, o
incremento do ângulo de ataque foi menor. A partir de α = 20o os incrementos foram
maiores. Devido à grande quantidade de casos simulados, a análise que se segue está
centralizada em apenas alguns valores do ângulo de ataque que são representativos dos
fenômenos físicos em estudo. Os resultados de todas as simulações realizadas com esses
parâmetros estão listados no Apêndice C.
Para ângulos de ataque baixos, inferiores a 10o, pode-se notar que, similarmente a
observações experimentais, as simulações geram um escoamento razoavelmente aderido à
superfície, com a formação de regiões de separação somente nas proximidades do bordo de
fuga do aerofólio. A Figura 7.8 exemplifica este tipo de comportamento, apresentando a
evolução temporal dos vórtices que modelam a camada limite e a esteira na simulação para
um ângulo de ataque de 6o. Através do retrato instantâneo da posição dos vórtices, pode-se
observar que a esteira formada é fina e com pequenas formações de pares de vórtices
responsáveis por pequenas oscilações temporais nos carregamentos sobre o aerofólio.
Como conseqüência, a distribuição de pressão sofre pequenas oscilações temporais. Os
gráficos apresentados na Figura 7.8 representam a média temporal do coeficiente de
pressão a cada unidade de tempo adimensional. Nesta mesma figura também é possível
observar o caráter oscilatório da nuvem de vórtices, sobretudo na região da esteira, que
aumenta suavemente à medida que o tempo evolui. A partir de t = 4,0 a esteira não é
apresentada por completo, para favorecer a visualização do escoamento nas proximidades
do aerofólio.
Os enrolamentos de vórtices em forma de cogumelo que se deslocam
significativamente na direção perpendicular ao escoamento, originando uma esteira mais
78
espessa que o esperado, são supostamente oriundos da bidimensionalidade das simulações
utilizando Método de Vórtices, pela ausência de difusão na terceira dimensão.
0 0.5 1-1
0
1
2
x
-cp
0 0.5 1-1
0
1
2
x
-cp
0 0.5 1-1
0
1
2
x
-cp
0 0.5 1-1
0
1
2
x
-cp
0 0.5 1-1
0
1
2
x
-cp
t = 1,0
t = 4,0
t = 3,0
t = 5,0
t = 2,0
Figura 7.8 (a) – Esteira ao longo da simulação e coeficiente de pressão médio aolongo do aerofólio NACA 0012, de t = 1,0 até t=5,0, para =α 6o. Re = 1,7 x 105
-Cp
X
-Cp
X
-Cp
X
-Cp
X
-Cp
X
79
0 0.5 1-1
0
1
2
x
0 0.5 1-1
0
1
2
x
-cp
0 0.5 1-1
0
1
2
x
-cp
0 0.5 1-1
0
1
2
x
-cp
0 0.5 1-1
0
1
2
x
-cp
t = 6,0
t =10,0
t = 9,0
t = 8,0
t = 7,0
Figura 7.8 (b) – Esteira ao longo da simulação e coeficiente de pressão médio aolongo do aerofólio NACA 0012, de t = 6,0 até t=10, para =α 6o. Re = 1,7 x 105
X
-Cp
-Cp
X
-Cp
X
-Cp
X
-Cp
X
80
Como o escoamento permanece razoavelmente aderido à superfície para o ângulo
de ataque de 6o, pode-se comparar a distribuição de pressão obtida na simulação com a
distribuição de pressão prevista por um modelo puramente potencial, obtido, neste caso,
via o Método dos Painéis com distribuição de vórtices lineares. A Figura 7.9 mostra uma
comparação da média temporal total da distribuição de pressão da simulação com o
modelo potencial. Como esperado, o pico de sucção no bordo de ataque é
significativamente reduzido, mas o comportamento global é semelhante, com uma menor
diferença de pressão ao longo do aerofólio, como esperado para um modelo viscoso.
O caráter oscilatório da esteira possui influência direta sobre o coeficiente de
sustentação, conforme pode ser visto na Figura 7.10. Neste caso, o coeficiente de
sustentação oscila em torno de um valor médio de cerca de 0,45, ligeiramente abaixo do
valor experimental de 0,56, apresentado em BLEVINS (1984) para o mesmo valor do
número de Reynolds. O coeficiente de arrasto, calculado apenas pela integração de
pressão, também é apresentado nesta mesma figura, apresentando porém, um valor de 0,05
muito maior que o obtido experimentalmente, que é da ordem de 0,01, segundo BLEVINS
(1984).
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioCp potencialContorno
Figura 7.9 – Comparação entre a distribuição de pressão média total (NACA0012, =α 6o , de t = 0 até 10), com o modelo potencial. Re = 1,7 x 105
-Cp
X
Cp Cp
81
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
Para ângulos de ataque entre 10o e 16o, a formação de estruturas de vórtices sobre a
superfície superior aumentam gradativamente, como pode ser observado pelo gráfico de
linhas de corrente para o escoamento para cada ângulo de ataque apresentadas no Apêndice
C. Quando o ângulo de 18o é atingido, o descolamento ocorre muito próximo ao bordo de
ataque, provocando uma queda súbita da sustentação, ou seja, o “estol” do aerofólio. A
Figura 7.11 apresenta a evolução da simulação para este ângulo de ataque crítico. Neste
caso, a esteira já assume uma configuração mais espessa e com desprendimento de grandes
estruturas de vorticidade, apresentando uma configuração similar à esteira de corpos
rombudos. Nesta mesma figura, pode-se observar que a distribuição de pressão apresenta
um comportamento altamente não-permanente. Inicialmente, ocorre uma diminuição da
pressão da parte superior do aerofólio. Esta região de sucção é justamente onde se forma
um grande vórtice logo após o bordo de ataque, podendo-se observar sua movimentação
desde t = 1,0 até 3,0, quando o vórtice formado chega ao centro do aerofólio e deixa sua
superfície. Logo em seguida, novos vórtices são formados periodicamente e deixam a
superfície em torno do mesmo ponto, caracterizando assim uma separação significativa do
escoamento próximo ao bordo de ataque. Este descolamento leva a uma queda brusca do
coeficiente de sustentação, entre t = 5,0 e 7,0, como pode ser visto na Figura 7.12.
Figura 7.10 – Carregamentos ao longo do tempo para o aerofólioNACA 0012, =α 6o e Re = 1,7 x 105.
Cl Cd
t
82
0 0.5 1-10123
x
-cp
0 0.5 1-10123
x
-cp
0 0.5 1-10123
x
-cp
0 0.5 1-10123
x
-cp
0 0.5 1-10123
x
-cp
t = 3,0
t = 5,0
t = 2,0
Figura 7.11 (a) – Esteira ao longo da simulação e coeficiente de pressão médio aolongo do aerofólio NACA 0012, de t = 1,0 até t=5,0, para =α 18o. Re = 1,7 x 105
-Cp
X
-Cp
X
-Cp
X
-Cp
X
X
-Cp
83
0 0.5 1-10123
x
-cp
0 0.5 1-10123
x
-cp
0 0.5 1-10123
x
-cp
0 0.5 1-10123
x
-cp
0 0.5 1-10123
x
-cp
t = 6,0
t =10,0
t = 9,0
t = 8,0
t = 7,0
Figura 7.11 (b) – Esteira ao longo da simulação e coeficiente de pressão médio aolongo do aerofólio NACA 0012, de t = 6,0 até t=10, para =α 18o. Re = 1,7 x 105
X
-Cp
X
-Cp
X
-Cp
X
-Cp
X
-Cp
84
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
Para ângulos de ataque ainda maiores, o escoamento torna-se altamente instável,
apresentando grandes flutuações com o tempo nos valores dos coeficientes de arrasto e
sustentação e uma esteira cada vez mais espessa e de grande riqueza de escalas. Os
resultados para os ângulos de 20o, 25o, 30o, 45o, 60o e 75o, estão mostrados no Apêndice C.
O resultado para 90o está apresentado separadamente no Apêndice D, pois devido à
simetria inicial do escoamento, foi necessária uma simulação mais longa para ocorrer o
desprendimento alternado de vorticidade. Enquanto que o coeficiente de sustentação oscila
em torno de zero, o coeficiente de arrasto médio temporal, a partir do início do
desprendimento alternado de vorticidade (t = 30) cresce significativamente. Este resultado
não é apresentado com os anteriores, pois utiliza parâmetros numéricos diferentes.
Com a finalidade de comparar os resultados obtidos com os resultados
experimentais, foi calculado, para cada ângulo de ataque, os coeficientes de arrasto e
sustentação médios entre t = 5,0 e t = 10,0, onde supôs-se que o transiente numérico e as
acomodações iniciais do escoamento já teriam terminado. Os resultados para Reynolds
1,7 x 105 para toda a faixa de simulação estão apresentados na Figura 7.13. A Figura 7.13a
está apresentado o coeficiente de sustentação e a Figura 7.13b, o coeficiente de arrasto.
Figura 7.12 – Carregamentos ao longo do tempo para o aerofólioNACA 0012, =α 18o e Re = 1,7 x 105.
Cl Cd
t
85
Juntamente com os resultados numéricos, estão apresentados os resultados experimentais
de CRITZOS et al. (1954), obtidos de LEWIS (1991) e HOERNER & BORST (1985).
Nestas figuras é notável a concordância qualitativa de ambos os resultados com os
experimentos.
α
dC
Figura 7.13 (b) – Comparação com resultados experimentais do coeficientede arrasto em função do ângulo de ataque para o aerofólio NACA 0012.
0
0,5
1
1,5
2
2,5
3
3,5
4
0 10 20 30 40 50 60 70 80 90 100
CRITZOS et al. (1954) - Re=1,8e06
Resultados - Re=1,7e05
dC
α
0
0,5
1
1,5
2
2,5
0 10 20 30 40 50 60 70 80 90
CRITZOS et al. (1954) - Re=1,8e06Resultados - Re=1,7e05
Figura 7.13 (a) – Comparação com resultados experimentais do coeficientede sustentação em função do ângulo de ataque para o aerofólio NACA 0012.
α
lC
86
Conforme já mencionado, existe um particular interesse na região de “estol” do aerofólio e,
por este motivo, são apresentados nas Figuras 7.14a e 7.14b os coeficientes de sustentação
e arrasto respectivamente, para a faixa de 0 a 20o, comparados agora com outros resultados
experimentais e inclusive com os resultados numéricos de AKBARI & PRICE (2003) para
Re = 104, que também empregam o Método dos Vórtices Discretos. Os presentes
resultados mostram uma melhora significativa no que tange ao coeficiente de sustentação.
Na Figura 7.14a os coeficientes de sustentação experimentais são apresentados para Re =
1,7 x 105 (BLEVINS, 1984) e Re = 3,0 x 106 (ABBOT & VON DOEHOFF, 1951). O valor
encontrado para o ângulo de sustentação máxima foi 16o com coeficiente de sustentação
correspondente de 1,15. A forte dependência do coeficiente de sustentação máximo com o
número de Reynolds não foi capturado nas simulações, inclusive o ângulo de “estol”
encontrado foi mais próximo ao valor correspondente a Reynolds 3,0 x 106 do que do
próprio valor correspondente a 1,7 x 105. Este fato se deve provavelmente às deficiências
do modelo de difusão de vorticidade pelo Método do Avanço Randômico, que se aproxima
da equação de difusão quanto maior for o número de Reynolds. Algumas simulações para
Reynolds de 1,7 x 106 estão listadas no Apêndice E, porém não demonstraram grandes
variações em relação aos demais resultados.
Em uma comparação quantitativa dos valores de arrasto e sustentação observa-se
que os valores de sustentação tendem a ser mais baixos que os experimentais, enquanto
que os valores de arrasto muito mais altos. O valor excessivamente alto do coeficiente de
arrasto e os baixos valores de sustentação originam-se principalmente da discretização do
campo de vorticidade, que na realidade é contínuo. Esta discretização provoca
instabilidades na esteira que, dependendo de sua escala, o modelo não é capaz de dissipá-
las por ação da viscosidade. Sendo assim, a esteira tende a ser mais espessa que o
observado experimentalmente, ocasionando um alto arrasto e reduzindo a sustentação. Os
efeitos turbulentos presentes no escoamento da esteira do corpo, que afetam o escoamento
da camada limite devem ser modelados em suas pequenas escalas ou caputurados na
simulação através do aumento do número de vórtices presentes na esteira e na própria
camada limite.
Após a região de “estol”, os carregamentos possuem um comportamento muito
instável e oscilam fortemente, sem apresentar um valor médio bem definido. Os valores
apresentados na Figura 7.13a foram calculados da mesma forma que os demais da região
87
antes do “estol”, porém seria adequado simular para tempos mais longos, buscando um
comportamento periódico mais bem definido para os carregamentos. Todavia, a simulação
de tempos ainda mais longos torna-se proibitiva do ponto de vista de custo computacional.
Um teste realizado para =α 20o é apresentado no Apêndice F, mostrando que o coeficiente
de sustentação tem uma tendência de decrescer ainda mais para tempos mais longos,
ficando seu valor médio mais próximo dos experimentais.
Figura 7.14 (a) – Comparação com resultados experimentais do coeficientede sustentação na região de “estol” para o aerofólio NACA 0012.
α
Figura 7.14 (b) – Comparação com resultados experimentais do coeficientede arrasto na região de “estol” para o aerofólio NACA 0012.
α
dC
lC
0,0
0,2
0,4
0,6
0,8
1,0
1,2
1,4
1,6
1,8
2,0
0 5 10 15 20
BLEVINS (1984) - Re=1,7e05ABBOT & VON DOEHOFF (1951) - Re=3,0e06AKBARI & PRICE (2003) - Re = 1,0e04Resultados - Re=1,7e05
0,00
0,05
0,10
0,15
0,20
0,25
0,30
0,35
0,40
0 2 4 6 8 10 12 14 16 18
BLEVINS (1984) - Re=1,7e05ABBOT & VON DOEHOFF (1951) - Re=3,0e06AKBARI & PRICE (2003) - Re = 1,0e04Resultados - Re=1,7e05
88
7.6 Aerofólio NACA 4415
Utilizando novamente os parâmetros descritos na seção 7.2, o escoamento ao redor
do aerofólio NACA 4415 foi simulado para diversos ângulos de ataque e Re = 1,7 x 105.
Este aerofólio foi escolhido devido à sua geometria com curvatura e sua similaridade com
aerofólios utilizados atualmente em aplicações comerciais. Este tipo de aerofólio foi
utilizado nos primeiros projetos de turbinas eólicas (MANWELL et al., 2001). Sua
curvatura lhe confere um “estol” mais suave, mesmo para altos valores do número de
Reynolds, pois a curvatura no bordo de ataque faz com que o escoamento, mesmo com
ângulo de ataque maior que zero, encontre uma geometria um pouco mais alinhada com
sua direção na parte dianteira do aerofólio.
Como esperado a geometria do aerofólio provoca uma assimetria do escoamento
mesmo com ângulo de ataque nulo, gerando uma força de sustentação positiva, conforme
pode ser visto na Figura 7.15. Os resultados gráficos para a posição dos vórtices na esteira
no último instante de tempo, para as linhas de corrente, distribuição de pressão e
carregamentos estão listados no Apêndice G para cada ângulo de ataque. No Apêndice G
estão todos os resultados, inclusive para ângulos após o “estol”, que não são apresentados
para comparação nesta seção devido à escassez de resultados experimentais disponíveis na
literatura. Além disso, não foram realizadas simulações com ângulos de ataque maiores
que 20o para este aerofólio, devido ao interesse de se determinar a região de “estol” apenas.
Um resumo dos resultados obtidos que compara os coeficientes de arrasto e
sustentação médios das simulações com os resultados experimentais é apresentado nas
Figuras 7.16a e b. No âmbito de avaliação global do modelo, as informações obtidas por
estes resultados são bem similares ao aerofólio anterior. O ângulo de “estol” é previsto
com boa precisão, o coeficiente de sustentação situa-se ligeiramente abaixo dos valores
Figura 7.15– Esteira ao final da simulação para um aerofólio NACA 4415 comângulo de ataque nulo e Re = 1,7 x 105.
89
experimentais de um valor de aproximadamente 0,2 , enquanto que o arrasto é muito mais
alto que o obtido nos experimentos.
Figura 7.16 (b) – Comparação com resultados experimentais do coeficientede arrasto na região de “estol” para o aerofólio NACA 4415.
α
dC
0,0
0,1
0,1
0,2
0,2
0,3
0,3
0,4
0 2 4 6 8 10 12 14 16 18
ABBOT & VON DOEHOFF (1951) - Re=3,0e06
Resultados (Re = 1,7e05)
Figura 7.16 (a) – Comparação com resultados experimentais do coeficientede arrasto na região de “estol” para o aerofólio NACA 4415.
α
lC
0,0
0,2
0,4
0,6
0,8
1,0
1,2
1,4
1,6
0 2 4 6 8 10 12 14 16 18
ABBOT & VON DOEHOFF (1951) - Re=3,0e06
Resultados (Re = 1,7e05)
90
7.7 Outros Aerofólios
O programa desenvolvido pode ser aplicado facilmente para outros tipos de
aerofólios, sobretudo da família NACA de 4 dígitos, para os quais rotinas já estão
implementadas e geram a geometria desejada sem necessidade de alterações no programa.
A família NACA de 5 dígitos, outros aerofólios e outras geometrias também podem
ser simulados com pequenas alterações no programa FORTRAN desenvolvido. Porém,
cuidados especiais devem ser tomados no caso de geometrias com superfícies muito finas,
onde os pontos de controle dos painéis no bordo de fuga, principalmente, encontram-se
muito próximos uns dos outros. O perfil NACA 63-415, por exemplo, possui o bordo de
fuga extremamente fino, conforme a Figura 7.17 ilustra. Por este motivo, a simulação deste
aerofólio apresenta erros numéricos no cálculo do campo de velocidade nesta região.
Nestes casos, são necessários tratamentos numéricos especiais, tais como a correção da
diagonal oposta do sistema ou um grande refinamento dos painéis nesta região.
Figura 7.17 – Detalhe do bordo de fuga do aerofólio NACA 63-415
91
CAPÍTULO 8: CONCLUSÕES
O presente trabalho apresenta um algoritmo numérico que emprega o Método de
Vórtices Discretos associado ao Método dos Painéis com correções para a simulação
numérica do escoamento incompressível, bidimensional, não permanente e de alto número
de Reynolds ao redor de aerofólios da família NACA. Dois aerofólios são estudados em
detalhe: NACA 0012 e NACA 4415. A vorticidade do escoamento é modelada como uma
nuvem de vórtices discretos, cujo movimento é calculado através da abordagem
lagrangiana a cada passo de tempo, de acordo com os processos de convecção e difusão
locais de vorticidade.
Na implementação do Método de Vórtices Discretos, utiliza-se o Método do
Avanço Randômico, para simular a difusão de vorticidade do escoamento, e o esquema de
Adams-Bashforth de segunda ordem, para o avanço temporal convectivo dos vórtices na
nuvem. As distribuições de vorticidade na superfície do corpo são calculadas de tal forma a
satisfazer à condição de impenetrabilidade e à conservação de circulação. A cada passo de
tempo são liberados novos vórtices na vizinhança da superfície do aerofólio, que
deslocam-se por ação dos efeitos de difusão e convecção, configurando a camada limite do
corpo e a esteira que se forma à jusante.
O Método de Vórtices Discretos associado ao Método dos Painéis com correções
mostrou-se capaz de descrever escoamentos qualitativamente complexos com alto número
de Reynolds ao redor de aerofólios. Conforme os resultados obtidos, foi possível
identificar características de comportamento do escoamento simulado semelhantes às
observações experimentais e demais resultados numéricos disponíveis na literatura
científica.
O comportamento particular do escoamento ao redor de aerofólios em função do
ângulo de ataque, principalmente quanto à evolução do coeficiente de sustentação, foi
capturado com sucesso nas simulações. Os resultados numéricos obtidos, em comparação
com os resultados de AKBARI & PRICE (2003) para o perfil NACA 0012, mostra uma
melhora significativa no cálculo do coeficiente de sustentação em função do ângulo de
ataque, quando comparados a resultados experimentais. Embora o número de Reynolds
92
utilizado neste trabalho tenha sido da ordem de 105 e os resultados de AKBARI & PRICE
(2003) sejam para o número de Reynolds de 104, a melhora é representativa principalmente
na região de baixo ângulo de ataque, onde o coeficiente de sustentação praticamente
independe do valor do número de Reynolds. O ângulo de sustentação máxima, onde se
inicia o “estol” sobre o aerofólio, foi previsto com boa concordância para os aerofólios
escolhidos.
A qualidade dos resultados comprova a superioridade do modelo de singularidades
tipo vórtice de distribuição linear em comparação coma distribuição constante, para o
Método dos Painéis, descrevendo satisfatoriamente o campo de velocidade ao redor do
aerofólio. Esta modelagem implica em uma formulação que inclui uma equação de
conservação de circulação e não gera um sistema de equações algébricas lineares sobre-
determinado. O sistema gerado pode, portanto, ser resolvido com maior rapidez e menores
erros numéricos.
As simulações para o aerofólio NACA 0012 mostram um grande conjunto de
resultados que cobrem toda a faixa de ângulos de ataque, desde 0 até 90o, mostrando desde
os escoamentos aderidos à superfície, com baixo ângulo de ataque, passando pela região de
“estol” e seguindo para complexos escoamentos completamente descolados, que
apresentam a formação de grandes estruturas de vorticidade.
Os resultados para o coeficiente de arrasto fornecem valores superestimados em
relação aos dados experimentais utilizados para comparação. As discrepâncias observadas
podem ser atribuídas a discretização e dinâmica da vorticidade na camada limite, que
necessita de um número maior de vórtices para melhor representação. Conforme observado
nos resultados obtidos, mesmo com o número elevado de vórtices utilizado, existe uma
tendência à separação prematura da camada limite sobre os aerofólios.
As simulações realizadas indicam que o modelo numérico desenvolvido possui
resultados tão melhores, quão maior for o número de vórtices gerados a cada passo de
tempo. A distância de geração dos vórtices possui grande influência na simulação. Os
resultados obtidos indicam que os coeficientes de arrasto e sustentação se aproximam dos
valores experimentais quando esta distância decresce, porém existe um limite inferior
numérico quando a simulação torna-se instável e desprovida de significado físico. O
93
incremento de tempo deve ser mantido o mais baixo possível, porém a redução do ∆t a
valores abaixo do utilizado (0,025) requer tempos de simulação muito grandes para o
mesmo tempo adimensional total final.
Os gráficos das linhas de corrente, apresentado em regiões próximas do aerofólio,
corroboram a conclusão de que a Física do escoamento está sendo adequadamente
representada pelo modelo numérico.
As simulações para o aerofólio NACA 4415 foram realizadas apenas até se
observar a ocorrência do “estol”, descrevendo o escoamento para esta geometria
essencialmente com a mesma precisão observada para o aerofólio NACA 0012. Este
aerofólio possui curvatura não nula e maior espessura que o anterior, e, por este motivo, é
de análise mais limitada quando utilizados os métodos puramente potenciais.
Conforme observado neste trabalho, e comentado por GUEDES (2003) e MUSTTO
(2004), as limitações da utilização do Método dos Vórtices Discretos são: o grande esforço
computacional envolvido, quando se calcula as velocidades induzidas pelos vórtices na
nuvem através da Lei de Biot-Savart; a dificuldade, mas não impossibilidade, de aplicação
para escoamento tridimensionais; as complicações práticas para modelagem de
escoamentos internos, pois o escoamento requer uma modelagem da vorticidade em todo o
domínio fluido, o que implica em um número de vórtices muito grande no escoamento; e,
principalmente, a grande dependência da qualidade das simulações com os parâmetros
numéricos escolhidos.
Por outro lado, o Método possui grandes vantagens, podendo-se destacar: não é
necessário dividir a região fluida em uma malha computacional, o que elimina todos os
problemas inerentes à geração de malha; a satisfação das condições de contorno com maior
precisão em regiões distantes do corpo; e a simulação do escoamento somente nas regiões
rotacionais, onde a viscosidade efetivamente atua.
Como sugestão para melhorias futuras no modelo, pode-se incluir uma modelagem
de turbulência específica, a exemplo do modelo utilizado por MUSTTO (2004), caso o
objetivo seja simular aerofólios com valores do número de Reynolds maior que 106, onde
geralmente ocorre transição da camada limite. A utilização do esquema de multipolos para
94
acelerar o algoritmo na etapa convectiva, onde ocorre o cálculo das velocidades induzidas
pela nuvem de vórtices, possibilita simular tempos maiores com um número de vórtices
maior no escoamento.
Outros aspectos da modelagem também oferecem oportunidade para melhoria, as
singularidades distribuídas sobre os painéis por exemplo, podem incluir fontes, de forma
que as condições de impenetrabilidade e aderência possam ser explicitamente empregadas.
Além disso, formulações mais recentes para o cálculo dos carregamentos ao redor do
corpo, tal como em SHINTANI & AKAMATSU (1994), podem reduzir significativamente
as oscilações numéricas dos coeficientes de arrasto e sustentação, de forma que as médias
obtidas podem ser mais representativas.
Um estudo mais detalhado dos parâmetros numéricos é necessário, sobretudo da
distância de geração de vórtices e do raio do núcleo, tendo em vista a grande dependência
que as simulações apresentam desses parâmetros. A substituição do modelo de Avanço
Randômico para a difusão de vorticidade por um modelo determinístico e adequado
também para Reynolds mais baixos, tal como o Método de Redistribuição de Vorticidade,
pode melhorar significativamente a precisão do cálculo do coeficiente de sustentação
máximo para diferentes números de Reynolds.
95
REFERÊNCIAS
ABBOTT, I. H. & VON DOENHOFF, A., 1959, Theory of wing sections, Dover
Publications, New York.
AKBARI, M. H. & PRICE, S. J., 2003, “Simulation of Dynamic Stall For a NACA
0012 Airfoil Using a Vortex Method”, Journal of Fluids and Structures, v. 17, pp. 855-874.
ALCÂNTARA PEREIRA, L. A., 1999, “Simulação Numérica do Escoamento em
Torno de um Corpo de Forma Arbitrária Utilizando o Método de Vórtice Discretos”, Tese
de M.Sc., EFEI, Itajubá, MG, Brasil.
ALCÂNTARA PEREIRA, L. A., 2002, “Simulação Numérica do Escoamento ao
Redor de Perfis Aerodinâmicos Montados em Grades Lineares de Turbomáquinas
Utilizando o Método de Vórtices Discretos com Modelagem Sub-Malha de Turbulência”,
Tese de D.Sc., EFEI, Itajubá, MG, Brasil.
ALCÂNTARA PEREIRA, L. A., HIRATA, M. H. & MAZANARES, N., 2004, “Wake
and Aerodynamics Loads in Multiple Bodies—Application to Turbomachinery Blade
Rows”, Journal of Wind Engineering and Industrial Aerodynamics, v. 92, pp. 477–491.
ALCÂNTARA PEREIRA, L. A., RICCI, J. E. R., HIRATA, M. H. & SILVEIRA
NETO, A., 2002, “Simulation of Vortex-Shedding Flow About a Circular Cylinder with
Turbulence Modelling”, Computational Fluid Dynamics, v. 11, pp. 315-322.
ANDERSON, J. D., 1991, Fundamentals of Aerodynamics, McGraw-Hill, 2nd edition.
AREF, H., 1983, “Integrable, Chaotic and Turbulent Vortex Motion in two-
dimensional Flows”, Annual Review of Fluid Mechanics, v. 15, pp. 345-389.
AREF, H. & KAMBE, T., 1988, “Report on the IUTAM symposium: Fundamental
Aspects of Vortex Motion”, Journal of Fluid Mechanics, v. 190, pp. 571-595.
96
BARNES, J. E. & HUT, P., 1986, “A Hierarchical O(N logN) Force-Calculation
Algorithm”, Nature, v. 324, pp. 446-449.
BATCHELOR, G. K., 1967, An Introduction to Fluid Dynamics, Cambridge
University Press.
BERTAGNOLIO, F., SORENSEN, N., JOHANSEN, J., FUGLSANG, P., 2001, Wind
Turbine Airfoil Catalog, Risoe National Laboratory, Risoe - R - 1280 (EN), Roskilde,
Dinamarca.
BERTIN, J. J. & SMITH, M. L., 1998, Aerodynamics for Engineers, Prentice Hall, 3rd
edition.
BLEVINS, R. D., 1984, Applied Fluid Dynamics Handbook, Van Nostrand Reinhold,
Co.
BREBBIA, C. A., TELLES, J. C. F., WROBEL, L. C., 1984, "Boundary Element
Techniques: Theory and Applications in Engineering", Springer-Verlag, Berlin-
Heidelberg.
CARREIRO, F. R., 2002, “Simulação do Escoamento Bidimensional, Incompressível,
de Alto Reynolds, em Torno de Corpos via Transformação Conforme e Método de
Vórtices”, Tese de M.Sc., COPPE/UFRJ, Rio de Janeiro, RJ, Brasil.
CARRIER, J., GREENGARD, L. & ROKHLIN, V., 1988, “Multipole Algorithm for
Particle Simulation”, SIAM Journal of Scientific Statistics and Computation, v. 9, n. 4, pp.
669-686.
CHORIN, A. J., 1973, “Numerical Study of Slightly Viscous Flows”, Jounal of Fluid
Mechanics, v. 57, pp. 785-796.
CHORIN, A. J. & BERNARD, P. S., 1973, “Discretization of a Vortes Sheet with an
Example of Roll up”, Journal of Computational Physics, v. 13, pp.423-429.
97
CHURCHILL, R. V. & BROWN, J. W., 1984, Complex Variable Theory, McGraw-
Hill, NY.
CLARKE, N. R. & TUTTY, 1993, O. R., “Construction and Validation of a Discrete
Vortex Method for the Two-dimensional Incompressible Navier-Stokes Equations”,
Computers Fluids, v. 23, pp. 751-783.
COTTET, G. H., KOUMOUTSAKOS, P. D., 2000, Vortex Methods: Theory and
Practice, Cambridge University Press, Cambridge, UK.
COTTET, G. H., KOUMOUTSAKOS, P. & OULD SALIHI, M. L., 2000, “Vortex
Methods with Spatially Varying Cores”, Journal of Computational Physics, v. 162, pp.
164–185.
CRITZOS, C. C., HEYSON, H. H., BOSWINKLE, R. W., 1954, “Aerodynamic
Characteristics of NACA 0012 Airfoil Section”, NACA Tech. 3361.
DAHL, K.S., FUGLSANG, P., ANTONIOU, I., 1999, “Experimental Verifications of
the New RISO-A1 Airfoil Family forWind Turbines”, Proc. of EWEC’99, pp.85-88.
DEWI - Deutsches Windenergie Institut (Instituto Alemão de Energia Eólica), Wind
Energy Report, 1998.
FONSECA, G. F., 1996, “Interação Bidimensional entre um Aerofólio e um Vórtice
Próximo a uma Superfície: Influência da Esteira”, Tese de M.Sc., COPPE/UFRJ, Rio de
Janeiro, RJ, Brasil.
GREENGARD, L., 1985, “The Core Spreading Vortex Method Approximates the
Wrong Equation”, Journal of Computational Physics, v.61, pp.345-348.
GREENGARD, L. & ROKHLIN, V., 1987, “A Fast Algorithm for Particle
Simulations”, Journal of Computational Physics, v. 73, pp. 325-348.
98
GUEDES, V. G., 2003, “Estudo Numérico do Escoamento ao Redor de Cilindros
Circulares e Retangulares Utilizando o Método de Vórtices”, Tese de D.Sc.,
COPPE/UFRJ, Rio de Janeiro, RJ, Brasil.
HESS, J. L., 1990, “Panel Methods in Computational Fluid Dynamics,” Annual
Review of Fluid Mechanics, Vol. 22, pp. 255-274.
HESS, J. L. & SMITH, A. M. O., 1966, “Calculation of Potencial Flow About
Arbitrary Bodies”, Progress in Aeronautical Sciences, Vol. 3, Pergamon Press, New York,
138p.
HIRATA, M. H., 2000, “O Método de Vórtices com Modelagem da Turbulência”,
Palestra Oficial, proferida no Congresso Norte-Nordeste de Engenharia Mecânica –
CONEN 2000, Natal, RN, Brasil.
HOERNER, S. F. & BORST, H. V., 1985, Fluid-Dynamic Lift: Practical Information
on Aerodynamic and Hydrodynamic Lift, Hoerner Fluid Dynamics.
KAMEMOTO, K., 1994, “Development of Vortex Methods for Grid-Free Lagrangian
Direct Numerical Simulation,” Proceedings of the 3th. JSME-KSME Fluids Engeneering
Conference, July 25-27, Sendai, Japan.
KAMEMOTO, K., 2004, “On Contribution of Advanced Vortex Element Methods
Toward Virtual Reality of Unsteady Vortical Flows in the New Generation of CFD”,
Invited Lecture - Proceedings of the 10th Brazilian Congress of Thermal Sciences and
Engineering, ENCIT, Rio de Janeiro, RJ, Brasil.
KARAMCHETTI, K., 1980, Principles of Ideal-Fluid Dynamics, R. E. Krieger
Publishing Co.
KATZ, J. & PLOTKIN, A., 2001, Low - Speed Aerodynamics, 2nd edition, Cambridge
University Press.
99
KOUMOUTSAKOS, P. & LEONARD, A., 1995, “High Resolution Simulations of the
Flow Around an Impulsively Started Cylinder Using Vortex Methods”, Journal of Fluid
Mechanics, v. 296, pp. 1-38.
KREYSZIG, E., 1999, Advanced Engineering Mathematics, 8th edition, J. Wiley.
KUNDU, P. K. & COHEN, M. C., 2002, Fluid Mechanics, 2nd edition, Academic
Press, San Diego.
KUWAHARA, K. & TAKAMI, H., 1973, “Numerical Studies of two-dimensional
Vortex Motion by a System of Point Vortices”, Journal of Physics Society of Japan, v. 34,
pp. 247-253.
LEONARD, A., 1980, “Vortex Method for Flow Simulation”, Journal of
Computational Physics, v. 37, pp. 289-335.
LEONARD, A., 1985, “Computing three-dimensional Inconpressible Flows with
Vortex Elements”, Annual Review of Fluid Mechanics, v. 17, pp. 523-559.
LEWIS, R. I., 1981, “Surface Vorticity Modeling of Separated Flows for Two-
Dimensional Bluff Bodies of Arbitrary Shape”, Journal of Mechanical Engineering and
Sciences, I. Mech. E., Vol. 23, No. 1.
LEWIS, R. I., 1991, Vortex Element Methods for Fluid Dynamic Analysis of
Engineering Systems, Cambridge, Cambridge University Press.
LEWIS, R. I., 1999, “Vortex Element Methods, the Most Natural Approach to Flow
Simulation – A Review of Methodology with Applications”, Proceedings of 1st International
Conference on Vortex Methods, Kobe, Japão, pp. 1–15.
MALISKA, C. R., 2004, Transferência de Calor e Mecânica dos Fluidos
Computacional, Rio de Janeiro, LTC.
100
MALTA, D., 1998, “Aplicação do Método de Vórtices ao Escoamento de Alto Número
de Reynolds em torno de um Cilindro Circular”, Tese de M.Sc., COPPE/UFRJ, Rio de
Janeiro, RJ, Brasil.
MANWELL, J.F., McGOWAN J.G., ROGERS, A.L., “Wind Energy Explained,
Theory, Design and Application”, John Wiley & Sons, West Sussex - England, 2001.
MARTENSEN, E., 1959, “The Calculation of the Pressure Distribution on a Cascade
of thick Airfoils by means of Fredholm Integral Equations of the Second Kind”,
Aerodynamics Experimental Station, Göttingen.
McCROSKEY, W. J., 1982, “Unsteady Airfoils”, Annual Review of Fluid Mechanics,
v. 14, pp. 285-311.
MILNE-THOMSON, L. M., 1955, Theoretical Hydrodynamics, MacMillan & Co, NY.
MILNE-THOMSON, L. M.,1958, Theoretical Aerodynamics, Dover Publications, NY.
MOOK, D. T. & DONG, B., 1994, “Perspective: Numerical Simulations of Wakes and
Blade-Vortex Interaction”, Journal of Fluids Engineering, Vol. 116, pp. 6-21.
MORAN, J., 1984, An Introduction to Theoretical and Computational Aerodynamics,
John Wiley & Sons, New York.
MUSTTO, A. A., 1998, “Simulação Numérica do Escoamento em Torno de um
Cilindro Circular com e sem Rotação Utilizando o Método de Vórtices Discretos”, Tese de
M.Sc., COPPE/UFRJ, Rio de Janeiro, RJ, Brasil.
MUSTTO, A. A., 2004, “Simulação Numérica do Escoamento Turbulento em Torno de
um Cilindro Circular via Método de Vórtices”, Tese de D.Sc., COPPE/UFRJ, Rio de
Janeiro, RJ, Brasil.
101
OGAMI, Y. & AKAMATSU, T., 1991, “Viscous Flow Simulation Using the Discrete
Vortex Model - The Diffusion Velocity Method”, Computers and Fluids, v. 19, pp. 433-441.
OGAMI, Y. & AYANO, Y., 1995, “Flows around a Circular Cylinder Simulated by
the Viscous Vortex Model: the Diffusion Velocity Method”, Computational Fluid
Dynamics Journal, v. 4, n. 3, (Oct), pp. 383-399.
OJIMA, A. & KAMEMOTO, K., 1999, “Numerical Simulation of UnsteadyFlow
Around a Sphere by a Vortex Method for Re Number from 300 to 1000”, Proceedings of 1st
International Conference on Vortex Methods, Kobe, Japão.
OJIMA, A. & KAMEMOTO, K., 2001, “Numerical Simulation of Unsteady Flow
Through a Horizontal Axis Wind Turbine”, Proceedings of the 2nd International Conference
on Vortex Methods, ITU Maçka Campus, Istanbul, Turquia, pp. 173-180.
PATANKAR, S. V., 1980, Numerical Heat Transfer and Fluid Flow, Nova Iorque,
McGraw-Hill.
PEREIRA, L. H. G. & BODSTEIN, G. C. R., 2004, “Método dos Painéis com
Distribuições de Singularidade Quadráticas Aplicados a Escoamentos Bidimensionais
sobre Aerofólios”, ENCIT – 10 Congresso Brasileiro de Ciências Térmicas, Rio de Janeiro,
RJ, Brasil.
PEREIRA, L. H. G., SILVA, D. F. C. & BODSTEIN, G. C. R., 2004, “Estudo
Numérico do Escoamento Potencial em Torno de um Aerofólio Utilizando o Método dos
Painéis”, ENCIT – 10 Congresso Brasileiro de Ciências Térmicas, Rio de Janeiro, RJ,
Brasil.
PLOUHAMS, P., WICKELMANS, G. S., SALMON, J. K., LEONARD, A. &
WARREN, M. S., 2002, “Vortex Methods for Direct Numerical Simulation of Three-
Dimensional Bluff Body Flows: Application to the Sphere at Re = 300, 500, and 1000”,
Journal of Computational Physics, v. 178, pp. 427-463.
102
PORTHOUSE, D. T. C. & LEWIS, R. I., 1981, “Simulation of Viscous Diffusion for
Extension of the Surface Vorticity Method to Boundary and Separated Flows”, Journal of
Mechanical Engineering and Sciences, I. Mech. E., Vol. 23, No. 3, pp. 157-167.
REPOWER SYSTEMS, Notícia sobre a construção de uma turbina eólica de 5MW de
potência, disponível em http://www.repower5m.com/index_flash.htm, Acesso em 20 de
jun. de 2004.
RICCI, J. E. R., 2002, “Simulação Numérica do Escoamento ao Redor de um Corpo de
Forma Arbitrária, Estacionado nas Imediações de uma Superfície Plana, com Emprego de
Método de Vórtices”, Tese de D.Sc., UNIFEI, Itajubá, MG, Brasil.
ROSENHEAD, L., 1931, “The Formation of Vortices from a surface of Discontinuity”,
Proceedings of Royal Society A, v. 134, pp. 170-192.
ROSENHEAD, L., 1963, Laminar Boundary Layers, Dover Publications, New York.
RUGGIERO, M. A. G., LOPES, V. L. R., 1997, Cálculo Numérico – Aspectos
Teóricos e Computacionais, São Paulo, Editora Makron Books.
SAFFMAN, P. G., BAKER, G. R., 1979, “Vortex Interactions”, Annual Review on
Fluid Mechanics, v. 11, pp. 95-112.
SARPKAYA, T., 1989, “Computational Methods with Vortices – The 1988 Freeman
Scholar Lecture”, Journal of Fluids Engineering, v.111, pp.5-52.
SARPKAYA, T., 1994, “Vortex Elment Methods for Flow Simulation”, Advances in
Applied Mechanics, v. 31, pp. 113-247.
SCHLINCHTING, H. & GERSTEN, K., 2000, Boundary Layer Theory, Springer-
Verlag, Berlim.
SHANKAR, S. & VAN DOMMELEN, L., 1996, “A New Diffusion Procedure for
Vortex Methods”, Journal of Computational Physics, v. 127, pp. 99-164.
103
SHINTANI, M & AKAMATSU, Y., 1994, "Investigation of Two-Dimensional
Discrete Vortex Method with Viscous Diffusion Model", Compuational Fluid Dynamics
Journal, vo. 3, no. 2, pp. 237-254.
SILVA, D. F. C., PEREIRA, L. H. G., BODSTEIN, G. C. R., 2004, Aerodynamic
Modeling of the Blade Geometry of Horizontal Axis Wind Turbines Using a Panel
Method”, World Wind Energy Conference, Pequim, China.
SILVEIRA NETO, A., 2004, Turbulência nos Fluidos Aplicada, apostila do curso
ministrado no programa de verão do LNCC, Petrópolis.
SPARLAT, P. R., LEONARD, A., 1981, “Computation of Separated Flows by a
Vortex Tracing Algorithm”, AIAA 14th Fluid and Plasma Dynamics Conference, Palo
Alto, Califórnia, AIAA-81-1246.
SPREITER, J. R. & SACKS, A. H., 1951, “The Rolling up of the Trailing Vortex and
its Effect on the Downwash behind Wings”, Journal of Aeronautical Sciences, v. 18, pp.
21-32.
TAKEDA, K. & TUTTY, O. R., 1997, “A Comparison of Four Viscous Models for the
Discrete Vortex Method”, 13rd AIAA Compuational Fluid Dynamics Conference,
Snowmass, CO, 29th June – 2nd July 1997, section No CFD-16, Paper No. 97-1977.
THWAITES, B., 1960, Incompressible Aerodynamics, Dover Publications, NY.
VAN DYKE, M., 1975, Perturbation Methods in Fluid Mechanics, The Parabolic
Press, Stanford.
VEZZA, M. & GALBRAITH, R. A. McD., 1985, “A Inviscid Model of Unsteady
Aerofoil Flow with Fixed Upper Surface Separation”, International Journal for Numerical
Methods in Fluids, Vol. 5, pp. 577-592.
104
XU, C., 1998, “Surface Vorticity Modeling of Flow Around Airfoils”, Computational
Mechanics, v. 21, pp.526-532.
XU, C., 1999, “A Vortex Method for Separated Flow Around an Airfoil with a
Detached Spoiler”, Computational Mechanics, v. 23, pp.271-278.
105
APÊNDICE A
SINGULARIDADES DE ORDEM SUPERIOR
NO MÉTODO DOS PAINÉIS
Como enfatizado, a modelagem desenvolvida neste trabalho utiliza o Método de
Painéis com distribuição de vorticidade linear sobre os painéis acoplado ao Método de
Vórtices para a simulação do escoamento ao redor de um aerofólio. Entretanto, vale
ressaltar que singularidades de ordem superior também podem ser utilizadas.
Para complementar a análise relativa ao Método dos Painéis, a Tabela A.1 apresenta o
potencial de velocidade e as componentes de velocidades nas direções x e z, no sistema de
coordenadas do painel, para diversos tipos de singularidades que podem ser empregados na
análise do escoamento em questão, da distribuição constante até a distribuição quadrática.
Estas expressões utilizam algumas variáveis auxiliares r1, r2, 1 2,θ θ , apresentadas na Figura
3.4b e dadas pelas expressões (3.10a,b) e (3.11a,b). Como pode ser claramente visto, a
utilização de distribuições de ordem superior, como a quadrática, aumenta em muito a
complexidade da formulação. Modelos mais elaborados como esse não necessariamente
correspondem ao melhor custo benefício para a simulação numérica do escoamento.
As vantagens e desvantagens da aplicação de cada um destes modelos pode ser vista
em PEREIRA et al. , (2004) e PEREIRA & BODSTEIN, (2004).
106
Ele
men
to
Constante
( ) 0xλ λ=
Linear
( ) 0 1x xλ λ λ= +
Quadrática
( ) 20 1 2x x xλ λ λ λ= + +
Font
e
( ) ( )
( )
01 1 2 2
2 1
[ ln ln2
]
x x r x x r
z
σφπ
θ θ
= − − − +
+ −
0 2
1
ln2
rur
σπ
= −
( )02 12
w σ θ θπ
= −
( ) ( )( ) ( )
2 2 2 2 2 211 1 2 2
2 1 2 1
[ ln ln4
2 ]
x x z r x x z r
xz x x x
σφπ
θ θ
= − − − − − +
+ − − −
( ) ( )1 22 1 2 1
1
ln2
ru z x x xr
σ θ θπ
= − − − −
( )1 22 1
1
ln2
rw z xr
σ θ θπ
= + −
( )( ) ( )( )( )
( )
2 2 2 221 2 1 2
2 2 3 22 1 2 2
23 2 2 2 21 1 2
1
[212
+2 3 log
log 9 log ]
x z x x x x x
z x z x r
rx r x x zr
σφπ
θ θ
= − − + − +
− − + +
− − −
( ) ( )
( )
2 221 2 1 2 2 1
22 2 2
21
[2 44
log ]
u x x x x x xz
rx zr
σ θ θπ
= − + − + − +
− −
( )( ) ( )2
2 22 22 1 2 1 2
1
[ log ]2
rw x z z x x xz
rσ
θ θπ
= − − + − +
Dip
olo
( )02 12
µφ θ θπ
= − −
02 2
1 22z zur r
µπ
= − −
0 1 22 2
1 22x x x xw
r rµπ
− −= −
( )1 22 1
1
ln2
rz xr
µφ θ θπ
= − + −
( )1 2 12 1 2 2
2 12x xu zr r
µ θ θπ
= − − − −
( ) ( )2 22 11 22 2
1 2 1
ln2
x x x z x x x zrwr r r
µπ
− + − += − + −
( ) ( )( )2
2 22 21 2 2 1 2
1
log2
rz x x x z xr
µφ θ θπ
= − − − − −
( )2 2 2
2 2 1 22 12 2 2
2 1 1
log 22
x x ru z xr r r
µ θ θπ
= − − − −
( ) ( ) ( )
( )
2 21 1 2 22
1 2 2 21 2
22
2 1 21
[22
2 log ]
x x x x x xw x x
r r
rz xr
µπ
θ θ
− −= − + − +
+ − −
Vór
tice
( ) ( )01 1 2 2[ ln
2rx x x x zγφ θ θ
π= − − − − −
( )02 12
u γ θ θπ
= −
0 2
1
ln2
rwr
γπ
=
( )2 2 2
1 2 22 1 2
12 2 2
11
[ ln2 2 2
]2
r x x zzxz x xr
x x z
γφ θπ
θ
− −= + − + +
− −−
( )1 22 1
1
ln2
ru x zr
γ θ θπ
= − +
( ) ( )1 22 1 2 1
1
ln2
rw x x z xr
γ θ θπ
= − − − +
( ) ( )
( )( ) ( )
2 2 3 322 1 2 1 1 1 2 2
22 2 2 2 2
2 1 21
[4 2 212
+2 3 3 log ]
xz x x z x x x x
rx x z z x zr
γφ θ θπ
θ θ
= − + − + − +
− − + −
( )( ) ( )2
2 22 22 1 2 1 2
1
log12
ru x z z x x xzr
γ θ θπ
= − − + − +
Tabela A.1 – Tipos de distribuições de singularidades.
107
APÊNDICE B
GEOMETRIA DE AEROFÓLIOS NACA DE 4 DÍGITOS
Através de um trabalho experimental e de mapeamento dos perfis mais comumente
utilizados, o NACA, antigamente um setor da NASA, desenvolveu vários perfis com
nomenclatura específica. A nomenclatura consta em uma série de 4, 5 ou 6 dígitos
representado as propriedades do aerofólio. As propriedades geométricas de um aerofólio
estão representadas na Figura B.1.
Figura B.1 – Parâmetros geométricos de um aerofólio NACA de 4 dígitos e curvas da linha média de curvatura e contorno. (ABBOT & VON DOEHOFF, 1959)
x
z
x=0 (bordo de ataque)
x=c (bordo de fuga)
corda c
linha da corda
linha média de curvatura
curvatura máxima (m)
espessura máxima (t)
Raio do bordo de ataque (rt)
posição x da espessura
máxima (p)
posição x da curvatura máxima
108
As equações utilizadas para calcular a linha média de curvatura e a distribuição de
espessura estão representadas na Tabela B.1. Para calcular o contorno do aerofólio, basta
adicionar estas duas contribuições.
Tabela B.1 - Equações para Determinar a geometria do aerofólio
Para a série de 4 dígitos, basta analisar os dígitos do perfil conforme a ilustrado na
Figura B.2 e substiruir os valores de p e m correspondentes das equações da tabela B.1.
.
Espessura máxima (100 x t )
posição x da curvatura máxima ( 10 x p )
Curvatura máxima ( 100 x m )
NACA X1 X2 X3 X4
Figura B.2 – Nomenclatura NACA de 4 dígitos
Linha Média de Curvatura
( ) ( )2
21 2 xpxpmxy −= 0 ≤ x ≤ p
( ) ( ) ( )2
22 2211
xpxpp
mxy −+−−
= p < x ≤ 1
Contorno do Aerofólio
( ) ( )432 10150,028430,035160,012600,029690,02,0
xxxxxtxyt −+−−±=
Raio do Bordo de Ataque
rt = 1,1019 t2
109
APÊNDICE C
RESULTADOS PARA O AEROFÓLIO NACA 0012
COM REYNOLDS = 1,7 X 105
Os resultados aqui apresentados consistem nas simulações numéricas para o
aerofólio NACA0012, com parâmetros definidos através da Tabela C.1. As simulações
foram todas realizadas com esquema de difusão de Avanço Randômico com a mesma
semente, esquema de Adams-Bashforth de segunda ordem para o avanço temporal e
reflexão dos vórtices que penetram no corpo. Os resultados listados compreendem as
seguintes figuras para cada ângulo de ataque estudado:
• Posição final da nuvem de vórtices ao final da simulação, ou seja, a esteira do
aerofólio;
• Linhas de corrente em uma região próxima do aerofólio, também ao final da
simulação. A proximidade das linhas de corrente não representa necessariamente
maiores velocidades, pois algumas regiões possuem linhas mais próximas apenas
para aumentar a resolução da visualização;
• Distribuição de pressão média sobre o aerofólio calculada ao longo de toda
simulação;
• Carregamentos de arrasto e sustentação ao longo do tempo.
Tabela C.1 – Parâmetros Numéricos Empregados
Parâmetro Símbolo Valor Empregado
Número de Reynolds Re 1,7 x 105
Incremento de tempo ∆t 0,025
Número de painéis N 300
Número de sub-painéis Nsub 5
Distância de geração ε 0,005
Raio do núcleo dos vórtices σ0 0,005
Número de Passos no tempo NT 400
110NACA 0012 - =α 0o
Reynolds = 1,7 x 105
Figura C.1 – Esteira ao final da simulação.
Figura C.2 – Linhas de corrente ao final da simulação.
111
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 0012 - =α 0o
Reynolds = 1,7 x 105
Figura C.3 – Distribuição de pressão média ao longo do aerofólio.
Figura C.4 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
112
NACA 0012 - =α 2o
Reynolds = 1,7 x 105
Figura C.5 – Esteira ao final da simulação.
Figura C.6 – Linhas de corrente ao final da simulação.
113
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 0012 - =α 2o
Reynolds = 1,7 x 105
Figura C.7 – Distribuição de pressão média ao longo do aerofólio.
Figura C.8 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
Cl Cd
Cp médio Contorno
-Cp
X
t
114
NACA 0012 - =α 4o
Reynolds = 1,7 x 105
Figura C.9 – Esteira ao final da simulação.
Figura C.10 – Linhas de corrente ao final da simulação.
115
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 0012 - =α 4o
Reynolds = 1,7 x 105
Figura C.11 – Distribuição de pressão média ao longo do aerofólio.
Figura C.12 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
116
NACA 0012 - =α 6o
Reynolds = 1,7 x 105
Figura C.13 – Esteira ao final da simulação.
Figura C.14 – Linhas de corrente ao final da simulação.
117
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 0012 - =α 6o
Reynolds = 1,7 x 105
Figura C.15 – Distribuição de pressão média ao longo do aerofólio.
Figura C.16 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
118
NACA 0012 - =α 8o
Reynolds = 1,7 x 105
Figura C.17 – Esteira ao final da simulação.
Figura C.18 – Linhas de corrente ao final da simulação.
119
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 0012 - =α 8o
Reynolds = 1,7 x 105
Figura C.19 – Distribuição de pressão média ao longo do aerofólio.
Figura C.20 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
Cl Cd
Cp médio Contorno
-Cp
X
t
120NACA 0012 - =α 10o
Reynolds = 1,7 x 105
Figura C.21 – Esteira ao final da simulação.
Figura C.22 – Linhas de corrente ao final da simulação.
121
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 0012 - =α 10o
Reynolds = 1,7 x 105
Figura C.23 – Distribuição de pressão média ao longo do aerofólio.
Figura C.24 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
122
NACA 0012 - =α 12o
Reynolds = 1,7 x 105
Figura C.25 – Esteira ao final da simulação.
Figura C.26 – Linhas de corrente ao final da simulação.
123
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 0012 - =α 12o
Reynolds = 1,7 x 105
Figura C.27 – Distribuição de pressão média ao longo do aerofólio.
Figura C.28 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
Cl Cd
Cp médio Contorno
-Cp
X
t
124
NACA 0012 - =α 14o
Reynolds = 1,7 x 105
Figura C.29 – Esteira ao final da simulação.
Figura C.30 – Linhas de corrente ao final da simulação.
125
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 0012 - =α 14o
Reynolds = 1,7 x 105
Figura C.31 – Distribuição de pressão média ao longo do aerofólio.
Figura C.32 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
126
NACA 0012 - =α 16o
Reynolds = 1,7 x 105
Figura C.33 – Esteira ao final da simulação.
Figura C.34 – Linhas de corrente ao final da simulação.
127
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 0012 - =α 16o
Reynolds = 1,7 x 105
Figura C.35 – Distribuição de pressão média ao longo do aerofólio.
Figura C.36 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
128
NACA 0012 - =α 18o
Reynolds = 1,7 x 105
Figura C.37 – Esteira ao final da simulação.
Figura C.38 – Linhas de corrente ao final da simulação.
129
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 0012 - =α 18o
Reynolds = 1,7 x 105
Figura C.39 – Distribuição de pressão média ao longo do aerofólio.
Figura C.40 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
130
NACA 0012 - =α 20o
Reynolds = 1,7 x 105
Figura C.41 – Esteira ao final da simulação.
Figura C.42 – Linhas de corrente ao final da simulação.
131
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9
0
0.5
1
1.5
2
2.5
3
t
ClCd
NACA 0012 - =α 20o
Reynolds = 1,7 x 105
Figura C.43 – Distribuição de pressão média ao longo do aerofólio.
Figura C.44 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cl Cd
Cp médio Contorno
X
t
132
NACA 0012 - =α 25o
Reynolds = 1,7 x 105
Figura C.45 – Esteira ao final da simulação.
Figura C.46 – Linhas de corrente ao final da simulação.
133
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9
0
0.5
1
1.5
2
2.5
3
t
ClCd
NACA 0012 - =α 25o
Reynolds = 1,7 x 105
Figura C.47 – Distribuição de pressão média ao longo do aerofólio.
Figura C.48 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
134
NACA 0012 - =α 30o
Reynolds = 1,7 x 105
Figura C.49 – Esteira ao final da simulação.
Figura C.50 – Linhas de corrente ao final da simulação.
135
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9
0
0.5
1
1.5
2
2.5
3
t
ClCd
NACA 0012 - =α 30o
Reynolds = 1,7 x 105
Figura C.51 – Distribuição de pressão média ao longo do aerofólio.
Figura C.52 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
136NACA 0012 - =α 45o
Reynolds = 1,7 x 105
Figura C.53 – Esteira ao final da simulação.
Figura C.54 – Linhas de corrente ao final da simulação.
137
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9
0
0.5
1
1.5
2
2.5
3
t
ClCd
NACA 0012 - =α 45o
Reynolds = 1,7 x 105
Figura C.55 – Distribuição de pressão média ao longo do aerofólio.
Figura C.56 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
138NACA 0012 - =α 60o
Reynolds = 1,7 x 105
Figura C.57 – Esteira ao final da simulação.
Figura C.58 – Linhas de corrente ao final da simulação.
139
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
0
1
2
3
4
5
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9
0
1
2
3
4
5
6
t
ClCd
NACA 0012 - =α 60o
Reynolds = 1,7 x 105
Figura C.59 – Distribuição de pressão média ao longo do aerofólio.
Figura C.60 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
140
NACA 0012 - =α 75o
Reynolds = 1,7 x 105
Figura C.61 – Esteira ao final da simulação.
Figura C.62 – Linhas de corrente ao final da simulação.
141
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
0
1
2
3
4
5
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9
0
1
2
3
4
5
6
t
ClCd
NACA 0012 - =α 75o
Reynolds = 1,7 x 105
Figura C.63 – Distribuição de pressão média ao longo do aerofólio.
Figura C.64 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
142
APÊNDICE D
RESULTADOS PARA O AEROFÓLIO NACA0012
COM REYNOLDS = 1,7 X 105 E ÂNGULO DE ATAQUE 90O
O resultado deste apêndice é para o aerofólio NACA0012, com parâmetros
definidos através da Tabela D.1. Para o ângulo de ataque de 90o foi necessário um
tratamento especial, requerendo uma simulação de um tempo maior, até que fosse
identificado o caráter oscilatório da esteira, similar ao comportamento de um corpo
rombudo. Como o tempo computacional tornava-se muito elevado, foi necessário também
aumentar o valor do incremento de tempo. A simulação foi realizada com esquema de
difusão de Avanço Randômico, esquema de Adams-Bashforth de segunda ordem para o
avanço temporal e reflexão dos vórtices que penetram no corpo.
Tabela D.1 – Parâmetros Numéricos Empregados
Parâmetro Símbolo Valor Empregado
Número de Reynolds Re 1,7 x 106
Incremento de tempo ∆t 0,05
Número de painéis N 300
Número de sub-painéis Nsub 5
Distância de geração ε 0,005
Raio do núcleo dos vórtices σ0 0,005
Número de Passos no tempo NT 800
143
Figura D.1 – Esteira ao final da simulação. Figura D.2 – Linhas de corrente ao final da simulação.
NACA 0012 - =α 90o
Reynolds = 1,7 x 105
144
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-2
-1
0
1
2
3
4
5
X
-cp
Cp medioContorno
0 5 10 15 20 25 30 35
-1
0
1
2
3
4
5
t
ClCd
NACA 0012 - =α 90o
Reynolds = 1,7 x 105
Figura D.3 – Distribuição de pressão média ao longo do aerofólio.
Figura D.4 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
145
APÊNDICE E
RESULTADOS PARA O AEROFÓLIO NACA0012
COM REYNOLDS = 1,7 X 106
Os resultados aqui apresentados consistem nas simulações numéricas para o
aerofólio NACA0012, com parâmetros definidos através da Tabela E.1. As simulações
foram todas realizadas com esquema de difusão de Avanço Randômico com a mesma
semente, esquema de Adams-Bashforth de segunda ordem para o avanço temporal e
reflexão dos vórtices que penetram no corpo. Os resultados listados compreendem as
seguintes figuras para cada ângulo de ataque estudado:
• Posição final da nuvem de vórtices ao final da simulação, ou seja, a esteira do
aerofólio;
• Linhas de corrente em uma região próxima do aerofólio, também ao final da
simulação. A proximidade das linhas de corrente não representa necessariamente
maiores velocidades, pois algumas regiões possuem linhas mais próximas apenas
para aumentar a resolução da visualização;
• Distribuição de pressão média sobre o aerofólio calculada ao longo de toda
simulação;
• Carregamentos de arrasto e sustentação ao longo do tempo.
Tabela E.1 – Parâmetros Numéricos Empregados
Parâmetro Símbolo Valor Empregado
Número de Reynolds Re 1,7 x 106
Incremento de tempo ∆t 0,025
Número de painéis N 300
Número de sub-painéis Nsub 5
Distância de geração ε 0,005
Raio do núcleo dos vórtices σ0 0,005
Número de Passos no tempo NT 400
146NACA 0012 - =α 16o
Reynolds = 1,7 x 106
Figura E.1 – Esteira ao final da simulação.
Figura E.2 – Linhas de corrente ao final da simulação.
147
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9
0
0.5
1
1.5
2
2.5
3
t
ClCd
NACA 0012 - =α 16o
Reynolds = 1,7 x 106
Figura E.3 – Distribuição de pressão média ao longo do aerofólio.
Figura E.4 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
148NACA 0012 - =α 18o
Reynolds = 1,7 x 106
Figura E.5 – Esteira ao final da simulação.
Figura E.6 – Linhas de corrente ao final da simulação.
149
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9
0
0.5
1
1.5
2
2.5
3
t
ClCd
NACA 0012 - =α 18o
Reynolds = 1,7 x 106
Figura E.7 – Distribuição de pressão média ao longo do aerofólio.
Figura E.8 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
150
APÊNDICE F
RESULTADO PARA O AEROFÓLIO NACA0012 COM REYNOLDS =
1,7 X 105 E TEMPO FINAL DE SIMULAÇÃO MAIS LONGO
Os resultados aqui apresentados consistem nas simulações numéricas para o
aerofólio NACA0012, com parâmetros definidos através da Tabela F.1. A simulação foi
realizada com esquema de difusão de Avanço Randômico, esquema de Adams-Bashforth
de segunda ordem para o avanço temporal, reflexão dos vórtices que penetram no corpo e
um número maior de passos de tempo que o utilizado anteriormente. Os resultados listados
consistem nas seguintes figuras:
• Posição final da nuvem de vórtices ao final da simulação, ou seja, a esteira do
aerofólio;
• Linhas de corrente em uma região próxima do aerofólio, também ao final da
simulação. A proximidade das linhas de corrente não representa necessariamente
maiores velocidades, pois algumas regiões possuem linhas mais próximas apenas
para aumentar a resolução da visualização;
• Distribuição de pressão média sobre o aerofólio calculada ao longo de toda
simulação;
• Carregamentos de arrasto e sustentação ao longo do tempo.
Tabela F.1 – Parâmetros Numéricos Empregados
Parâmetro Símbolo Valor Empregado
Número de Reynolds Re 1,7 x 105
Incremento de tempo ∆t 0,025
Número de painéis N 300
Número de sub-painéis Nsub 5
Distância de geração ε 0,005
Raio do núcleo dos vórtices σ0 0,005
Número de Passos no tempo NT 600
151NACA 0012 - =α 20o
Reynolds = 1,7 x 105
Figura F.1 – Esteira ao final da simulação.
Figura F.2 – Linhas de corrente ao final da simulação.
152
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 2 4 6 8 10 12 14
0
0.5
1
1.5
2
2.5
3
t
ClCd
NACA 0012 - =α 20o
Reynolds = 1,7 x 105
Figura F.3 – Distribuição de pressão média ao longo do aerofólio.
Figura F.4 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
153
APÊNDICE G
RESULTADOS PARA O AEROFÓLIO NACA 4415
COM REYNOLDS = 1,7 X 105
Os resultados aqui apresentados consistem nas simulações numéricas para o
aerofólio NACA4415, com parâmetros definidos através da Tabela G.1. As simulações
foram todas realizadas com esquema de difusão de Avanço Randômico com a mesma
semente, esquema de Adams-Bashforth de segunda ordem para o avanço temporal e
reflexão dos vórtices que penetram no corpo. Os resultados listados compreendem as
seguintes figuras para cada ângulo de ataque estudado:
• Posição final da nuvem de vórtices ao final da simulação, ou seja, a esteira do
aerofólio;
• Linhas de corrente em uma região próxima do aerofólio, também ao final da
simulação. A proximidade das linhas de corrente não representa necessariamente
maiores velocidades, pois algumas regiões possuem linhas mais próximas apenas
para aumentar a resolução da visualização;
• Distribuição de pressão média sobre o aerofólio calculada ao longo de toda
simulação;
• Carregamentos de arrasto e sustentação ao longo do tempo.
Tabela F.1 – Parâmetros Numéricos Empregados
Parâmetro Símbolo Valor Empregado
Número de Reynolds Re 1,7 x 105
Incremento de tempo ∆t 0,025
Número de painéis N 300
Número de sub-painéis Nsub 5
Distância de geração ε 0,005
Raio do núcleo dos vórtices σ0 0,005
Número de Passos no tempo NT 400
154NACA 4415 - =α 0o
Reynolds = 1,7 x 105
Figura G.1 – Esteira ao final da simulação.
Figura G.2 – Linhas de corrente ao final da simulação.
155
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 4415 - =α 0o
Reynolds = 1,7 x 105
Figura G.3 – Distribuição de pressão média ao longo do aerofólio.
Figura G.4 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
156NACA 4415 - =α 2o
Reynolds = 1,7 x 105
Figura G.5 – Esteira ao final da simulação.
Figura G.6 – Linhas de corrente ao final da simulação.
157
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 4415 - =α 2o
Reynolds = 1,7 x 105
Figura G.7 – Distribuição de pressão média ao longo do aerofólio.
Figura G.8 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
158NACA 4415 - =α 4o
Reynolds = 1,7 x 105
Figura G.9 – Esteira ao final da simulação.
Figura G.10 – Linhas de corrente ao final da simulação.
159
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 4415 - =α 4o
Reynolds = 1,7 x 105
Figura G.11 – Distribuição de pressão média ao longo do aerofólio.
Figura G.12 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
160
NACA 4415 - =α 6o
Reynolds = 1,7 x 105
Figura G.13 – Esteira ao final da simulação.
Figura G.14 – Linhas de corrente ao final da simulação.
161
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 4415 - =α 6o
Reynolds = 1,7 x 105
Figura G.15 – Distribuição de pressão média ao longo do aerofólio.
Figura G.16 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
162NACA 4415 - =α 8o
Reynolds = 1,7 x 105
Figura G.17 – Esteira ao final da simulação.
Figura G.18 – Linhas de corrente ao final da simulação.
163
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 4415 - =α 8o
Reynolds = 1,7 x 105
Figura G.19 – Distribuição de pressão média ao longo do aerofólio.
Figura G.20 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
164
NACA 4415 - =α 10o
Reynolds = 1,7 x 105
Figura F.21 – Esteira ao final da simulação.
Figura F.22 – Linhas de corrente ao final da simulação.
165
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 4415 - =α 10o
Reynolds = 1,7 x 105
Figura G.23 – Distribuição de pressão média ao longo do aerofólio.
Figura G.24 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
166
NACA 4415 - =α 12o
Reynolds = 1,7 x 105
Figura G.25 – Esteira ao final da simulação.
Figura G.26 – Linhas de corrente ao final da simulação.
167
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 4415 - =α 12o
Reynolds = 1,7 x 105
Figura G.27 – Distribuição de pressão média ao longo do aerofólio.
Figura G.28 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
168
NACA 4415 - =α 14o
Reynolds = 1,7 x 105
Figura G.29 – Esteira ao final da simulação.
Figura G.30 – Linhas de corrente ao final da simulação.
169
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 4415 - =α 14o
Reynolds = 1,7 x 105
Figura G.31 – Distribuição de pressão média ao longo do aerofólio.
Figura G.32 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
170
NACA 4415 - =α 16o
Reynolds = 1,7 x 105
Figura G.33 – Esteira ao final da simulação.
Figura G.34 – Linhas de corrente ao final da simulação.
171
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
t
ClCd
NACA 4415 - =α 16o
Reynolds = 1,7 x 105
Figura G.35 – Distribuição de pressão média ao longo do aerofólio.
Figura G.36 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
172
NACA 4415 - =α 18o
Reynolds = 1,7 x 105
Figura G.37 – Esteira ao final da simulação.
Figura G.38 – Linhas de corrente ao final da simulação.
173
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9
0
0.5
1
1.5
2
2.5
3
t
ClCd
NACA 4415 - =α 18o
Reynolds = 1,7 x 105
Figura G.39 – Distribuição de pressão média ao longo do aerofólio.
Figura G.40 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
174
NACA 4415 - =α 20o
Reynolds = 1,7 x 105
Figura G.41 – Esteira ao final da simulação.
Figura G.42 – Linhas de corrente ao final da simulação.
175
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.5
0
0.5
1
1.5
2
2.5
3
X
-cp
Cp medioContorno
0 1 2 3 4 5 6 7 8 9
0
0.5
1
1.5
2
2.5
3
t
ClCd
NACA 4415 - =α 20o
Reynolds = 1,7 x 105
Figura G.43 – Distribuição de pressão média ao longo do aerofólio.
Figura G.44 – Coeficientes de sustentação (Cl ) e arrasto (Cd) em função do tempo.
-Cp
Cp médio Contorno
Cl Cd
X
t
