Upload
others
View
3
Download
0
Embed Size (px)
Citation preview
Construção de uma Ferramenta Computacional para aSimulação Aerodinâmica de um conjunto Asa e Flap
Jorge André Fernandes Baginha
Dissertação para a obtenção de Grau de Mestre em
Engenharia Aeroespacial
Júri
Presidente: Professor Doutor João Manuel Lage de Miranda LemosOrientador: Professor Doutor Fernando José Parracho LauCo-orientador: Professor Doutor Filipe Szolnoky Ramos Pinto CunhaVogais: Professor Doutor João Manuel Melo de Sousa
Outubro 2013
Agradecimentos
O espaco limitado desta seccao de agradecimentos, seguramente, nao me permite agradecer, como
devia, a todas as pessoas que, ao longo do mestrado em Engenharia Aeroespacial e no processo de
elaboracao desta tese me ajudaram, directa ou indirectamente, a concretizar mais esta etapa da minha
formacao academica.
Desta forma, deixo apenas algumas palavras, poucas, mas um profundo sentimento de agradeci-
mento.
Aos Professores Fernando Lau e Filipe Cunha (Orientador e Co-orientador, respectivamente),
agradeco a disponibilidade, o incentivo, a sabedoria e os ensinamentos constantes em todo o processo
de orientacao cientıfica desta dissertacao.
A Optimal Structural Solutions, em especial, ao Engº Antonio Reis, agradeco toda a disponibilidade
demonstrada e o apoio prestado.
Aos meus amigos, em especial ao Miguel e a Marta, pelas constantes manifestacoes de interesse
e encorajamento, pela preocupacao, pelos momentos de descontraccao e companheirismo.
A Sara, minha namorada, ouvinte atenta de algumas duvidas, inquietacoes, desanimos e sucessos,
pelo apoio e transmissao de confianca, pela dedicacao e carinhos diarios, um agradecimento especial.
A minha famılia, em especial aos meus Pais, um enorme obrigada por acreditarem sempre em
mim e por todos os ensinamentos de vida, pela compreensao inestimaveis, pelos diversos sacrifıcios
suportados. Espero que esta etapa, que agora termino, possa, de alguma forma, retribuir e compensar
todo o carinho, apoio e dedicacao que, constantemente, me oferecem. A eles, dedico todo este trabalho.
ii
Resumo
O presente trabalho tem como principal objectivo a construcao de uma ferramenta computacional para
a previsao das caracteristicas aerodinamicas de um conjunto asa e flap.
Utilizou-se o metodo dos paineis 3D, com a condicao de fronteira de Dirichlet, para a solucao do
escoamento invıscido. Neste escoamento invıscido e possivel contabilizar o efeito solo atraves do
metodo das imagens.
Para a contabilizacao dos efeitos viscosos a geometria e dividida em seccoes bidimensionais onde
por sua vez sao analisadas individualmente pelo XFOIL. O acoplamento das duas solucoes e feito por
intermedio de uma velocidade de transpiracao.
As solucoes do escoamento invıscido e viscoso foram comparadas com resultados experimentais
obtidos para dois tipos de asas diferentes e um conjunto asa e flap.
Os resultados demonstraram-se satisfatorios na medida em que os erros relativos do metodo com-
putacional relativamente ao procedimento experimental foram pequenos.
Concluımos que o codigo desenvolvido pode constituir, numa primeira aproximacao, uma valiosa
ferramenta computacional a ser utilizada no projecto de velas rigidas de catamaras
Palavras-chave: Metodo dos Paineis 3D, Condicao de fronteira de
Dirichlet, Efeito Solo, Acoplamento dos Efeitos Viscosos
iv
Abstract
The main purpose of the present work was to develop a computational tool used to predict the aeorody-
namic characteristics of a wing-flap.
The 3D panels method using the Dirichlet-type boundary condition was adopted to obtain the inviscid
flow. The latter allowed the consideration of the ground effect through the use of the images method.
In order to consider the viscid effects, the geometry was discretized/divided into bidimensional cross
sections. Each of them was individualy analysed by XFOIL program/routine. The coupling of both
solutions was done by means of a transpiration velocity.
The inviscid and viscous flows solutions were compared with experimental data obtained for two
different types of wings and for a wing-flap aggregate.
The results were satisfactory in the sense that the relative erros between the computational and
experimental data were small.
We concluded that the developed code may constitute, as a first approach, a valuable computational
tool to be used in the design of rigid sails.
Keywords: 3D Panel Method, Dirichlet Boundary Condition, Ground
Effect, Viscous Effects
vi
Conteudo
Agradecimentos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Resumo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi
Lista de Sımbolos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii
1 Introducao 1
1.1 Objectivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Historia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Estrutura da Dissertacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2 Metodo dos Paineis 3D 6
2.1 Formulacao Teorica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Condicoes de Fronteira . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.1 Condicao de Fronteira de Dirichlet . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.3 Forma e intensidade da esteira . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.4 Procedimento Numerico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4.1 Discretizacao da superfıcie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.5 Calculo das Forcas Aerodinamicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3 Efeito Solo 20
3.1 Aplicacao do metodo das imagens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2 Camada Limite Atmosferica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4 Acoplamento dos Efeitos Viscosos 25
4.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2 Esquemas de acoplamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.2.1 Metodo Directo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2.2 Metodo Inverso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2.3 Metodo Semi-Inverso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.2.4 Metodo Simultaneo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.3 Interaccao com o XFOIL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.3.1 Asa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
viii
4.3.2 Asa com Flap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.4 Procedimento Numerico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.5 Calculo das Forcas Aerodinamicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
5 Validacao e Discussao dos Resultados 37
5.1 Estudo de Convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.2 Efeito Solo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.3 Asa na Vertical . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.4 Asa com Flap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.5 Comparacao com Resultados Experimentais . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.5.1 Asa Rectangular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.5.2 Asa com Afilamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.5.3 Asa RAE com Flap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6 Conclusoes e Trabalhos Futuros 63
Bibliografia 64
ix
Lista de Tabelas
3.1 Rugosidade da superfıcie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.1 Coeficiente de sustentacao e resistencia para as quatro malhas . . . . . . . . . . . . . . 39
5.2 Variacao do coeficiente de sustentacao e resistencia para uma asa rectangular com uma
superfıcie solida junto a uma das extremidades . . . . . . . . . . . . . . . . . . . . . . . . 48
5.3 Resultado do coeficiente de sustentacao experimental e invıscido com o respectivo erro
relativo para um Reynolds 3× 106 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.4 Resultado experimental, invıscido e viscoso acoplado ao invıscido para um Reynolds
1.92× 106 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.5 Dados da asa com afilamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.6 Resultado experimental, invıscido e viscoso acoplado ao invıscido para um Reynolds
4.4× 106 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.7 Parametros relativamente a distancia do flap ao bordo de fuga da asa para dois angulos
de deflexao, 10° e 25° . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.8 Resultados experimentais e numericos para uma deflexao do flap de 10° . . . . . . . . . 60
5.9 Resultados experimentais e numericos para uma deflexao do flap de 25° . . . . . . . . . 61
x
Lista de Figuras
1.1 O barco actual da Team Cascais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.1 Escoamento potencial em torno de um corpo solido. . . . . . . . . . . . . . . . . . . . . . 7
2.2 Aplicacao da condicao de Kutta. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Discretizacao de um corpo em paineis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.4 Relacao entre o painel da esteira e os dois paineis do bordo de fuga. . . . . . . . . . . . 12
2.5 Construcao do painel plano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.6 Aparecimento de fugas na construcao dos paineis. . . . . . . . . . . . . . . . . . . . . . . 15
2.7 Painel a funcionar como a) fonte b) dipolo. . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.8 Esquema de paineis para o calculo da velocidade no painel k. . . . . . . . . . . . . . . . 17
2.9 Esquema de paineis numa extremidade para o calculo da velocidade no painel k. . . . . 17
2.10 Esquema do aparecimento da resistencia induzida para uma seccao transversal de uma
asa tri-dimensional. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1 Esquema da velocidade induzida pelos vortices das extremidades da imagem da asa . . 21
3.2 Asa e imagem, com a representacao da sequencia da aplicacao da condicao de fronteira 22
3.3 Esquema de uma superfıcie sustentadora a operar proxima do solo com uma camada
limite atmosferica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.1 Decomposicao do escoamento em duas regioes: Invıscida e Viscosa . . . . . . . . . . . 25
4.2 Perfil de velocidade para um escoamento a) real e b) invıscido. . . . . . . . . . . . . . . . 26
4.3 Espessura de deslocamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.4 Alteracao nas linhas de corrente, com o objectivo de simular os efeitos viscosos, atraves
de uma velocidade de transpiracao. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.5 Metodo directo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.6 Metodo inverso. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.7 Metodo simultaneo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.8 As duas primeiras seccoes para a analise viscosa . . . . . . . . . . . . . . . . . . . . . . 32
4.9 Seccao do Flap para analise no XFOIL . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.10 Deslocamento da linha de corrente atraves de uma velocidade normal a superficie . . . 35
xi
5.1 Quatro malhas para o estudo da convergencia: a) M1=10x5, b) M2=20x11, c) M3=40x21,
d) M4=80x41 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.2 Coeficientes de pressao Cp em funcao da percentagem de corda %c: a) α = 0° , b)
α = 5° , c) α = 10° , d) α = 15° . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.3 a) erro relativo CL Vs Nº Paineis b) erro relativo CD Vs Nº Paineis c) Tempo Vs Nº Paineis 40
5.4 Representacao da asa e imagem para simular o efeito solo pelo metodo dos paineis 3D . 41
5.5 Comparacao dos resultados obtidos pelo metodo dos paineis com resultados STAR-
CCM+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.6 CL Vs α. a) STAR-CCM+, b) Metodo dos paineis . . . . . . . . . . . . . . . . . . . . . . . 43
5.7 Distribuicao do coeficiente de pressao a meia envergadura. a) α = 0 e h/c = 1 b) α = 0
e h/c = 0.1 c) α = 8 e h/c = 1 d)α = 8 e h/c = 0.1 . . . . . . . . . . . . . . . . . . . . . . 44
5.8 CL Vs CD. a) STAR-CCM+, b) Metodo dos paineis . . . . . . . . . . . . . . . . . . . . . . 45
5.9 Asa na vertical . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.10 CL Vs h/c, a) α = 2, b) α = 4, c) α = 6, d) α = 8 . . . . . . . . . . . . . . . . . . . . . . . . 47
5.11 Coeficiente de pressao nas duas configuracoes para a) asa e b) flap . . . . . . . . . . . 49
5.12 Ilustracao do conjunto asa e flap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.13 Distribuicao de Cp na asa e flap para a configuracao 1 . . . . . . . . . . . . . . . . . . . . 50
5.14 Distribuicao de Cp na asa e flap para a configuracao 2 . . . . . . . . . . . . . . . . . . . . 51
5.15 Perfil NACA0012 utilizado na construcao da asa rectangular . . . . . . . . . . . . . . . . 51
5.16 Malha na asa rectangular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.17 Coeficientes de pressao para duas seccoes a um angulo de ataque de 4.64º com Rey-
nolds igual a 3× 106 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.18 Coeficientes de pressao para duas seccoes a um angulo de ataque de 10.97º com Rey-
nolds igual a 3× 106 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.19 CL vs alpha para uma asa rectangular a um Reynolds igual a 3× 106 . . . . . . . . . . . 54
5.20 CL vs alpha para uma asa rectangular a um Reynolds igual a 1.92× 106 . . . . . . . . . . 55
5.21 Perfil NACA65-210 utilizado na construcao da asa com afiamento . . . . . . . . . . . . . 55
5.22 Esquema da asa com afilamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.23 Malha na asa com afilamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.24 a) CL Vs alpha, b) CL Vs CD. Para asa com afilamento . . . . . . . . . . . . . . . . . . . 57
5.25 Representacao das variaveis que definem a distancia do flap a asa. . . . . . . . . . . . . 58
5.26 Conjunto asa flap para simulacao, com a respectiva imagem . . . . . . . . . . . . . . . . 59
5.27 Conjunto asa flap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.28 Coeficientes de sustentacao e resistencia para um angulo de deflexao 10° para o flap. a)
CL vs α, b) CL vs CD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.29 Coeficientes de sustentacao e resistencia para um angulo de deflexao 25° para o flap. a)
CL vs α, b) CL vs CD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
xii
Lista de Sımbolos
α Angulo de ataque geometrico
αef Angulo de ataque efectivo
αi Angulo de ataque induzido
∆Sk Area do painel k
δ∗ Espessura de deslocamento
Γ Circulacao do vortice
∆Fk Forca no painel k
n Vector unitario normal a superficie
nk Vector unitario normal ao painel k
Q∞ Vector velocidade do escoamento de aproximacao
Qk Velocidade total no painel k
µ Distribuicao de intensidade do dipolo
µk Intensidade do dipolo do painel k do corpo
µl Intensidade do dipolo no painel l da esteira
Φ Potencial de velocidade de perturbacao
Φ∗ Potencial de velocidade exterior
Φ∞ Potencial de velocidade do escoamento de aproximacao
Φd Potencial de velocidade induzido por um painel a comportar-se como dipolo
Φf Potencial de velocidade induzido por um painel a comportar-se como fonte
Φ∗i Potencial de velocidade interior
σ Distribuicao de intensidade da fonte
σk Intensidade da fonte do painel k do corpo
xiii
Bk Coeficiente de influencia da fonte
CD Coeficiente de Resistencia
Ck Coeficiente de influencia do dipolo
CL Coeficiente de Sustentacao
Cpk Coeficiente de pressao no painel k
Dind Resistencia induzida
L Sustentacao
N Numero de paineis na superficie do corpo
Nw Numero de paineis na superficie da esteira
ql Velocidade tangencial ao painel na direccao l
qm Velocidade tangencial ao painel na direccao m
qn Velocidade normal ao painel
SB Superficie do corpo
U∞ Componente x do vector velocidade do escoamento de aproximacao
V∞ Componente y do vector velocidade do escoamento de aproximacao
W∞ Componente z do vector velocidade do escoamento de aproximacao
xiv
Capıtulo 1
Introducao
O objetivo principal deste estudo passa pela construcao de uma ferramenta computacional capaz de
prever as caracterısticas aerodinamicas de uma superfıcie sustentadora. A finalidade ultima da ferra-
menta sera a sua insercao no projecto de velas rıgidas (Wing Sails) de catamaras dedicados a alta
competicao. De facto, este trabalho surge no seguimento de uma parceira desenvolvida entre a Opti-
mal Structural Solution [1] e a Tony Castro [2] com vista a construcao do primeiro catamara totalmente
portugues a participar na International C-Class Catamaran Competition (ICCCC). Para o efeito, a Op-
timal Structural Solutions sugeriu ao Instituto Superior Tecnico a elaboracao, por parte de mestrandos
do curso de engenharia aeroespacial, de uma tese com o tema “Development of simulation tool for
catamaran solid sails”.
A ICCCC e uma competicao onde cada equipa participante desenvolve o seu proprio catamara,
sendo que, apesar de existirem algumas restricoes geometricas ao comprimento do catamara, largura
e area da vela, e conferido grande grau de liberdade no que diz respeito ao seu design, o que por sua
vez incentiva a procura e desenvolvimento de tecnologias avancadas incorporando domınios cientıficos
como a aerodinamica, hidrodinamica e materiais. Neste tipo de aparelho velico – velas rıgidas-, que
surgiu nos anos 80, o perfil e espesso contrariamente ao que acontece nas velas convencionais. A
existencia de espessura garante uma maior forca gerada e, deste modo, melhores prestacoes quando
comparadas com as obtidas pelas velas tradicionais. As velas rıgidas sao constituıdas por duas su-
perfıcies (Fig. 1.1): o elemento principal e o flap. Este ultimo serve para incluir curvatura na vela,
aumentando a forca resultante.
Dadas as condicoes de operacao da vela, considerou-se oportuno incluir na ferramenta computa-
cional quer o efeito da superfıcie da agua, quer a variacao da velocidade no interior da camada limite
atmosferica. E neste enquadramento que o presente trabalho sera desenvolvido.
1.1 Objectivo
Um metodo computacional que permita prever as caracterısticas aerodinamicas de uma forma economica
e eficiente numa asa, e uma ferramenta essencial nos dias de hoje.
1
Figura 1.1: O barco actual da Team Cascais
Apesar de existir uma enorme variedade de programas de Dinamica dos Fluidos Computacional
(Computational Fluid Dynamics - CFD) capazes de obter solucoes muito proximas da realidade, os
seus custos computacionais torna-os pouco atractivos para o utilizador, pelo menos para uma fase
inicial do projecto.
Contrariamente aos codigos CFD que discretizam todo o domınio para o desenvolvimento das
equacoes de Navier-Stokes, o metodo dos paineis permite a simulacao de um escoamento em torno
de um corpo discretizando apenas a sua superfıcie. Esta simplificacao faz do metodo dos paineis um
instrumento rapido na obtencao da solucao final. Naturalmente que os resultados nao terao um grau de
precisao tao elevado como sucederia perante um calculo em CFD. Porem, numa fase inicial de projecto,
em que nao e tao relevante o resultado exacto mas antes a tendencia de variacao de uma determinada
caracterıstica, a sua utilidade e inegavel. Ou seja, a aplicacao do metodo dos paineis e uma boa opcao
para uma analise preliminar, quando o custo computacional e um dos factores mais importantes. Ora,
uma vez que a ferramenta que sera desenvolvida ao longo deste trabalho sera utilizada numa fase
inicial de projecto, fica justificada a escolha do metodo dos paineis 3D.
Por outro lado, e considerando que a vela e composta por duas superficies sustentadoras - ele-
mento principal e flap -, o codigo a desenvolver tera de ser capaz de simular dois corpos no mesmo
escoamento.
Por fim, e para a simulacao da influencia da superfıcie de agua nas caracterısticas da vela, e tendo
2
em conta que o metodo utilizado e puramente invıscido, sera necessario o recurso ao metodo das
imagens.
O codigo a desenvolver sera implementado em linguagem MATLAB©(R2008b).
1.2 Historia
A primeira implementacao computacional do metodo dos paineis foi desenvolvida por Hess e Smith [3]
em 1962. Este metodo baseava-se na discretizacao do corpo em paineis planos a comportarem-se
como fontes com intensidades constantes e tinha a capacidade de simulacao de escoamentos tridi-
mensionais em torno de corpos complexos nao sustentadores. No ano 1972, Hess [4] adiciona os
elementos do tipo dipolo ao seu metodo, tornando possıvel o calculo de escoamentos tridimensionais
em torno de corpos sustentadores.
Em 1980 foi desenvolvido por Bristow [5] o codigo MCAIR no qual os paineis utilizados para modelar
o corpo eram considerados planos mas as intensidades das singularidades por painel, contrariamente
aos codigos anteriores, podiam variar. Foi o codigo pioneiro neste tipo de sigularidades de ordem
elevada. Com o surgimento deste codigo passou a ser possıvel a simulacao dos efeitos viscosos
atraves de um metodo iterativo: a cada iteracao a geometria era modificada com o objectivo de se
incluir a espessura de deslocamento da camada limite formada na superfıcie do corpo.
Um ano volvido, a Boing desenvolveu um codigo para a NASA com o nome de PAN AIR [6]. Neste
codigo o corpo era discretizado em paineis planos, cada qual dividido em cinco subpaineis nos quais
era permitida a variacao da intensidade da singularidade. Neste sentido era considerado um metodo
de singularidade de ordem elevada. Alem disso, o codigo era capaz de calcular escoamentos quer
subsonicos quer supersonicos.
Ate ao inıcio da decada de 80, a maioria dos codigos era dedicada a empresas ligadas a industria
espacial. No entanto, a rapida evolucao da tecnologia computacional permitiu as pequenas empresas
ligadas a outras industrias comecarem a usar este tipo de codigos.
O primeiro codigo comercial diponıvel as pequenas empresas surgiu em 1982 atraves de Maskew
com o nome VSAERO [7, 8]. Utilizava uma formulacao em que os paineis eram considerados planos
e as intensidades das singularidades constantes. Este codigo possibilitava o acoplamento dos efeitos
viscosos a solucao inviscida atraves de uma velocidade de transpiracao ao inves da tradicional alteracao
da geometria. Tambem foi adicionado ao codigo uma rotina para o calculo do enrolamento da esteira.
Em 1987 surgiu o PMARC [9], desenvolvido por Ashby. Trata-se de um codigo muito semelhante ao
VSAERO mas capaz de proceder ao calculo de escoamentos nao estacionarios.
Em 1988 Mark Drela desenvolveu um codigo denominado XFOIL [10], ainda hoje frequentemente
utilizado por ser uma ferramenta fiavel no que concerne a previsao de escoamentos invıscidos e vis-
cosos em torno de perfis bidimensionais. O facto de ter sido implementado em linguagem FORTRAN
facilita a interligacao deste codigo com ferramentas de optimizacao.
Nos anos seguintes os esforcos foram direcionados na melhoria do pre e pos processamento
(geracao da malha e apresentacao dos resultados) .
3
Com o aparecimento de diferentes tipos de codigos, estes ficam sujeitos a comparacoes. Em [11],
os autores efectuaram um estudo com o objectivo de compararem codigos que utilizam singularidades
de baixa e elevada ordem. A conclusao foi a de que os codigos desenvolvidos com a utilizacao de
singularidades de ordem elevada apresentam, para o mesmo numero de paineis, um maior grau de
precisao. Por outro lado, os codigos que utilizam as singularidades de baixa ordem revelaram-se, os
mais rapidos na obtencao da solucao.
Uma vez alcancado um estado de maturidade no que diz respeito as questoes teoricas do metodo
dos paineis, os esforcos comecaram a concentrar-se na procura de formas que permitissem a sua
integracao em processos iterativos de acoplamento dos efeitos viscosos.
Em 1997, Milewski [12] apresenta no seu trabalho um acoplamento de uma camada limite tridimen-
sional ao metodo dos paineis 3D de baixa ordem, o qual como se vera e feito atraves de um esquema
simultaneo.
Cabeci e Besnard [13], em 1998, apresentam um metodo de acoplamento, conhecido como semi-
inverso, de uma camada limite tridimensional a uma solucao obtida a partir do metodo dos paineis. Os
resultados apresentados foram satisfatorios tanto para uma unica asa imersa num escoamento assim
como para um conjunto de asa e flap.
Em 2000, Veldman et al., [14] apresentam um esquema de acoplamento, baseado no trabalho de
Veldman [15], denominado por quasi-simultaneo, sendo considerado um dos melhores metodos na
relacao qualidade da solucao/tempo computacional.
Em 2003, Andre Deperrois [16] construiu um codigo, XFLR5, com uma interface grafica user friendly,
que calcula escoamentos bidimensionais e tridimensionais em torno de perfis e asas, respectivamente.
No primeiro caso, o calculo e feito com recurso ao codigo do XFOIL. No segundo, o utilizador pode
escolher entre a teoria da linha sustentadora, metodo da malha de vortices ou metodo dos paineis 3D.
Neste codigo o acoplamento dos efeitos viscosos e efectuado atraves de uma analise bidimensional
das seccoes ao longo da envergadura do corpo.
Xwing [17], programa desenvolvido emMATLAB© e apresentado em 2008, possui a capacidade de
calculo de um escoamento potencial 3D com o acoplamento dos efeitos viscosos atraves da alteracao
da geometria do corpo com o objectivo de simular a espessura de deslocamento. Os autores assumi-
ram que o desenvolvimento da camada limite tridimensional podia ser aproximado a diversas seccoes
bidimensionais ao longo da envergadura. Para o calculo dos parametros das camadas limites bidimen-
sionais foi utilizado o programa XFOIL.
1.3 Estrutura da Dissertacao
A presente dissertacao encontra-se dıvida em seis capıtulos.
No capıtulo 2 sera feita uma explicacao detalhada do metodo dos paineis 3D, desde a sua formulacao
teorica ao calculo das forcas aerodinamicas geradas num corpo imerso num escoamento invıscido.
No capıtulo 3 sera exposta a via que foi seguida para introduzir o efeito de uma superfıcie solida na
proximidade de um corpo sustentador. Este efeito e conhecido como efeito solo e, para a sua simulacao,
4
utilizar-se-a o metodo das imagens.
O capıtulo 4 incidira sobre o processo necessario para se conseguir acoplar os efeitos viscosos a
um escoamento inviscido, utilizando uma velocidade de transpiracao.
O capıtulo 5 sera dedicado a validacao dos resultados obtidos pelo metodo dos paineis 3D assim
como o efeito solo e o acoplamento dos efeitos viscosos.
No ultimo capıtulo serao apresentadas as conclusoes dos resultados obtidos e fornecidas algumas
sugestoes para trabalhos futuros.
5
Capıtulo 2
Metodo dos Paineis 3D
O estudo de um escoamento de fluido real em torno de um corpo solido pode ser aproximado a um es-
coamento de fluido perfeito caso a interaccao viscosa/invıscida seja fraca. Para que tal aconteca, num
escoamento de fluido real, as espessuras de deslocamento δ∗ tem que ser pequenas comparativa-
mente a uma dimensao caracterıstica do corpo [18]. Sendo o metodo dos paineis um metodo numerico
puramente invıscido, este pode ser utilizado para prever o comportamento do escoamento exterior.
Neste capıtulo serao apresentadas todas as consideracoes tidas em conta para o calculo de um
escoamento inviscido em torno de um corpo sustentador pelo metodo dos paineis 3D. Para alem da
formulacao teorica do metodo serao abordadas questoes como o tipo de condicao de fronteira utilizada,
a modelacao da esteira de vortices libertada pelo sistema sustentador, discretizacao e, por fim, a forma
de calculo das forcas aerodinamicas resultantes no corpo.
Nesta seccao, o escoamento sera considerado invıscido e, portanto, sem a presenca da camada
limite. Posteriormente, e como se vera no capıtulo do Acopolamento dos Efeitos Viscosos, o efeito da
espessura de deslocamento sera tido em conta.
2.1 Formulacao Teorica
A explicacao que se segue da formulacao teorica do metodo dos paineis 3D tera como base as re-
ferencias [5, 8, 9].
Considerando um corpo em que e conhecida a sua fronteira SB sujeito a um escoamento potencial,
como mostra a Fig. 2.1, o escoamento exterior (regiao V) e incompressıvel, irrotacional e invıscido
fazendo com que o potencial de velocidade Φ∗ respeite a equacao de Laplace
O2Φ∗ = 0 (2.1)
Este potencial Φ∗ pode ser decomposto numa componente devido ao escoamento de aproximacao,
Φ∞, e num potencial de perturbacao, Φ, devido ao corpo.
Φ∗ = Φ + Φ∞ (2.2)
6
Figura 2.1: Escoamento potencial em torno de um corpo solido. [5]
A solucao da Eq. 2.1 pode ser construida, segundo Green [5, 18], pela soma de singularidades do
tipo fonte σ e dipolo µ com uma determinada distribuicao continua ao longo da superfıcie SB :
Φ∗(x, y, z) = − 1
4π
ˆ
SB
[σ
(1
r
)− µn · ∇
(1
r
)]dS + Φ∞ (2.3)
onde n e o vector unitario normal a superfıcie SB com a direccao para o interior do corpo, r a distancia
entre um ponto no espaco e o elemento dS e Φ∞ o potencial de velocidade do escoamento de aproximacao.
Φ∞ = U∞x+ V∞y +W∞z (2.4)
onde U∞, V∞ e W∞ sao as componentes x,y e z da velocidade do escoamento de aproximacao.
Em princıpio, e observando a Eq. 2.3, o escoamento em torno de um corpo sustentador pode ser
construıdo a partir de um numero infinito de combinacoes entre as distribuicoes de fontes e dipolos.
De forma a tornar essa combinacao unica para um determinado tipo de problema, e necessario definir
uma uma condicao de fronteira que exija que a velocidade normal a superfıcie do corpo seja igual a
zero ∂Φ∗/∂n = 0, tornando-a impermeavel. Porem, uma vez que a aplicacao da condicao de fronteira
nao e suficiente para a resolucao de problemas de escoamentos tridimensionais em torno de corpos
sustentadores, torna-se necessario introduzir condicoes adicionais que irao depender da fısica do pro-
blema. Estas condicoes estao relacionadas com a modelacao da esteira de vortices a jusante do corpo
sustentador e a relacao entre as intensidades dos dipolos na esteira com os que estao localizados na
superfıcie do corpo (condicao de Kutta).
Tendo em conta que as fontes sao usadas para simular o efeito da espessura do corpo e os dipolos
para a geracao de circulacao no sistema sustentador, e natural que a esteira seja modelada apenas
com dipolos uma vez que se considera que tem espessura desprezavel. Assim, a Eq. 2.3 pode ser
reescrita da seguinte forma:
Φ∗(x, y, z) =1
4π
ˆ
corpo+esteira
µn · ∇(
1
r
)dS − 1
4π
ˆ
corpo
σ
(1
r
)dS + Φ∞ (2.5)
7
Note-se que a formulacao na Eq. 2.5 tambem poderia ser feita sem atender a distribuicao de dipolos.
No entanto, tal impediria a simulacao de um corpo sustentador . A primeira implementacao numerica
desta formulacao, com verdadeiro sucesso, foi conseguida pelo matematico Hess e o aerodinamicista
Smith [3]. No entanto, so em 1972, quando Hess [4] adicionou a distribuicao de dipolos, e que se tornou
possıvel calcular escoamentos em torno de corpos sustentadores.
2.2 Condicoes de Fronteira
A condicao de fronteira, que garante que o corpo no seio do escoamento seja impermeavel, para a
Eq.2.1 pode ser aplicada de duas formas diferentes. Por um lado, especificando a velocidade normal
∂Φ∗/∂n = 0 a superfıcie SB igual a zero. Neste caso, onde a condicao de fronteira e aplicada de forma
directa, diz-se que se esta perante um problema de Neumann. Na segunda forma, conhecida por
problema de Dirichlet, o mesmo efeito e obtido por via da especificacao de um potencial de velocidade
na fronteira do corpo. A utilizacao desta condicao de fronteira e mais vantajosa a nıvel computacional
uma vez que se trabalha com uma quantidade escalar em vez de um vector velocidade.
Este factor, aliado a maior simplicidade da implementacao a nivel de programacao, levou a que a
condicao de fronteira escolhida neste trabalho fosse a de Dirichlet.
2.2.1 Condicao de Fronteira de Dirichlet
O potencial de velocidade no interior da superfıcie SB tambem respeita a equacao de Laplace [9] e
portanto a solucao para um ponto P qualquer no interior do corpo e igual a Eq. 2.5.
A condicao de que a velocidade normal a superfıcie do corpo tem que ser igual a zero em termos
de potencial de velocidade ∇ (Φ∗i ) · n = 0 resulta em Φ∗i = (Φ + Φ∞)i = const. Isto significa que, para
um corpo fechado, o potencial de velocidade no seu interior nao se altera.
Φ∗i (x, y, z) =1
4π
ˆ
corpo+esteira
µ∂
∂n
(1
r
)dS − 1
4π
ˆ
corpo
σ
(1
r
)dS + Φ∞ = const (2.6)
A Eq. 2.6 e a base dos problemas resolvidos pela condicao de fronteira Dirichlet.
Em [7], Maskew comparou os resultados das distribuicoes de pressoes e coeficientes de sustentacao
e resistencia, para um perfil simetrico NACA 0012 a um angulo de 10° com o escoamento de aproximacao,
para dois valores diferentes do potencial de velocidade no interior do corpo. Considerou-o igual a zero e
igual ao potencial de velocidade a infinito Φ∞. Para estes dois valores Maskew verificou que a diferenca
nas distribuicoes de pressoes e coeficientes de sustentacao e resistencia nao eram relevantes para
uma distribuicao simetrica (mesmo numero de paineis no extradorso e intradorso) de 38 e 78 paineis.
Por outro lado, para uma distribuicao simetrica de apenas 8 paineis, os resultados obtidos com a cons-
tante igual ao potencial de velocidade a infinito Φ∞ apresentam valores mais proximos dos resultados
com 78 paineis. A utilizacao do potencial de velocidade no interior do corpo igual ao potencial a infinito
tambem mostrou-se superior numa situacao em que a distribuicao de paineis nao era simetrica. Isto
8
especialmente na zona do bordo de fuga.
Subjacente a esta melhoria de comportamentos esta o facto da solucao (µ) ter uma intensidade
menor quando comparada com a formulacao que considera o potencial de velocidade no interior do
corpo igual a zero e, portanto, acredita-se que numericamente e mais estavel [8].
Neste trabalho foi utilizada a formulacao que considera o potencial de velocidade no interior do corpo
constante e igual ao potencial de velocidade do escoamento de aproximacao Φ∞. A Eq. 2.6 pode ser
reduzida para uma forma mais simples :
1
4π
ˆ
corpo+esteira
µ∂
∂n
(1
r
)dS − 1
4π
ˆ
corpo
σ
(1
r
)dS = 0 (2.7)
Neste tipo de formulacao, para a Eq. 2.7 ser valida, a intensidade da fonte tem que ser igual a
componente normal a superfıcie do vector velocidade do escoamento de aproximacao Q∞ [5].
σ = n ·Q∞ (2.8)
Assim, ao fixarem-se as intensidades das fontes, as incognitas passam a ser as intensidades dos
dipolos ao longo da superfıcie do corpo e esteira.
2.3 Forma e intensidade da esteira
As hipoteses simplificativas adoptadas sao semelhantes as que sao feitas na teoria da linha sus-
tentadora [lifting line theory ] ou metodo da malha de vortices [vortex lattice method ] na modelacao
de um corpo sustentador. A esteira considera-se plana, indeformavel e paralela ao escoamento de
aproximacao. Com esta simplificacao lineariza-se um problema que e nao linear, ja que a forma da
esteira depende do campo de velocidades induzido pelos vortices arrastados e este campo de veloci-
dades induzido depende da forma da esteira. Esta situacao pode ser resolvida com a introducao de um
processo iterativo onde os pontos que formam os paineis da esteira sao ajustados tendo em conta a
velocidade neles induzida pelos vortices arrastados. Este procedimento e conhecido por relaxacao da
esteira [wake relaxation]. Ao longo das iteracoes a esteira vai-se tornando cada vez mais enrolada junto
as extremidades uma vez que as velocidades induzidas nessas zonas sao provocadas por vortices de
maior intensidade. Trata-se do fenomeno denominado por enrolamento da esteira [wake roll-up]. Em
[19, 20, 21, 5] podem ser consultados calculos de forma iterativa tendo em conta o enrolamento da es-
teira. Uma vez que este processo iterativo requer elevados recursos computacionais - de facto, a cada
iteracao e necessario fazer um novo calculo de escoamento em torno do corpo e esteira, atendendo a
que cada iteracao apresenta-se com formato diferente - optou-se por nao incluı-lo no presente trabalho.
Por outro lado, e no que diz respeito as intensidades dos dipolos que modelam a esteira de vortices
libertada pelo corpo sustentador, as mesmas podem ser escritas em funcao das intensidades presentes
no corpo, atraves da condicao de Kutta aplicada ao longo do bordo de fuga. Esta condicao indica que
o escoamento deve deixar o bordo fuga suavemente. Por outras palavras, garante que as velocidades
9
Figura 2.2: Aplicacao da condicao de Kutta. [5]
tanto para o intradorso como o extradorso, na zona do bordo de fuga, sao iguais. A aplicacao desta
condicao no bordo de fuga e representada pela seguinte expressao,
µU − µL − µW = 0 (2.9)
onde µU e µL correspondem as intensidades dos dipolos no bordo de fuga, para o painel do extradorso
e intradorso, respectivamente.
Por exemplo, a especificacao da condicao de Kutta no bordo de fuga de uma asa, considerando
as intensidades dos dipolos constantes ao longo do painel, pode ser visualizada na Fig. 2.2. Note-se
que um painel a funcionar como um dipolo com intensidade µ e equivalente a este mesmo painel com
quatro vortices de intensidade Γ = µ em cada aresta, conhecido por anel de vortice [vortex ring]. Na
linha que define o bordo de fuga, e considerando que o sentido positivo do dipolo aponta para o interior
do corpo, a aresta do painel superior tem intensidade −ΓU , a aresta do painel inferior +ΓL e para o
painel da esteira +ΓW . O balanco das intensidades dos vortices na linha do bordo de fuga e entao :
− ΓU + ΓL + ΓW = 0 (2.10)
ou, como na Eq. 2.9
ΓW = ΓU − ΓL (2.11)
2.4 Procedimento Numerico
Para a aplicacao numerica, a superfıcie SB e dividida em N paineis no corpo e em Nw na esteira
como mostra a Fig. 2.3. Cada painel que faca parte do corpo comporta-se como um dipolo com
uma intensidade µ e uma fonte com intensidade σ. Diferentemente, um painel que pertenca a esteira
comporta-se so como um dipolo. A condicao de fronteira e imposta no centro de cada painel, nos
chamados pontos de controlo [control points]. A Eq. 2.7 para cada ponto de controlo toma a seguinte
forma :
10
Figura 2.3: Discretizacao de um corpo em paineis. [5]
N∑k=1
1
4π
ˆ
Painel−corpo
µn · ∇(
1
r
)dS −
N∑k=1
1
4π
ˆ
Painel−corpo
σ
(1
r
)dS+
+
Nw∑l=1
1
4π
ˆ
Painel−esteira
µn · ∇(
1
r
)dS = 0. (2.12)
Para cada ponto de controlo sao necessarias as influencias de todos os k paineis do corpo e de
todos os l paineis da esteira. A integracao, ao contrario das seccoes anteriores em que era feita ao
longo de toda a superfıcie do corpo e, agora, feita para cada painel individualmente e, considerando as
intensidades das singularidades ( µ e σ ) unitarias, depende apenas da geometria do painel:
Ck =1
4π
ˆ
Painel−corpo
∂
∂n
(1
r
)dS, (2.13)
Bk = − 1
4π
ˆ
Painel−corpo
(1
r
)dS. (2.14)
As Equacoes 2.13 e 2.14 representam os coeficientes de influencia para o dipolo e fonte, respecti-
vamente, e podem ser determinados se as coordenadas dos vertices do painel k e a distancia r forem
conhecidas. Estes coeficientes traduzem o modo como um determinado painel k influencia um ponto P
a uma distancia r do seu centro, em termos de potencial de velocidade Φ. Por outras palavras, repre-
senta o potencial de velocidade induzido num ponto P a uma distancia r do seu centro considerando a
intensidade da singularidade igual a 1.
Tendo em conta estes coeficientes de influencia, que serao desenvolvidos numa seccao posterior, a
Eq. 2.12 pode ser escrita da seguinte forma:
11
Figura 2.4: Relacao entre o painel da esteira e os dois paineis do bordo de fuga. [5]
N∑k=1
Ckµk +
N∑k=1
Bkσk +
Nw∑l=1
Clµl = 0. (2.15)
A aplicacao da Eq. 2.15 no ponto de controlo de um determinado painel na superficie do corpo,
garante que o somatorio dos potenciais de velocidade induzidos por todos os paineis (dipolos e fontes)
e igual a 0. Esta e a expressao numerica da aplicacao da condicao de fronteira.
As intensidades das fontes σk sao obtidas pela Eq. 2.8 onde n, agora, e a normal ao painel k,
σk = nk ·Q∞, (2.16)
e deste modo, o termoN∑k=1
Bkσk fica conhecido e pode passar para o lado direito da Eq. 2.15.
As intensidades dos dipolos da esteira µl podem ser relacionadas com as intensidades dos dipolos
do corpo localizados no bordo de fuga pela condicao de Kutta, como apresentado na seccao anterior.
Por exemplo, na Fig. 2.4 a intensidade do dipolo da esteira µt pode ser expressa em funcao das duas
intensidades µr e µs dos paineis que se encontram no bordo de fuga,
µt = µr − µs. (2.17)
O potencial de velocidade induzido pelo painel da esteira vem dado por
Ctµt = Ct (µr − µs) . (2.18)
Assim, e eliminada a incognita da intensidade do painel da esteira µt, bastando saber o coeficiente
de influencia do painel Ct. Para a Eq. 2.15 este coeficiente e somado aos coeficientes de influencia
do bordo de fuga quando k = r, e subtraıdo quando k = s. Por exemplo, a aplicacao da condicao de
fronteira, utilizando a Eq. 2.15, no ponto de controlo numero 1 resulta na expressao:
12
C11µ1 + . . .+ (C1r + C1t)µr + (C1s − C1t)µs + . . .+ C1NµN +
N∑k=1
Bkσk = 0, (2.19)
onde por exemplo C1kµk e potencial de velocidade induzido pelo painel k a funcionar como dipolo com
uma intensidade µk no ponto de controlo numero 1.
A Eq. 2.15 pode, assim, ser escrita numa forma ainda mais simples :
N∑k=1
Akµk = −N∑k=1
Bkσk, (2.20)
onde Ak = Ck se k nao pertencer a um painel do bordo de fuga, ou Ak = Ck ± Ct se k pertencer a um
painel do bordo de fuga.
A aplicacao da Eq. 2.20 a todos os pontos de controlo de todos os paineis do corpo, de modo a que
a condicao de fronteira seja satisfeita em toda a superfıcie, resulta num sistema de N equacoes a N
incognitas (µk):
a11 a12 · · · a1N
a21 a22 · · · a2N
· · · · · · · · · · · ·
aN1 aN2 · · · aNN
µ1
µ2
...
µN
=
RHS1
RHS2
...
RHSN
, (2.21)
onde os RHSk sao os chamados termos do lado direito da equacao [Right-Hand Side] e podem ser
determinados uma vez que as intensidades das fontes para cada painel σk sao conhecidas:
RHS1
RHS2
...
RHSN
= −
b11 b12 · · · b1N
b21 b22 · · · b2N
· · · · · · · · · · · ·
bN1 bN2 · · · bNN
σ1
σ2
...
σN
. (2.22)
As matrizes aij e bij contem as influencias de todos os paineis, a funcionarem como dipolos e fontes,
respectivamente, em todos os pontos de controlo. Sao denominadas por matrizes dos coeficientes de
influencia [influence coefficient matrices].
As incognitas (µk) sao obtidas pela solucao do sistema de equacoes lineares. Com esta solucao,
juntamente com as intensidades das fontes, garante-se que no interior do corpo o potencial de veloci-
dade total e igual ao potencial de velocidade do escoamento de aproximacao.
2.4.1 Discretizacao da superfıcie
No presente trabalho, optou-se pela formulacao mais basica relativamente a discretizacao da geometria
e especificacao da forma como as intensidades das singularidades variam em cada painel. Considerou-
se que as distribuicoes continuas das intensidades do dipolo e fonte ao longo da superfıcie do corpo
podem ser aproximadas por uma distribuicao de varios dipolos e fontes com intensidades constantes
por painel. Ou seja, cada painel e considerado plano (sem curvatura) e as intensidades do seu dipolo e
13
Figura 2.5: Construcao do painel plano
fonte constantes. Metodos que utilizem este tipo de formulacao sao conhecidos como metodos de baixa
ordem [low-order methods]. Uma formulacao mais complexa, e que requer mais poder computacional,
passa por admitir curvatura no painel e variacao das intensidades das singularidades ao longo do
painel. Metodos que utilizem esta formulacao sao conhecidos como metodos de elevada ordem [High-
order methods].
Em [11], os autores fazem uma comparacao de resultados obtidos em cinco programas que utilizam
o metodo dos paineis 3D, tanto de baixa como elevada ordem. Os cinco programas utilizados no estudo
apresentam bons resultados quando comparados com dados experimentais. Embora a simulacao de
ordem elevada se tenha revelado superior em situacoes em que o numero de paineis era o mesmo,
tal apenas foi conseguido a custa de um elevado tempo e custo computacional. Na verdade, para o
mesmo tempo de computacao, o metodo de baixa ordem demonstrou ser superior.
Tendo em conta o objectivo deste trabalho, o aumento da precisao da solucao nao compensaria o
tempo de espera da mesma. Tambem, o poder computacional disponıvel por qualquer utilizador, hoje
em dia, ja permite a construcao de uma malha relativamente refinada ao ponto de se obter solucoes
bastante satisfatorias com uma formulacao de ordem inferior.
Apos a discretizacao da superfıcie do corpo em multiplos pontos, torna-se necessario construir os
paineis de modo a que se aproximem o mais possıvel da superfıcie inicial, uma vez que, raramente,
os quatro pontos para a construcao do painel fazem parte do mesmo plano. Assim, e necessario
construir um painel que passe o mais proximo desses quatro pontos mantendo-se plano como mostra,
por exemplo, a Fig. 2.5. O aparecimento de fugas [leakage] entre os paineis sao inevitaveis e ha
que reduzi-las ao maximo com o metodo escolhido para a construcao do painel. Estas fugas sao
apresentadas na Fig. 2.6.
A metodologia utilizada para a construcao do painel plano, inicialmente seguia a linha descrita em
[4] mas, por apresentar maior simplicidade, optou-se pela descrita em [8]. A maior diferenca entre estas
duas reside na forma de calculo do ponto de controlo e do referencial local ao painel. Enquanto que
no primeiro o ponto de controlo e calculado no centroide do painel, no ultimo e obtido pela media dos
quatro vertices do painel.
14
Figura 2.6: Aparecimento de fugas na construcao dos paineis. [5]
Figura 2.7: Painel a funcionar como a) fonte b) dipolo. [5]
Os integrais dos coeficientes de influencia nas equacoes 2.13 e 2.14, considerando painel plano
e intensidade da singularidade constante (Fig. 2.7), sao desenvolvidos em funcao dos pontos das
extremidades do painel por Hess e Smith em [22] ou apresentados numa forma mais simples em [5]: o
potencial de velocidade induzido num ponto qualquer P por um painel a trabalhar como fonte e assim
dado por :
Φf =−σ4π{[aux112 + aux123 + aux134 + aux141] + | z | [aux212 + aux223 + aux234 + aux241]} (2.23)
enquanto que para o dipolo :
Φd =µ
4π[aux212 + aux223 + aux234 + aux241] (2.24)
com,
aux1ij =(x− xi) (yj − yi)− (y − yi) (xj − xi)
dijln
(ri + rj + dijri + rj − dij
), (2.25)
15
aux2ij = arctan
(mijei − hi
zri
)− arctan
(mijej − hj
zrj
), (2.26)
onde,
dij =
√(xj − xi)2
+ (yj − yi)2, (2.27)
mij =yj − yixj − xi
, (2.28)
e,
rk =
√(x− xk)
2+ (y − yk)
2+ z2, (2.29)
ek = (x− xk)2
+ z2, (2.30)
hk = (x− xk) (y − yk) , (2.31)
em que k = 1, 2, 3 e 4.
As equacoes 2.23 e 2.24 sao utilizadas para o calculo das matrizes dos coeficientes de influencia
aij e bij , contidas nas equacoes 2.21 e 2.22, onde o ponto P e substituido pelo ponto de controlo i
no referencial local ao painel j. As construcoes destas matrizes fazem parte da fase mais pesada em
termos computacionais do metodo dos paineis 3D.
2.5 Calculo das Forcas Aerodinamicas
Uma vez conhecidas as intensidades dos dipolos (µk), as velocidades nos referenciais locais aos
paineis (l,m, n) podem ser obtidas por diferenciacao.
As velocidades tangentes a superfıcie sao dadas por:
ql =∂µ
∂l, (2.32)
qm =∂µ
∂m, (2.33)
e a componente normal,
qn = −σ. (2.34)
Para as velocidades nas direccoes tangenciais, os resultados podem variar consoante o modelo
numerico escolhido para a diferenciacao. Normalmente e escolhido o metodo das diferencas finitas
16
Figura 2.8: Esquema de paineis para o calculo da velocidade no painel k. [8]
Figura 2.9: Esquema de paineis numa extremidade para o calculo da velocidade no painel k. [8]
centrais. Por exemplo, para direccao l, temos
ql =µk+1 − µk−1
∆l, (2.35)
onde ∆l e a distancia entre os pontos de controlo dos paineis k + 1 e k − 1.
Neste estudo, o metodo de obtencao da velocidade e semelhante ao descrito em [23], onde e feita
uma interpolacao das intensidades dos dipolos numa curva de segunda ordem e a obtencao do declive
no ponto desejado. Tal e possıvel com a funcao polyfit do MATLAB© que encontra os tres coeficientes
da funcao polinomial de segunda ordem. Se os dados de entrada para a funcao forem cuidadosamente
inseridos de modo a que o referencial seja considerado no painel central (k), o declive neste painel
(velocidade do escoamento no referencial local ao painel k) e igual ao segundo coeficiente obtido da
funcao polyfit. Isto e feito na direccao l e m caso o painel k tenha um vizinho imediatamente a frente e
atras nestas duas direccoes, como mostra a Fig. 2.8.
17
Nas situacoes em que seja necessario obter a velocidade do escoamento nas extremidades da asa
ou no bordo de fuga, o facto do painel K so ter um painel vizinho numa das direccoes implica que se
tenha que arranjar um terceiro painel imediatamente depois desse vizinho, conforme descrito em [8] e
ilustrado na Fig. 2.9.
A velocidade total no painel k e entao a soma da velocidade local, obtida pela diferenciacao das
intensidades dos dipolos, com a velocidade do escoamento de aproximacao,
Qk = (Q∞l, Q∞m
, Q∞n)k + (ql, qm, qn)k. (2.36)
Repare-se que a componente normal da velocidade ao painel Q∞n+ qn e igual a zero.
Depois de obtidas as velocidades, e como se trata de um escoamento invıscido, a passagem para
coeficiente de pressao e imediata,
Cpk = 1− Q2k
Q2∞
(2.37)
onde Q2k =
(Q2kl
+Q2km
+Q2kn
)e Q2
∞ =(Q2∞l
+Q2∞m
+Q2∞n
).
A contribuicao de cada painel para a forca aerodinamica e dada por,
∆Fk = Cpk
(1
2ρQ2∞
)∆Sknk. (2.38)
onde ρ corresponde a densidade do fluido e ∆Sk a area do painel k.
A forca resultante do escoamento invıscido em torno do corpo e alcancada somando os ∆Fk para
todos os paineis k (k = 1→ N ). Esta forca resultante pode ser decomposta na direccao perpendicular
ao escoamento de aproximacao dando origem a um termo conhecido como sustentacao, L. Na direccao
do escoamento, tambem ha uma componente da forca resultante conhecida por resistencia induzida,
Dind. Esta resistencia, que em 2D e igual a zero [24, 18, 5], esta presente pelo facto de se tratar de
um corpo tridimensional. Os vortices libertados pelas extremidades do sistema sustentador induzem
velocidades descendentes (w) no escoamento de aproximacao (Q∞), e portanto, cada seccao da asa
fica a trabalhar num escoamento com uma velocidade efectiva igual a soma vectorial de Q∞ com w,
resultando numa alteracao do angulo de ataque de uma quantidade,
αi = arctan
(w
Q∞
)(2.39)
designada por angulo de ataque induzido [induced angle of attack ]. Cada seccao, com angulo de
ataque geometrico α relativamente ao escoamento de aproximacao Q∞ fica, na realidade, a operar a
um angulo de ataque inferior. O angulo de ataque efectivo a que determinada seccao transversal opera
e definido como,
αef = α− αi (2.40)
e encontra-se representado na Fig. 2.10. E de notar que o angulo de ataque induzido αi e tanto maior
18
Figura 2.10: Esquema do aparecimento da resistencia induzida para uma seccao transversal de umaasa tri-dimensional. [25]
quanto mais proximo estiver das extremidades do corpo sustentador.
Neste momento, considerando escoamento invıscido, os coeficientes de sustentacao e de resistencia
podem ser obtidos por:
CL =L
(0.5ρQ2∞)S
, (2.41)
CD =Dind
(0.5ρQ2∞)S
(2.42)
onde S e igual a uma area de referencia, normalmente utiliza-se a area da asa com vista de topo.
19
Capıtulo 3
Efeito Solo
Um corpo sustentador a operar junto ao solo apresenta uma alteracao nas caracterısticas aerodinamicas
quando comparado com uma situacao em que o solo nao esteja presente. Na aproximacao do corpo
ao solo, as caracterısticas aerodinamicas sao alteradas devido a dois fenomenos: o efeito solo e o
efeito Venturi. Enquanto que no primeiro ha um aumento da sustentacao e diminuicao da resistencia,
no segundo a sustentacao diminui e a resistencia aumenta significativamente [26, 27].
No efeito solo o aumento da sustentacao e diminuicao da resistencia esta relacionado com a
alteracao da distancia horizontal entre os vortices libertados pelas extremidades da asa, o que que
leva a uma diminuicao da sua influencia no escoamento de aproximacao [27]. Traduzindo-se numa
menor gama de angulos induzidos nas seccoes transversais ao longo da envergadura da asa.
Ja no efeito Venturi, a proximidade do intradorso da asa com o solo e semelhante a uma passagem
convergente-divergente de uma tubeira, o que conduz a que as velocidades na zona convergente se-
jam bastante elevadas, provocando uma forca descendente superior a que se verifica quando o solo
nao esta presente. Este efeito e bastante explorado em desportos motorizados como por exemplo na
Formula 1 [26].
Uma vez que no presente trabalho a construcao da ferramenta computacional e feita com a utilizacao
do metodo dos paineis, a simulacao do efeito de uma superfıcie solida junto ao corpo passa pela
utilizacao do metodo das imagens. E atraves deste metodo que, em fluido perfeito, se modela qualquer
escoamento que se processe na presenca de uma fronteira solida, em que o escoamento de um dos
lados da fronteira pode ser sempre interpretado como a imagem, nessa fronteira, do escoamento real
ocorrendo do outro lado [18]. Com auxılio da Fig. 3.1 e possıvel concluir que os vortices libertados
pelas extremidades da imagem vao induzir ao escoamento de aproximacao uma velocidade ascendente
contrariando a velocidade descendente induzida pelos vortices da asa, resultando num angulo induzido
inferior comparativamente a um caso em que o solo nao esteja presente.
20
Figura 3.1: Esquema da velocidade induzida pelos vortices das extremidades da imagem da asa
3.1 Aplicacao do metodo das imagens
Esta seccao sera dedicada a implementacao do metodo das imagens no codigo do metodo dos paineis
3D.
Atraves do metodo das imagens procedeu-se a simulacao de um escoamento potencial em torno de
um corpo a operar junto a uma superfıcie solida sendo que, para a sua implementacao, foi necessario
comecar por construir a imagem do respectivo corpo. A imagem e de construcao simples na medida em
que a diferenca para corpo verifica-se apenas no sinal da componente perpendicular ao solo. O metodo
dos paineis 3D descrito no capıtulo anterior foi, nesta fase, aplicado ao corpo e imagem inseridos no
mesmo escoamento e, portanto, a interaccao entre estes teve que ser contabilizada.
Na Fig. 3.2 apresenta-se uma asa com a sua imagem para a simulacao de uma superficie na
proximidade de uma das extremidades, onde tambem foi incluıda a sequencia, adoptada neste trabalho,
de aplicacao da condicao de fronteira (Eq. 2.20) para cada ponto de controlo ao longo de todo o corpo
e imagem, resultando na Eq. 2.21. Com a presenca da imagem, as matrizes dos coeficientes de
influencia aij e bij passam a ter uma dimensao de 2N × 2N , ja que cada ponto de controlo vai ser
induzido, em termos de potencial de velocidade, por todos os paineis do corpo e imagem. Porem,
repare-se que, em qualquer ponto de controlo na superfıcie do corpo, o potencial de velocidade induzido
por todos os paineis da imagem e o mesmo que numa situacao inversa, ou seja, a imagem deste ponto
de controlo induzido por todos os paineis do corpo. Por exemplo, ainda na Fig. 3.2, o potencial de
velocidade induzido no painel numero 2 e igual ao potencial induzido no painel numero N + 2. Tendo
em conta este factor, as construcoes das matrizes dos coeficientes de influencia aij e bij podem ser
executadas de forma relativamente simples e sem apresentarem um custo computacional acrescido,
quando comparado com uma situacao em que o solo nao esteja presente, inerente a aplicacao das
2N equacoes que satisfacam a condicao de fronteira ao longo da superfıcie do corpo e imagem. A
matriz dos coeficientes de influencia aij , que contem a influencia de todos os paineis a comportarem-
se como dipolos em todos os pontos de controlo, e apresentada da seguinte forma e com a respectiva
simplificacao,
aij =
aasa−asa aasa−imagem
aimagem−asa aimagem−imagem
=
aasa−asa aasa−imagem
aasa−imagem aasa−asa
(3.1)
21
Figura 3.2: Asa e imagem, com a representacao da sequencia da aplicacao da condicao de fronteira
onde aasa−asa e aimagem−asa sao as matrizes dos coeficientes de influencia de todos os paineis da
asa em todos os pontos de controlo da asa e imagem, respectivamente. Por sua vez, aasa−imagem
e aimagem−imagem sao as matrizes dos coeficientes de influencia de todos os paineis da imagem em
todos os pontos de controlo da asa e imagem, respectivamente. Cada uma destas matrizes tem uma
dimensao N ×N .
A construcao da matriz dos coeficientes de influencia bij , em que neste caso os paineis comportam-
se como fontes, foi feita de forma analoga.
Na capıtulo de Validacao e Discussao dos Resultados, e com objectivo de demonstrar os efeitos nas
caracterısticas aerodinamicas de uma asa com a presenca de uma superfıcie plana, serao comparados
resultados obtidos pelo programa comercial STAR-CCM+ com os calculados pelo metodo dos paineis
3D.
3.2 Camada Limite Atmosferica
Uma vez que o objectivo do presente trabalho passa pela previsao das caracterısticas aerodinamicas
de uma vela rıgida de um catamaran, seria interessante que a simulacao fosse o mais proximo de uma
situacao real. Assim como o efeito solo foi introduzido para simular a superfıcie da agua junto a uma
22
Descricao da superfıcie z0 [mm]
Mar calmo 0.20Mar agitado 0.50
Tabela 3.1: Rugosidade da superfıcie
das extremidades da vela, seria interessante simular tambem o efeito da camada limite atmosferica. Ou
por outras palavras, a variacao da velocidade do escoamento de aproximacao.
Uma forma de aproximar o perfil de velocidades da camada limite atmosferica passa pela utilizacao
de uma funcao logarıtmica [28]:
u(z) = u(zr)ln( zz0 )
ln( zrz0 ), (3.2)
onde zr e uma altura de referencia, u(zr) velocidade de referencia, z0 a rugosidade da superfıcie e
z a altura acima da linha media da superfıcie da agua. A rugosidade da superfıcie z0 varia com as
condicoes do mar, e os valores tıpicos estao na Tab. 3.1. A altura de referencia zr para uma velocidade
de referencia u(zr) e tipicamente 10 metros [28].
A introducao do efeito da camada limite atmosferica no calculo das caracterısticas aerodinamicas
pelo metodo dos paineis 3D, e feita quando aplicada a Eq. 2.16 na obtencao das intensidades das
fontes, e na Eq. 2.36 no calculo da velocidade do escoamento num determinado painel.
Para cada painel k, e calculada a altura z a que se encontra do solo e a velocidade u(z) correspon-
dente a essa altura na camada limite atmosferica, Eq. 3.2, com a velocidade de referencia u(zr) igual
a velocidade do escoamento de aproximacao Q∞. Portanto, para ter em conta os efeitos da camada li-
mite atmosferica a velocidade do escoamento de aproximacao Q∞, utilizada nas equacoes 2.16 e 2.36,
e substituıda por u(z) onde varia de painel para painel consoante a sua altura ao solo. Os paineis mais
proximos do solo estarao sujeitos a uma velocidade inferior a que se faz sentir a alturas superiores. Com
efeito, as caracterısticas aerodinamicas da asa no interior da camada limite atmosferica sao alteradas
uma vez que a solucao da Eq. 2.21, agora, tem em conta a variacao da velocidade do escoamento de
aproximacao. Na Fig. 3.3 esta representado esquematicamente o que foi descrito anteriormente.
23
Figura 3.3: Esquema de uma superfıcie sustentadora a operar proxima do solo com uma camada limiteatmosferica
24
Capıtulo 4
Acoplamento dos Efeitos Viscosos
Perante um caso de escoamento de fluido real em torno de um corpo solido, assiste-se a presenca de
duas regioes onde o fluido assume caracterısticas diferentes. A primeira, formada em contacto com a
superfıcie onde os gradientes de velocidade sao bastante elevados e os efeitos viscosos nao podem
ser desprezados. Por outro lado, na segunda, os gradientes de velocidades sao pequenos o suficiente
para que os efeitos viscosos sejam desprezados aproximando-se portanto a um escoamento de fluido
ideal. Ate a presente seccao assumiu-se que a primeira regiao era de pequena espessura podendo
entao ser desprezada. Na Fig. 4.1 estao representadas estas duas regioes.
A primeira regiao, onde os efeitos viscosos nao podem ser desprezados, e conhecida como camada
limite. No seu interior, o perfil de velocidade varia de zero na superficie, devido a condicao de nao
escorregamento, ate a uma velocidade proxima da velocidade exterior (Fig. 4.2 a)). Esta variacao da
velocidade esta relacionada com as tensoes de corte que se fazem sentir entre os elementos de fluido.
A espessura da camada limite δ e definida pela distancia desde a superfıcie do corpo ate a uma posicao
onde a velocidade seja igual a 99% da velocidade exterior.
Considerando a Fig. 4.2, na qual esta representada uma distribuicao de velocidade, na proximidade
de uma superfıcie, para os escoamentos real e invıscido, e possıvel retirar, sobretudo, duas conclusoes.
Por um lado, no escoamento real, a velocidade varia continuamente de 0 U(y = 0) = 0 ate uma
velocidade proxima da velocidade exterior U(y = δ) ≈ Ue, enquanto que no invıscido a velocidade
mantem-se constante. Por outro lado, ao considerar-se o fluido real, verifica-se um deficit de caudal
associado a presenca da camada limite. Sendo Udy o caudal volumetrico atraves de uma seccao
Figura 4.1: Decomposicao do escoamento em duas regioes: Invıscida e Viscosa
25
Figura 4.2: Perfil de velocidade para um escoamento a) real e b) invıscido.
Figura 4.3: Espessura de deslocamento
dy×1, a diferenca de caudais volumetricos escoados numa seccao δ×1 para um escoamento de fluido
perfeito e real e dada por :´ δ0Uedy −
´ δ0Udy.
Este deficit de caudal deve-se ao facto de se estar a considerar fluido real. Em fluido perfeito, e
possivel simular o efeito da camada limite ao garantir-se que o caudal escoado e o mesmo que num
caso de fluido real. Para o efeito, e necessario deslocar a superficie do corpo numa distancia δ∗, como
mostra a Fig. 4.3, de modo a que o caudal escoado a uma velocidade Ue em condicao de fluido
perfeito seja o mesmo que num caso de escoamento real. Esta distancia, designada por espessura de
deslocamenteo [displacement thickness], e escrita do seguinte modo:
δ∗Ue =
ˆ δ
0
(Ue − U) dy ⇔ δ∗ =
ˆ δ
0
(1− U
Ue
)dy (4.1)
Este parametro tambem pode ser interpretado como o deslocamento que as linhas de corrente sofrem
de forma a que a conservacao de massa seja satisfeita ao longo de um escoamento de fluido real[25].
De modo semelhante ao que foi efectuado anteriormente para o deficit de caudal volumetrico,
tambem e possıvel para a quantidade de movimento. Portanto, o parametro ligado ao deficit de quan-
tidade de movimento, associado ao facto de se considerar fluido real, e chamado de espessura de
quantidade de movimento [momentum thickness] e e igual a :
26
θ =
ˆ δ
0
U
Ue
(1− U
Ue
)dy (4.2)
4.1 Introducao
Considera-se que os efeitos viscosos manifestam-se so no interior da camada limite. Esta camada
desenvolve-se na superficie do corpo desde o ponto de estagnacao prolongando-se para jusante sob
a forma de esteira. Em todo o restante campo o escoamento comporta-se como fluido perfeito e a
solucao pode ser obtida pelo metodo dos paineis.
As equacoes da camada limite 2D sao dadas por,
U∂U∂x + V ∂U
∂x = − 1ρdpedx + ν ∂
2U∂y2
∂U∂x + ∂V
∂y = 0
(4.3)
sistema de duas equacoes a duas incognias, onde U e V sao as componentes horizontal e vertical,
respectivamente, da velocidade do escoamento. Se se pretender obter informacao apenas quanto aos
efeitos globais, ao longo da superfıcie do corpo, produzidos por um escoamento de fluido real, e mais
economica a utilizacao da equacao integral de von-Karman,
dθ
dx+ θ
H + 2
Ue
dUedx
=Cf2
(4.4)
onde θ corresponde a espessura de quantidade de movimento, H = δ∗/θ o factor de forma [form factor ],
Ue velocidade exterior e Cf = τw/0.5ρU2e o coeficiente de tensao de corte superficial [skin friction coeffi-
cient ]. O desenvolvimento da Eq. 4.3 ate a Eq. 4.4 pode ser encontrado, por exemplo, em [18, 25]. De
realcar que a Eq. 4.4 e o ponto de partida para os metodos simples como de Thwaites e Heads para
escoamentos laminares e turbulentos, respectivamente. Estes metodos hoje em dia ainda sao muito
utilizados, ja que apenas precisam das distribuicoes de velocidade ao longo da superficie do corpo,
enquanto que a Eq. 4.3 para alem da distribuicao anterior ainda precisa da distribuicao de velocidade
na direccao prependicular a superfıcie.
A resolucao da Eq. 4.4 pressupoe que se conheca a distribuicao de velocidades exteriores Ue pelo
que, antes de mais, e necessario que se faca uma analise invıscida para a obtencao desta distribuicao.
O calculo da distribuicao de velocidades ira, posteriormente, permitir avaliar o comportamento da ca-
mada limite ao longo da superfıcie do corpo. Uma vez que a camada limite interage com o escoa-
mento exterior - atendendo a que, devido a condicao de nao escorregamento, as linhas de corrente sao
deslocadas em uma distancia δ∗- a distribuicao de velocidades do escoamento invıscido sofrera uma
alteracao relativamente a primeira solucao obtida. Deste modo, uma vez conhecida a distribuicao de
δ∗ ao longo da superfıcie do corpo e necessario fazer um novo calculo de fluido perfeito mas, agora,
em torno de um corpo fictıcio com a forma do corpo original mais a espessura do deslocamento (corpo
solido + δ∗). Depois de calculada esta nova distribuicao de velocidades Ue(x) na superfıcie do corpo
fictıcio procede-se ao calculo de um novo comportamento da camada limite de onde se obtem um novo
27
Figura 4.4: Alteracao nas linhas de corrente, com o objectivo de simular os efeitos viscosos, atraves deuma velocidade de transpiracao. [29]
δ∗ [18]. Os passos anteriores sao repetidos e estamos perante um processo iterativo. Repare-se que
a medida que a espessura de deslocamento aumenta o corpo fictıcio fica a operar a um menor angulo
de ataque efectivo do que o corpo solido original [17, 18].
Nao e computacionalmente eficiente alterar a geometria do corpo sempre que se pretenda fazer uma
nova analise invıscida depois de calculada a espessura de deslocamento. E mais economico utilizar,
a cada iteracao, a mesma geometria do corpo inicial e simular o efeito do deslocamento das linhas
de corrente com uma velocidade de transpiracao [transpiration velocity ]. Esta velocidade pode ser
simulada alterando a condicao de fronteira na superfıcie do corpo, passando de superfıcie impermeavel
(Vw = 0) a uma superfıcie com sopro (Vw > 0), como ilustrado para um perfil 2D na Fig. 4.4, e com uma
intensidade [18, 5]:
Vw =d
dx(Ueδ
∗) (4.5)
O processo aqui descrito revela, num modo geral, o modo como se processa o acoplamento dos
efeitos viscosos a um escoamento invıscido. Porem, a solucao nem sempre converge, como sucede nos
casos de escoamentos que apresentam separacao da camada limite. Com efeito, e para a obtencao
da solucao desejada, foi necessario desenvolver algumas tecnicas de acoplamento. De seguida sao
apresentados quatro esquemas basicos de interaccao entre as regioes viscosa e invıscida.
4.2 Esquemas de acoplamento
As proximas seccoes serao dedicadas a apresentacao de quatro esquemas - os mais basicos -de
interaccao viscosa-invıscida, atraves dos quais e possıvel combinar, de forma iterativa, as solucoes
das equacoes, obtidas separadamente, que caracterizam o escoamento no interior e exterior a camada
limite de forma a proporcionar uma solucao global do escoamento.
Matematicamente a interaccao entre a regiao invıscida e viscosa pode ser escrita da seguinte forma:
28
Figura 4.5: Metodo directo.[32]
~ue = E [δ∗]
~ue = B [δ∗]
(4.6)
representando um sistema nao linear [30]. Na Eq. 4.6 ~ue e o vector velocidade no topo da camada
limite, δ∗espessura de deslocamento, E e B o conjunto das equacoes para o escoamento externo e
interno a camada limite, respectivamente.
4.2.1 Metodo Directo
O esquema de interaccao mais simples e intuitivo e o chamado metodo directo. A solucao do escoa-
mento invıscido e dada pela velocidade ao longo da superficie da asa. Esta velocidade e considerada
a velocidade exterior a camada limite ~ue e portando, as equacoes que definem o escoamento interior
a camada limite usam essa velocidade para o calculo da espessura de deslocamento δ∗, entre ou-
tros parametros da camada limite. Posteriormente, esta espessura de deslocamento e utilizada para a
construcao de uma nova geometria (corpo + δ∗), ou, mais facilmente, simulada com uma velocidade de
transpiracao. O metodo continua por iteracao entre o codigo invıscido e o da camada limite [29, 31].
Este e o metodo mais simples porque a unica alteracao a ser feita e, a cada iteracao, a alteracao da ge-
ometria do corpo ou a implementacao da velocidade de transpiracao no metodo invıscido (neste caso,
metodo dos paineis). No entanto, o metodo directo apresenta grandes limitacoes. Para escoamentos
com uma grande interaccao entre a parte inviscida e viscosa, como por exemplo em caso de separacao,
os resultados nao convergem [14]. A expressao que descreve este metodo e a seguinte:
~u(n)e = E
[δ∗(n−1)
]δ∗(n) = B−1
[~u
(n)e
] (4.7)
onde n e o numero da iteracao. A Fig. 4.5 mostra graficamente o que esta representado na Eq. 4.7.
4.2.2 Metodo Inverso
Com o objectivo de contornar o problema de convergencia no caso de separacao da camada limite no
esquema anterior, as equacoes que definem o escoamento exterior e interior a camada limite serao,
agora, resolvidas de forma inversa [30, 32]. No entanto, com este esquema ainda nao e possivel a
resolucao de escoamentos com separacoes de grandes dimensoes e a convergencia para a solucao
desejada pode ser lenta [14, 32]. A expressao matematica para o metodo inverso e a seguinte:
29
Figura 4.6: Metodo inverso.[32]
δ∗(n) = E−1
[~u
(n−1)e
]~u
(n)e = B
[δ∗(n)
] (4.8)
ou de forma grafica na Fig. 4.6.
O modo de utilizacao deste esquema e a forma de calculo dos parametros da camada limite de
forma inversa podem ser consultados, com maior detalhe, em [25, 33].
4.2.3 Metodo Semi-Inverso
Neste metodo, a solucao do escoamento exterior e calculada de forma semelhante a apresentada no
metodo directo, enquanto que a do escoamento interior continua a ser calculada como apresentada
no metodo inverso, evitando a nao convergencia no caso de haver separacao da camada limite. A
expressao que descreve o presente metodo e a seguinte:
~u
(n)e = E
[δ∗(n−1)
]~u
(n)v = B
[δ∗(n−1)
]δ∗(n) = δ∗(n−1) + ω
(u
(n)v − u(n)
e
) (4.9)
onde ue e a velocidade obtida da formulacao do escoamento exterior, uv a velocidade na extremidade
da camada limite obtida pelas equacoes da regiao viscosa e ω e o parametro de relaxacao. A solucao
converge quando ue = uv.
Este metodo e utilizado por Cabeci e Besnard em [13] para o calculo do escoamento em torno de
uma asa com a presenca de um flap. Os resultados sao bastante satisfatorios para casos bidimensi-
onais, permitindo prever com precisao o angulo para o coeficiente de sustentacao maximo. Ja para
casos tridimensionais, os resultados sao satisfatorios mas a previsao do angulo para o coeficiente de
sustentacao maximo nao e tao precisa como no caso bidimensional.
4.2.4 Metodo Simultaneo
Este metodo baseia-se no facto de nao considerar a hierarquia entre as solucoes do escoamento ex-
terior e interior a camada limite. As solucoes desses dois conjuntos de equacoes, que caracterizam a
zona viscosa e invıscida do escoamento, sao obtidas em simultaneo (Fig. 4.7) pelo metodo iterativo de
Newton.
30
Figura 4.7: Metodo simultaneo.[32]
De entre os esquemas ja apresentados, este e o mais robusto e o que apresenta melhores resul-
tados em termos de convergencia. Porem, apresenta tambem um maior custo computacional devido
a inversao da matriz Jacobiana, necessaria para a linearizacao do conjunto de equacoes nao lineares
[29, 30, 12].
Este e o esquema utilizado no software XFOIL [10] para o acoplamento dos efeitos viscosos a
solucao invıscida.
4.3 Interaccao com o XFOIL
O XFOIL e um programa que esta disponıvel de forma gratuita em [34], desenvolvido por Mark Drela e
que permite calcular as caracteristicas aerodinamicas de um perfil bidimensional tanto para um escoa-
mento invıscido como viscoso.
Numa analise invıscida, o calculo e realizado com recurso ao metodo dos paineis. Para tal, a
superficie do perfil e discretizada em paineis planos com singularidades do tipo fonte e vortice. A
singularidade do tipo fonte tem uma intensidade constante por painel enquanto que a do tipo dipolo
varia de forma linear.
No caso de uma analise viscosa, Drela basiou-se no seu trabalho em [35] para o acoplamento
das equacoes integrais da camada limite bidimensional com a formulacao invıscida do metodo dos
paineis. Este conjunto de equacoes que definem as duas regioes, viscosa e inviscida, sao resolvidas
em simultaneo pelo metodo de Newton. No final do calculo, para alem dos esforcos a que o perfil fica
sujeito, ha um conjunto de parametros ligados a camada limite que ficam conhecidos (i.g. espessura da
camada limite δ, espessura de deslocamento δ∗ e factor de forma H). Uma descricao mais detalhada
do modo de funcionamento do XFOIL pode ser encontrada em [10, 36].
O XFOIL e utilizado neste trabalho para o calculo do parametro da camada limite δ∗, importante
para o acoplamento dos efeitos viscosos a solucao invıscida atraves da Eq. 4.5.
Uma vez que o XFOIL e uma ferramenta de analise de perfis bidimensionais, e estando perante
um escoamento tridimensional, assumiu-se que a camada limite desenvolvida na superficie do corpo
pode ser aproximada por varias camadas limites bidimensionais ao longo da envergadura da asa, como
nos trabalhos [37, 38, 17, 9, 8]. A asa e dividida em seccoes bidimensionais, como mostra a Fig. 4.8,
as quais, posteriormente, serao analisadas de forma independente pelo XFOIL. Com recurso a uma
funcao construida por Gus Brow [39], torna-se possıvel e simples a interaccao entre MATLAB e XFOIL.
31
Figura 4.8: As duas primeiras seccoes para a analise viscosa
4.3.1 Asa
No caso em que so esteja presente uma asa no seio do escoamento, os dados necessarios para uma
analise viscosa pelo XFOIL, para cada seccao, sao os seguintes:
• Perfil
• Numero de Reynolds
• Angulo de ataque efectivo (αef )
O perfil (representado a vermelho na Fig. 4.8) e construido com o conjunto de coordenadas dos pontos
de controlo da respectiva seccao. O numero de Reynolds e calculado em relacao a corda do perfil. E por
fim, o angulo de ataque efectivo e obtido atraves da soma do vector da velocidade do escoamento de
aproximacao com o vector da velocidade induzida, na respectiva seccao, pelos vortices que modelam
a esteira. Uma vez que a analise invıscida ja foi efectuada, as intensidades µ dos paineis da esteira
sao conhecidas e portanto a velocidade induzida por estes pode ser obtida atraves da derivada nas
tres direccoes da Eq. 2.24 (Velocidade induzida =(∂Φd
∂x ,∂Φd
∂y ,∂Φd
∂z
)). O resultado desta derivada, que
corresponde a velocidade induzida por um painel a comportar-se como um dipolo, pode ser consultada
em [5].
Depois de analisada a uma determinada seccao, um dos outputs ligados a caracterizacao da ca-
mada limite devolvido pelo XFOIL - o parametro δ∗-, e utilizado para o calculo da velocidade de transpiracao
pela Eq. 4.5.
Na presenca de malhas pouco refinadas na direccao da corda, os pontos de controlo para a construcao
do perfil podem ser de tal forma escassos que nao permitam ao XFOIL a convergencia dos seus re-
sultados. Neste caso, para a seccao em questao, os efeitos viscosos nao sao acoplados a solucao
invıscida.
Outro parametro da camada limite que e obtido pela analise a cada seccao da asa e o coeficiente
de resistencia devido a tensao de corte na superfıcie do corpo CDτw . Este parametro e utilizado para o
32
Figura 4.9: Seccao do Flap para analise no XFOIL
calculo da resistencia total.
4.3.2 Asa com Flap
Numa analise de um escoamento em torno de dois elementos, conjunto asa e flap, a metodologia
utilizada para a analise da seccao e ligeiramente alterada. A analise e realizada separadamente para
cada elemento. Na analise da asa, o metodo e identico ao explicado na seccao anterior. No que diz
respeito a analise no flap, e tendo em conta a dificuldade na previsao do angulo de ataque efectivo a
que determinada seccao deste elemento fica sujeito, a funcao de Gus Brow foi alterada de modo a que
o input para o XFOIL seja o seguinte:
• Perfil
• Numero de Reynolds
• Coeficiente de sustentacao da seccao (Cl)
Neste caso, o XFOIL devolve os parametros de uma camada limite desenvolvida no perfil a um de-
terminado numero de Reynolds com um angulo de ataque de modo a que produza um coeficiente de
sustentacao Cl (Fig. 4.9). Ou seja, o XFOIL “procura” um angulo de ataque que satisfaca o Cl de
entrada.
Repare-se que tanto o angulo efectivo αef como o coeficiente de sustentacao da seccao do flap
Cl sao obtidos a partir da solucao do escoamento invıscido pelo que todo o processo iterativo, para o
calculo dos parametros da camada limite, fica a cargo do XFOIL.
No proximo capıtulo procurar-se-a explicar a forma atraves da qual a espessuras de deslocamento,
obtidas para cada seccao, sao acopladas ao modelo invıscido atraves de uma velocidade de transpiracao
na superfıcie do corpo.
33
4.4 Procedimento Numerico
Uma forma precisa para a representacao do escoamento total, que tem em conta a influencia da ca-
mada limite e a solucao invıscida, e conseguida com um termo adicional nas intensidades das fontes
[12, 37, 38].
Assim, a formulacao apresentada no capıtulo Metodo dos Paineis 3D e alterada para contemplar a
influencia deste termo. A Eq. 2.7 vem agora com σ = σinv + σvis,
1
4π
ˆ
corpo+esteira
µ∂
∂n
(1
r
)dS − 1
4π
ˆ
corpo
σinv(
1
r
)dS − 1
4π
ˆ
corpo+esteira
σvis(
1
r
)dS = 0 (4.10)
onde σvis e o novo termo da fonte para adicionar a espessura de deslocamento. Repare-se que este
termo e adicionado na superfıcie do corpo e na esteira. Aplicou-se condicao de fronteira (Eq. 4.10) a
todos os pontos de controlo da superfıcie do corpo, como feito na seccao Procedimento Numerico do
capıtulo Metodo dos Paineis 3D, resultando na seguinte equacao,
aµ = −bσinv − bvisσvis (4.11)
onde a e b sao as matrizes dos coeficientes de influencia para os paineis a funcionarem como dipolos
e fontes, respectivamente; bvis e a matriz dos coeficientes de influencia dos paineis do corpo e da
esteira nos pontos de controlo do corpo. Assim, a matriz bvis e o vector σvis apresentam uma dimensao
(N × (N +Nw)) e (N +Nw)× 1, respectivamente.
A solucao do escoamento total pode ser considerada como um resultado da soma entre as solucoes
invıscida e viscosa,
µ = µinv + µvis. (4.12)
com,
µinv = a−1(−bσinv) (4.13)
e,
µvis = a−1(−bvisσvis) (4.14)
onde o termo µinv e obtido pela solucao da Eq. 2.21.
Por exemplo, a expressao geral para a intensidade do dipolo no ponto de controlo do painel i e a
seguinte:
µi =
N∑j=1
a−1ij
(−
N∑k=1
bjkσinvk
)+
N+Nw∑l=1
N∑k=1
(a−1ik
(−bviskl
))σvisl (4.15)
34
Figura 4.10: Deslocamento da linha de corrente atraves de uma velocidade normal a superficie
onde o primeiro termo do segundo membro corresponde a solucao invıscida,
µinvi =
N∑j=1
a−1ij
(−
N∑k=1
bjkσinvk
)(4.16)
e o segundo termo corresponde a correccao viscosa,
µvisi =
N+Nw∑l=1
N∑k=1
(a−1ik
(−bviskl
))σvisl . (4.17)
A Eq. 4.15 pode ser escrita da seguinte forma:
µi = µinvi +
N+Nw∑l=1
Hilσvisl (4.18)
A solucao geral vem unicamente em funcao do novo termo adicionado as fontes σvis e pode ser
escrita da seguinte foma:
µ = µinv +Hσvis (4.19)
com H = a−1(−bvis). A unica incognita na Eq. 4.19 e o vector coluna σvis. Este vector contem
as intensidades de cada painel necessarias para o deslocamento das linhas de corrente de modo a
simular a espessura de deslocamento (Fig. 4.10). Uma vez que, neste momento, ja sao conhecidas as
evolucoes das espessuras de deslocamento para cada seccao da asa, a intensidade σvisk para o painel
k e determinada pela Eq. 4.5.
Como ja foi referido anteriormente, o XFOIL pode nao conseguir devolver uma solucao. Neste caso,
os paineis que pertencem a esta seccao vem com σvis = 0. Isto significa que, para a seccao em que o
XFOIL nao forneca uma solucao, aquela comporta-se como imersa num escoamento invıscido.
4.5 Calculo das Forcas Aerodinamicas
Depois de se adicionar o termo que simula o efeito da espessura de deslocamento σvisl , a solucao
final µ e alterada em relacao a solucao µinv e portanto, a forca resultante do corpo sustentador (Seccao
35
2.5 do Capıtulo 2), para alem de contabilizar a resistencia induzida, agora aparece um novo termo
denominado por resistencia de pressao. Esta resistencia de pressao resulta da alteracao do campo
de pressoes associada a alteracao da forma do corpo real ao introduzir-se o efeito da espessura de
deslocamento no corpo. A soma destes dois termos com a resistencia de atrito, devido a tensao de
corte τw, resulta na componente total de forca na direccao do escoamento de aproximacao,
DT = (Dind +Dδ∗) +Dτw (4.20)
onde Dδ∗ e a resistencia devido a introducao da espessura de deslocamento δ∗ e Dτw a resistencia
devido a tensao de corte na superficie do corpo.
36
Capıtulo 5
Validacao e Discussao dos
Resultados
O objectivo desta seccao passa pela validacao do codigo desenvolvido. Para tal, primeiro, e apresen-
tado um estudo de convergencia onde pode-se perceber como e que a solucao varia em funcao do
numero de paineis e, em segundo, os resultados obtidos pelo metodo dos paineis 3D sao comparados
com resultados experimentais para tres asas de geometrias diferentes. A primeira comparacao e feita
para numa asa rectangular de geometria simples. De seguida, e realizada a comparacao com uma asa
onde as cordas nas extremidades sao diferentes da corda na raiz e, com a introducao de um pequeno
angulo de diedro. Por ultimo, a comparacao de um conjunto asa e flap com afilamento e angulo de
flecha.
5.1 Estudo de Convergencia
Esta seccao e dedicada ao estudo de convergencia dos resultados obtidos pelo metodo dos paineis 3D.
Para o efeito, analisou-se o modo como os coeficientes de pressao ao longo da corda e os coeficientes
de sustentacao e resistencia variam em funcao do numero de paineis.
Num primeiro momento, procedeu-se a construcao de quatro malhas com diferentes numeros de
paineis para a modelacao de uma asa rectangular com 1m de corda, 6m de envergadura e um perfil
NACA0012. Enquanto que na direccao da envergadura o espacamento e considerado constante, na
direcao da corda o mesmo varia segundo uma funcao conseno, de modo a permitir uma maior densi-
dade de paineis na zona do bordo de ataque e fuga. A funcao que define a posicao x dos pontos das
extremidades dos paineis que modelam a asa e a seguinte :
xi =c
2(1 + cosβi) , (5.1)
βi = iπ
ndiv, i = 0, 1, 2...ndiv (5.2)
37
Figura 5.1: Quatro malhas para o estudo da convergencia: a) M1=10x5, b) M2=20x11, c) M3=40x21, d)M4=80x41
onde c corresponde a corda do perfil e ndiv o numero de divisoes na direccao da corda, ou por outras
palavras, o numero de paineis.
Na Fig. 5.1 a) esta representada uma malha com 50 paineis com 10 na direccao da corda (5
extradorso + 5 intradorso) e 5 na envergadura. Fig. 5.1 b) 220 paineis, 20 na direccao da corda e 11 na
envergadura. Fig. 5.1 c) 840 paineis, 40 na direccao da corda e 21 na envergadura. Fig. 5.1 d) 3280
paineis, 80 na direccao da corda e 41 na envergadura. A partir deste momento estas quatro malhas
serao designadas por M1, M2, M3 e M4, respectivamente.
As simulacoes das asas da Fig.5.1, imersas no seio de um escoamento invıscio, foram realizadas
com recurso ao metodo dos paineis 3D desenvolvido neste trabalho para diferentes angulos de ataque.
Na Fig. 5.2 estao representados os coeficientes de pressao, a meia envergadura, em funcao da
percentagem de corda do perfil a diferentes angulos de ataque para as quatro malhas M1, M2, M3
e M4. As solucoes obtidas com a malha menos refinada, M1, revelam-se adequadas para todos os
angulos de ataque para posicoes proximas do bordo de fuga. No que diz respeito aos coeficientes de
pressao numa posicao proxima do bordo de ataque estes demonstram-se imprecisos no que concerne
a previsao do valor do pico de succao e a sua posicao. Para as malhas M2 e M3, os valores dos
coeficientes de pressao sao bastante similares aos obtidos pela malha M4, e apresentam uma boa
previsao do valor do pico de succao.
No que diz respeito aos coeficientes de sustentacao e resistencia, os resultados para as quatro
malhas sao apresentados na Tab. 5.1. Os resultados obtidos com a malha M4 servem de referencia
para a comparacao dos valores calculados pelas outras malhas. Na Fig. 5.3 a) e b) sao apresentados os
38
Figura 5.2: Coeficientes de pressao Cp em funcao da percentagem de corda %c: a) α = 0° , b) α = 5° ,c) α = 10° , d) α = 15°
αCL CD
M1 M2 M3 M4 M1 M2 M3 M45 0.4359 0.4058 0.3957 0.3933 0.0181 0.0106 0.0090 0.008410 0.8634 0.8072 0.7878 0.7833 0.0653 0.0380 0.0333 0.032215 1.2742 1.1998 1.1729 1.1666 0.1412 0.0827 0.0731 0.071320 1.6612 1.5796 1.5476 1.5403 0.2415 0.1429 0.1273 0.1247
Tabela 5.1: Coeficiente de sustentacao e resistencia para as quatro malhas
comportamentos dos erros relativos, para o coeficiente de sustentacao e resistencia respectivamente,
em funcao do numero de paineis. O erro relativo para o coeficiente de sustentacao e resistencia foi
definido como :
εCL =
∣∣∣∣CLM4 − CLMi
CLM4
∣∣∣∣ , (5.3)
εCD =
∣∣∣∣CDM4 − CDMi
CDM4
∣∣∣∣ . (5.4)
com i=1,2 e 3.
Em relacao aos erros relativos do coeficiente de sustentacao, observa-se que para todos os angulos
de ataque ha um decrescimo acentuado na passagem da malha M1 para M2. Para a malha M1, os erros
relativos para os quatro angulos de ataque apresentam-se numa gama entre aproximadamente 8% a
11%; Para a malha M2 situam-se entre 2.5% a 3% e, por fim, na malha M3, sao proximos de 0.5%.
Para o erro relativo do coeficiente de resistencia, o comportamento e o mesmo que no coeficiente
39
Figura 5.3: a) erro relativo CL Vs Nº Paineis b) erro relativo CD Vs Nº Paineis c) Tempo Vs Nº Paineis
de sustentacao. Isto significa que ha um descrescimo acentuado do erro relativo da malha M1 para a
M2. Os resultados para os quatro angulos de ataque passam de uma gama entre 95% a 115% para
15% a 25%. Para a malha M3 os resultados sao inferiores a 10%.
Na Fig. 5.3 c) apresenta-se o tempo em segundos para a obtencao de uma solucao em funcao do
numero de paineis. O tempo computacional aumenta com a melhoria da precisao, que e conseguida
com o aumento do numero de paineis na malha.
As malhas utilizadas nas proximas seccoes foram escolhidas de modo a que o numero de elementos
fique situado entre o numero de elementos da malha M3 e M4, onde o coeficiente de sustentacao e
resistencia apresentam valores bastante proximos dos resultados para uma malha bastante refinada.
5.2 Efeito Solo
De modo a confirmar se os fenomenos apresentados no capıtulo do Efeito solo foram simulados de
forma correcta, os resultados obtidos pelo metodo dos paineis foram comparados com resultados do
programa comercial STAR-CCM+. Os dados e os resultados do problema podem ser consultados em
[26]. A simulacao foi feita para uma asa rectangular com um perfil simetrico simples NACA0015, com
aspect ratio de 6.6 e um numero de Reynolds 1.5× 106.
Para a simulacao onde se utilizou o metodo dos paineis 3D, varios angulos de ataque foram anali-
sados para diferentes distancias h do bordo de fuga ao solo , como e mostrado na Fig. 5.4 e com os
respectivos resultados nas Figuras 5.5 e 5.6. Na primeira figura, esta representada o esquema da asa
com a sua imagem. Na segunda, Fig. 5.5, apresenta-se uma sobreposicao dos resultados obtidos pelo
metodo dos paineis e os do STAR-CCM+ onde o coeficiente de sustentacao varia em funcao de h/c,
40
Figura 5.4: Representacao da asa e imagem para simular o efeito solo pelo metodo dos paineis 3D
com c igual a corda da asa. Na ultima figura (Fig. 5.6) o coeficiente de sustentacao surge em funcao do
angulo de ataque α. Nas figuras 5.5 e 5.6 e possıvel observar que os resultados obtidos pelo metodo
dos paineis seguem a mesma tendencia que os obtidos pelo software STAR-CCM+.
A maior diferenca, no que diz respeito aos resultados do coeficiente de sustentacao obtidos por am-
bos os metodos, ocorre quando a asa se encontra proxima do solo com um baixo angulo de ataque. Isto
deve-se ao facto da solucao do metodo dos paineis, contrariamente a do STAR-CCM+, nao contabilizar
a viscosidade. Assim, as velocidades desenvolvidas entre o solo e o intradorso sao bastante elevadas,
originando um coeficiente de sustentacao muito negativo (forca descendente) como mostra a Fig. 5.5
para um angulo de ataque igual a 0° . Neste caso em especıfico o coeficiente de sustentacao previsto
pelo metodo dos paineis e de -3.3 enquanto que para o STAR-CCM+ e aproximadamente 0.94.
Na Fig. 5.7 apresentam-se as distribuicoes dos coeficientes de pressao a meia envergadura da asa
para dois angulos de ataque, α = 0° e α = 8° , a duas distancias do solo, h/c = 1 e h/c = 0.1. Na
Fig. 5.7 a) os resultados sao bastante proximos dos obtidos pelo STAR-CCM+, o que era de esperar
uma vez que o efeito solo ainda nao se faz sentir a um h/c = 1, e tambem porque para um angulo
de ataque α = 0° as espessuras de deslocamento δ∗ ao longo da asa sao pequenas fazendo com
que os efeitos viscosos sejam quase inexistentes. Na Fig. 5.7 b) , nota-se que a distribuicao dos
coeficientes de pressao no extradorso seguem os do STAR-CCM+, enquanto que no intradorso, onde
ocorre o pico de succao, os resultados nao sao proximos. Os coeficientes de pressao mınimos sao
aproximadamente -4 e -11 para o STAR-CCM+ e metodos dos paineis, respectivamente. E ocorrem na
mesma zona, aproximadamente 30% da corda, o que ja era espectavel uma vez que e a posicao da
espessura maxima de um perfil NACA de 4 digitos e portanto, onde ocorre a distancia minima entre o
solo e a superfıcie da asa. Nas outras duas Figuras, 5.7 c) e 5.7 d), os resultados conseguidos pelo
metodo dos paineis aproximam-se bem dos do STAR-CCM+. Com a aproximacao ao solo observa-se
um aumento da area entre as distribuicoes de pressoes para o extradorso e intradorso. Traduzindo-se
num aumento do coeficiente de sustentacao.
E de realcar que para esta asa com um perfil simetrico NACA0015, o efeito solo e o efeito Venturi
podem ser observados dependendo do angulo de ataque a que a asa se encontra. Para angulos de
41
Figura 5.5: Comparacao dos resultados obtidos pelo metodo dos paineis com resultados STAR-CCM+
ataque inferiores a 4 o efeito Venturi e dominante em relacao ao efeito solo. Por outro lado, para angulos
superiores a 4° , o coeficiente de sustentacao aumenta a medida que a distancia ao solo diminui, sendo
o efeito solo o fenomeno dominante. Para um angulo de ataque igual a 4° verifica-se que a medida que
a asa se aproxima do solo os dois efeitos estao presentes. Ate a um valor de h/c = 0.3 o coeficiente de
sustentacao aumenta ligeiramente, enquanto que a partir deste valor o efeito Venturi torna-se dominante
e, consequentemente, o coeficiente de sustentacao diminui.
Na Fig. 5.6, para um h/c = 1 a variacao do coeficiente de sustentacao com o angulo de ataque e
de forma linear, semelhante a um caso onde o solo nao esta presente. Isto significa que a influencia do
solo na asa, para este h/c, nao e significativa.
No que diz respeito ao coeficiente de resistencia, este diminui quando se esta perante o fenomeno
de efeito solo. Por outro lado, quando o efeito Venturi e o fenomeno dominante o coeficiente de re-
sistencia aumenta. Na Fig. 5.8 estao os graficos CL Vs CD para os resultados obtidos pelo STAR-
CCM+ e pelo metodo dos paineis. A diferenca entre estes graficos e significativa ja que, pelo metodo
dos paineis, a unica resistencia que e tida em conta e a induzida; no STAR-CCM+ para alem da in-
duzida e, ainda, contabilizada a resistencia devida a tensao de corte que se desenvolve ao longo da
superfıcie da asa e a resistencia de pressao que esta relacionada com a formacao da espessura de
deslocamento. Nao obstante, os resultados do metodo dos paineis seguem a mesma tendencia que os
obtidos pelo STAR-CCM+. Por exemplo, para um angulo de ataque igual a zero α = 0° e a medida que
a asa se aproxima do solo, o coeficiente de resistencia aumenta significativamente; para um angulo de
ataque mais elevado, onde se faca sentir o efeito solo, e.g. α = 8° , o coeficiente de resistencia diminui
com a aproximacao ao solo.
42
Figura 5.7: Distribuicao do coeficiente de pressao a meia envergadura. a) α = 0 e h/c = 1 b) α = 0 eh/c = 0.1 c) α = 8 e h/c = 1 d)α = 8 e h/c = 0.1
44
Figura 5.9: Asa na vertical
Assim, o metodo das imagens, incorporado no codigo do metodo dos paineis 3D desenvolvido neste
trabalho, simula com uma precisao adequada a alteracao das caracterısticas aerodinamicas de uma asa
devido a presenca de uma superfıcie solida.
5.3 Asa na Vertical
De forma a perceber o modo como a agua influencıa as caracterısticas aerodinamicas de uma vela foi
efectuada uma nova simulacao mas, desta vez, com a asa colocada na vertical e com a superfıcie junto
a uma das extremidades (Fig. 5.9).
Analogamente ao que foi feito na seccao anterior, executaram-se varias simulacoes para varios
angulos de ataque da asa e para diferentes alturas ao solo. Os resultados estao representados nos
graficos na Fig. 5.10, nos quais constam, tambem, os valores dos coeficientes de sustentacao caso
o solo nao estivesse presente, dando uma boa nocao do quao esta caracterıstica e alterada com a
aproximacao ao solo.
Repare-se que apesar da diferenca entre o coeficiente de sustentacao com e sem presenca do solo
aumentar com o angulo de ataque, em termos de ganho percentual e sempre o mesmo. Por exemplo,
para um h/c = 0.1 e para um angulo de ataque igual a 2° , o ganho em termos de coeficiente de
sustentacao com a aproximacao ao solo e de(CL/CLref
)α=2
= 0.1777/0.1622 = 1.0956. E, para um
angulo de ataque igual a 6 e o mesmo h/c,(CL/CLref
)α=6
= 0.5321/0.4857 = 1.0956. Onde CLrefe o
coeficiente de sustentacao para o caso em que o solo nao esta presente.
Os resultados dos coeficientes de sustentacao e resistencia em relacao as suas referencias, para
os diferentes valores de h/c, estao tabelados em 5.2. Para esta asa, com a diminuicao da distancia
ao solo, o coeficiente de sustentacao pode aumentar ate aproximadamente 9.5% enquanto que para o
46
h/c CL/CLrefCDi
/CDiref
10 1.0011 0.99765 1.0040 0.99181 1.0264 0.9506
0.5 1.0435 0.92300.3 1.0585 0.90130.2 1.0715 0.88400.1 1.0956 0.8542
Tabela 5.2: Variacao do coeficiente de sustentacao e resistencia para uma asa rectangular com umasuperfıcie solida junto a uma das extremidades
coeficiente de resistencia sofre um decrescimo de quase 15%.
Conclui-se que mesmo com a asa na vertical a influencia do solo continua a fazer-se sentir de forma
significativa nos coeficientes de sustentacao e resistencia.
5.4 Asa com Flap
Para demonstrar que o codigo desenvolvido e capaz simular um escoamento invıscido em torno de
duas superficies sustentadoras, serao apresentados os resultados, para um conjunto asa e flap, de
duas situacoes distintas. Na primeira, as superfıcies foram colocadas a uma distancia elevada entre si,
enquanto que na segunda, o conjunto asa e flap operam proximos um do outro.
Considerou-se duas asas rectangulares exactamente iguais com uma corda de 1m, envergadura 6m
e um perfil NACA0012. Estas asas foram simuladas com duas configuracoes diferentes. Na primeira
o flap foi colocado a uma distancia de 50 cordas da asa, enquanto que na segunda, a distancia foi de
apenas 1 corda.
Na Fig. 5.11 estao representados os coeficientes de pressao na asa e no flap para as duas
configuracoes consideradas. Isto para um angulo de ataque do escoamento de aproximacao igual
a 10° .
Como esperado, a uma distancia de 50 cordas (Configuracao 1 na Fig 5.11), a asa e o flap comportaram-
se como se estivessem a operar isoladamente no escoamento, uma vez que apresentam coeficientes
de pressao identicos. Estas superfıcies foram colocadas a uma distancia tao grande uma da outra que
a influencia dos paineis de um corpo nao se faz sentir no outro. A distribuicao do coeficiente de pressao
para estas duas superficies e identica a distribuicao apresentada no capitulo Estudo de Convergencia
(Fig. 5.2 c))
Na segunda configuracao (Configuracao 2 na Fig 5.11) o flap foi colocado a 1 corda de distancia do
bordo de ataque da asa, com um desvio de 0.2c na direccao vertical para que nao houvesse a hipotese
de interseccao com a esteira, como ilustrado na Fig. 5.12 . No que diz respeito aos resultados, observa-
se que a distribuicao dos coeficientes de pressao na asa e flap sao alterados significativamente para
os dois casos considerados. Na asa, na seccao a meia envergadura o pico de succao tornou-se mais
intenso passando de -4 para aproximadamente -6, enquanto que para o flap, de -4 para -0.5.
A distribuicao do coeficiente de pressao para a asa e flap nas configuracoes 1 e 2 estao ilustradas
48
Figura 5.11: Coeficiente de pressao nas duas configuracoes para a) asa e b) flap
Figura 5.12: Ilustracao do conjunto asa e flap
49
Figura 5.13: Distribuicao de Cp na asa e flap para a configuracao 1
nas Figuras 5.13 e 5.14, respectivamente. Note-se que nestas figuras a distancia do flap relativamente
a asa nao corresponde a real distancia utilizada nas respectivas simulacoes, com isto se pretendendo
uma melhor visualizacao.
Fica assim demonstrado que o codigo desenvolvido contempla a interaccao entre as duas superfi-
cies imersas num escoamento invıscido.
5.5 Comparacao com Resultados Experimentais
Nesta seccao serao apresentadas algumas comparacoes de resultados experimentais com resultados
obtidos pelo metodo dos paineis 3D com e sem acoplamento dos efeitos viscosos.
5.5.1 Asa Rectangular
Procedeu-se a comparacao de resultados experimentais com os numericos para uma asa rectangular
com 6 metros de envergadura, 1 metro de corda e com um perfil NACA 0012 (Fig. 5.15). A activi-
dade experimental esta descrita em [40] no qual consta a distribuicao de pressoes ao longo de varias
estacoes na direccao da envergadura para uma gama de angulos de ataque de −6° a 20° a um Rey-
nolds igual a 3 × 106 com base na corda da asa. O referido documento inclui, ainda, um grafico do
50
Figura 5.14: Distribuicao de Cp na asa e flap para a configuracao 2
Figura 5.15: Perfil NACA0012 utilizado na construcao da asa rectangular
51
Figura 5.16: Malha na asa rectangular
α CLexp CLinv erro relativo [%]2 0.1562 0.1590 1.796 0.4573 0.4762 4.1310 0.7580 0.7908 4.3215 1.1115 1.1776 5.9520 1.4296 1.5543 8.72
Tabela 5.3: Resultado do coeficiente de sustentacao experimental e invıscido com o respectivo errorelativo para um Reynolds 3× 106
coeficiente de sustentacao em funcao do angulo de ataque, para varios numeros de Reynolds.
Primeiramente compararam-se os resultados invıscidos dos coeficientes de pressoes para uma
seccao a meio da envergadura e outra junto a extremidade, assumindo-se um angulo de ataque igual
4.64° e outro a 10.97° . A malha utilizada, de onde se obteve a solucao numerica pelo metodo dos
paineis 3D, possui 60 paineis na direccao da corda e 20 na direccao da envergadura como mostra a
Fig. 5.16.
Note-se que os resultados obtidos com o metodo dos paineis no que diz respeito a distribuicao de
pressao, tanto a meia envergadura como junto a uma das extremidades, para os dois angulos de ataque
considerados, 4.64° e 10.97° (Fig. 5.17 e 5.18), sao bastante proximos dos resultados experimentais.
Para os dois angulos de ataque, os resultados obtidos para a seccao a meia envergadura apresentam
uma correspondencia quase perfeita com os resultados experimentais. Permitindo prever, com um ele-
vado grau de precisao, o valor e a posicao do pico de succao. Na seccao mais proxima da extremidade
da asa, para o angulo de ataque igual a 10.97° , o valor do pico de succao continua a ser bem previsto.
Por outro lado, para o angulo de ataque de 4.64° , os valores dos coeficientes de pressao na zona do
bordo de ataque ficam aquem dos experimentais. Para valores proximos do bordo de fuga os resulta-
dos numericos, tanto para 10.97° como 4.64° , afastam-se dos experimentais em termos de intensidade,
mas no que diz respeito a evolucao da distribuicao de pressao demonstram ser muito semelhantes ao
caso experimental.
A Fig. 5.19 representa a evolucao do coeficiente de sustentacao em funcao do angulo de ataque
cujos valores, com o respetivo erro relativo, constam da Tab 5.3, de onde e possıvel observar que os
valores comparados sao bastante proximos. Por outro lado, verifica-se que os erros relativos aumentam
a medida que o angulo de ataque tambem aumenta, nunca ultrapassando os 9%.
52
Figura 5.17: Coeficientes de pressao para duas seccoes a um angulo de ataque de 4.64º com Reynoldsigual a 3× 106
Figura 5.18: Coeficientes de pressao para duas seccoes a um angulo de ataque de 10.97º com Rey-nolds igual a 3× 106
53
Figura 5.19: CL vs alpha para uma asa rectangular a um Reynolds igual a 3× 106
α CLexp CLinv erro[%] CLinv+vis erro[%]
2 0.1464 0.1590 8.60 0.1435 1.986 0.4338 0.4762 9.77 0.4276 1.4310 0.7307 0.7908 8.22 0.7367 0.8215 1.0790 1.1776 9.13 1.0834 0.4018 1.2558 1.4157 12.73 1.2399 1.2620 1.1428 1.5543 36 1.3142 15
Tabela 5.4: Resultado experimental, invıscido e viscoso acoplado ao invıscido para um Reynolds 1.92×106
A justificacao para os baixos erros relativos assenta no facto de o resultados experimentais apre-
sentarem uma evolucao com o angulo de ataque muito proxima da linearidade (Fig. 5.19) o que, por
sua vez, deve-se ao elevado numero de Reynolds.
Na comparacao dos resultados numericos - onde sao tidos em conta os efeitos viscosos - com
os resultados experimentais considerou-se um numero de Reynolds igual a 1.92× 106 do procedimento
experimental. Tal escolha assentou no facto de ser o menor Reynolds da atividade experimental e, como
tal, e o que se traduz em maior influencia viscosa no escoamento. Os resultados, tanto numericos como
experimentais, encontram-se na Tab. 5.4 e ilustrados na Fig. 5.20.
Numa primeira analise da Tab. 5.4, e comparando com o caso anterior, observa-se que o erro rela-
tivo aumentou numa comparacao entre o escoamento inviscido e experimental. Este aumento deve-se
ao facto da curva CLvsα da actividade experimental para este novo Reynolds ter diminuido de declive.
O que assenta na circunstancia de os efeitos viscosos serem mais intensos quando comparados com
o escoamento a um Reynolds 3× 106
Nas duas ultimas colunas os resultados para um escoamento viscoso calculado numericamente,
onde os efeitos viscosos sao acoplados a solucao inviscida, mostram-se bastante proximos dos resulta-
dos da actividade experimental. Ate ao angulo da sustentacao maxima, 18° , os resultados apresentam
um erro relativo maximo de 2%; depois deste angulo, a diferenca de resultados torna-se evidente.
54
Figura 5.20: CL vs alpha para uma asa rectangular a um Reynolds igual a 1.92× 106
Figura 5.21: Perfil NACA65-210 utilizado na construcao da asa com afiamento
5.5.2 Asa com Afilamento
Nesta seccao serao comparados os resultados obtidos numa actividade experimental, descrita em [41],
com a solucao obtida numericamente quer para um caso invıscido quer viscoso.
Os dados da geometria da asa utilizada na actividade experimental e na simulacao numerica podem
ser consultados na Tab. 5.5 e ilustrados na Fig. 5.22. A asa foi modelada com uma malha com 60
paineis na direccao da corda e 20 na direccao da envergadura (Fig. 5.23).
Os resultados do coeficiente de sustentacao para a actividade experimental juntamente com a
solucao numerica para um escoamento invıscido e viscoso, com o respectivo erro em relacao aos
primeiros, sao apresentados na Tab 5.6 e, em formato de grafico na Fig. 5.24 a).
No que diz respeito aos coeficientes de sustentacao do escoamento invıscido verifica-se uma boa
aproximacao aos coeficientes de sustentacao experimentais, com erros relativos, para angulos de ata-
Perfil NACA 65-210 (Fig. 5.21)corda na raiz cr 0.726 m
corda na extremidade ct 0.290 mrazao de afilamento= ct
cr0.4
angulo de diedro 0.3ºReynolds 4.4× 106
Tabela 5.5: Dados da asa com afilamento
55
Figura 5.22: Esquema da asa com afilamento
Figura 5.23: Malha na asa com afilamento
α CLexp CLinv erro relativo[%] CLinv+vis erro relativo[%]2 0.2937 0.3375 14.91 0.3233 10.086 0.6148 0.7053 14.72 0.6623 7.73
10 0.9697 1.0712 10.46 0.9977 2.8912 1.1228 1.2530 11.59 1.1776 4.8814 1.0060 1.4338 42.52 1.2958 28.81
Tabela 5.6: Resultado experimental, invıscido e viscoso acoplado ao invıscido para um Reynolds 4.4×106
56
Figura 5.24: a) CL Vs alpha, b) CL Vs CD. Para asa com afilamento
que ate aproximadamente 13º, a nao ultrapassarem os 15%. Para o angulo de ataque 14° , e visısvel
que a asa ja se encontra em perda, e portanto, o valor da solucao invıscida encontra-se bastante afas-
tada da experimental.
Quando os efeitos viscosos sao acoplados a solucao inviscida os resultados melhoram conside-
ravelmente ate ao angulo de sustentacao maxima. Os erros relativos para os angulos de ataque 2°
e 6° passam de 14.91% e 14.72%, numa situacao de escoamento inviscido, para 10.08% e 7.73%,
respectivamente. Depois do angulo de ataque igual a 6° e ate aproximadamente 13° , onde ocorre a
separacao do escoamento, os erros mantem-se baixos nao sendo superiores a 5%. Para angulos de
ataque superiores a 13° o coeficiente de sustentacao diminui de forma acentuada para os resultados
experimentais, representando a entrada em perda da asa, enquanto que o resultado numerico nao con-
segue acompanhar esse descrescimo originando um erro relativamente elevado mas, bastante inferior
que no caso invıscido.
Na Fig. 5.24 b) estao ilustradas as evolucoes do coeficiente de resistencia para a actividade experi-
mental, para o escoamento invıscido e para a solucao viscosa acoplada ao escoamento invıscido.
Repare-se que no calculo do coeficiente de resistencia para o escoamento invıscido, apenas e
contabilizada a resistencia induzida. Daı a sua anulacao para uma situacao em que a sustentacao e
proxima de zero, ja que nao ha libertacao de vortices das extremidades da asa que sao os responsaveis
pela inducao de velocidades descendentes no escoamento de aproximacao.
No escoamento em que sao acoplados os efeitos viscosos as resistencias devido ao facto de se
considerar a espessura de deslocamento (Dδ∗ ) e a tensao de corte na superficie da asa (Dτw ) sao
contabilizadas, resultando numa maior concordancia da curva CL Vs CD com a esperimetal.
5.5.3 Asa RAE com Flap
A asa RAE testada em [42], e construida com um perfil com uma expessura maxima de 11.7% a uma
posicao de 37.5% e com uma curvatura maxima a ocorrer a 75%, percentagens em relacao a corda. Os
57
Figura 5.25: Representacao das variaveis que definem a distancia do flap a asa.[42]
10° 25°g/c 0.01300 0.01864h/c 0.06295 0.03530
Tabela 5.7: Parametros relativamente a distancia do flap ao bordo de fuga da asa para dois angulos dedeflexao, 10° e 25°
pontos utilizados para a construcao deste perfil, assim como o do flap, estao disponıveis no documento
da actividade experimental. A asa tem uma envergadura de 2.14 m com uma razao de afilamento de
0.35 e nao possui angulo de torcao nem de diedro mas tem um angulo de flecha de 28° com a linha
que define 1/4 de corda. O flap apresenta uma corda de 34% da corda local da asa. A actividade
experimental foi realizada a um Reynolds de 1.35× 106 definido em relacao a corda media, 0.26 m.
Os resultados experimentais com os quais vao ser comparados os resultados numericos foram
obtidos para duas situacoes: na primeira, o flap apresenta um angulo de 10° e, na segunda, angulo de
25° . Para estes dois angulos, o flap apresenta diferentes distancias em relacao ao bordo de fuga da
asa. Na Fig. 5.25 estao representadas as variaveis que, em conjunto, vao definir a distancia do flap ao
bordo de fuga da asa para as duas situacoes, onde g e a distancia do ponto mais alto do extradorso do
flap ao bordo de fuga da asa e h do bordo de ataque do flap ao bordo de fuga da asa. Na Tab. 5.7 sao
apresentados estes valores para os dois angulos, 10° e 25° .
Atendendo a que a asa e simetrica procedeu-se a simulacao de apenas metade da geometria, sendo
a outra metade simulada pelo metodo das imagens. Ou seja, e utilizado o mesmo metodo desenvolvido
na seccao do Efeito Solo mas, agora, com uma altura ao solo de h = 0, como ilustrado na Fig. 5.26.
A malha utilizada na asa e flap e contruıda com 100 paineis na direccao da corda (50 extradorso +
50 intradorso) e 20 na direccao da envergadura. Nesta simulacao foi necessario uma malha com um
maior numero de paineis na direccao da corda do que o habitual. Na verdade, para uma malha menos
refinada, o XFOIL nao devolvia uma solucao que tivesse convergido e portanto, o acoplamento dos
efeitos viscosos a solucao inviscida nao podia ser realizada. Com o aumento do numero de paineis
na direccao da corda, o XFOIL recebe um maior numero de pontos de controlo e consequentemente o
perfil a ser analisado tem uma geometria mais “suave”. Este aumento de pontos para a construcao do
perfil traduz-se na obtencao de uma solucao final convergente por parte do XFOIL.
Por outro lado, repare-se que nas Fig. 5.26 e 5.27 esta ilustrado o conjunto asa e flap com uma ma-
lha bastante inferior a utilizada na simulacao numerica, com o que se pretendeu melhorar a visualizacao
da geometria.
Na Tab. 5.8 estao representados os resultados experimentais, inviscidos e viscosos para a situacao
58
Figura 5.26: Conjunto asa flap para simulacao, com a respectiva imagem
Figura 5.27: Conjunto asa flap
59
αFlap 10°
CLexp CLinv e [%] e[%] CLinv+vis e [%] e[%]0 0.5453 0.6153 12.8
11.2
0.5750 5.4
6.72 0.7079 0.8044 13.6 0.7737 9.34 0.8983 0.9922 10.5 0.9621 7.16 1.0774 1.1780 9.3 1.1398 5.88 1.2377 1.3617 10.0 1.3119 6
10 1.2986 1.5427 18.8 - 1.4944 15 -11 1.1726 1.6321 39.2 - 1.5618 33.2 -
Tempo - ≈ 20min - - ≈ 1h20min - -
Tabela 5.8: Resultados experimentais e numericos para uma deflexao do flap de 10°
Figura 5.28: Coeficientes de sustentacao e resistencia para um angulo de deflexao 10° para o flap. a)CL vs α, b) CL vs CD
em que o flap apresenta uma deflexao de 10° . Estes resultados estao ilustrados na Fig. 5.28 a).
Para o flap com um angulo de 10° , os resultados experimentais apresentam uma evolucao apro-
ximadamente linear, ate a um angulo de ataque igual a 8° . A partir daı, o coeficiente de sustentacao
aumenta de uma forma menos acentuada ate a um angulo de ataque de 10° altura em que a asa
comeca a entrar em perda provocando uma diminuicao do coeficiente de sustentacao.
Os resultados da simulacao de um escoamento invıscido em torno do conjunto asa e flap apresentam-
se satisfatorios tendo em conta a complexidade da geometria analisada. Para angulos de ataque ate
aproximadamente 8° os erros relativos apresentam-se na sua maioria superiores a 10% e nao ultrapas-
sam os 13.6%. Com a entrada da asa em perda os resultados inviscidos, tendo uma evolucao linear,
afastam-se rapidamente e apresentam erros relativos de 18.8% e 39.2% para os angulos de 10° e 11° ,
respectivamente.
Na simulacao em que os efeitos viscosos sao tidos em consideracao, a evolucao dos resultados
e similar a solucao invıscida mas, por apresentar valores inferiores, os erros relativos sao menores.
Ate a um angulo de ataque igual a 8° os erros tomam valores que nao ultrapassam os 9.3%. Uma
vez que a tendencia dos resultados em que sao contemplados os efeitos viscosos e a mesma que
a dos resultados inviscidos, depois de um angulo de ataque igual a 8° os erros relativos aumentam
consideravelmente. Em suma, a introducao dos efeitos viscosos na simulacao permite uma previsao da
60
αFlap 25°
CLexp CLinv e [%] e[%] CLinv+vis e [%] e[%]0 1.2566 1.3881 10.5
10.5
1.3287 5.7
6.42 1.4458 1.5730 8.8 1.5120 4.64 1.5793 1.7555 11.2 1.7020 7.86 1.7369 1.9352 11.4 1.8647 7.48 1.6128 2.1116 30.9 - 2.0426 26.6 -
Tempo - ≈ 17min - - ≈ 1h - -
Tabela 5.9: Resultados experimentais e numericos para uma deflexao do flap de 25°
Figura 5.29: Coeficientes de sustentacao e resistencia para um angulo de deflexao 25° para o flap. a)CL vs α, b) CL vs CD
curva CL vs α com maior precisao.
Na Fig. 5.28 b) esta ilustrada a curva CL vs CD para os resultados experimentais, inviscidos e
com o acoplamento dos efeitos viscosos. E visivel o melhoramento introduzido na previsao desta curva
apos a introducao dos efeitos efeitos viscosos a solucao invıscida. Os resultados para o coeficiente de
resistencia depois da asa entrar em perda nao sao bem previstos.
Na Tab. 5.9 estao representados os resultados experimentais, inviscidos e viscosos para a situacao
em que o flap tem uma deflexao de 25° . Estes resultados estao ilustrados na Fig. 5.29 a).
Tal como no caso anterior, os erros relativos para os resultados de um escoamento inviscido sao, na
sua maioria, superiores a 10% com um tecto maximo de 11.4%. Depois de um angulo de ataque igual
a 6° o coeficiente de sustentacao diminui a sua razao de crescimento e atinge o seu valor maximo a 7°
apos o qual entra em perda. Assim, depois de um angulo de ataque igual a 6° o erro relativo para o
caso inviscido aumenta significativamente.
Para a simulacao em que sao acoplados os efeitos viscosos a solucao inviscida, e ate a um angulo
de ataque de 6° , os erros relativos correspondem, em termos mınimos, a 4.6% e nao ultrapassam os
7.8%. Ja para um angulo de ataque igual a 8° , onde a asa ja se encontra em perda, o erro passa de
um valor de 30.9%, no caso inviscido, para 26.6%.
Na Fig. 5.29 b) esta ilustrada a curva CL vs CD para os resultados experimentais, inviscidos e com
o acoplamento dos efeitos viscosos. Mais uma vez, a solucao do escoamento viscoso aproxima-se com
uma boa precisao dos dos resultados experimentais
61
Para a situacao em que o flap apresenta uma deflexao de 25° os resultados, tanto invıscido como
os que consideram os efeitos viscosos, mantem, em media, aproximadamente os mesmos valores
de erros ate ao angulo de ataque em que o coeficiente de sustentacao, experimental, respeite uma
evolucao linear. Ou seja, em media o erro relativo para um escoamento, a um angulo de ataque no
interior da gama de angulos de ataque que respeitem uma evolucao linear por parte do coeficiente de
sustentacao experimental, inviscido e viscoso sao de 10.5% e 6.4%, respectivamente. E para a situacao
anterior com flap a 10° , 11.2% e 6.7%.
Concluindo, para as duas situacoes analisadas, 10° e 25° , o aumento da precisao da previsao dos
coeficientes de sustentacao com a utilizacao do acoplamento dos efeitos viscosos e evidente e, ainda
maior no que diz respeito a previsao do coeficiente de resistencia. Isto para uma situacao em que a asa
nao apresenta sinais de separacao. Quando a asa comeca a apresentar sinais de entrada em perda
os resultados nao sao bem previsto. Este melhoramento da solucao e acompanhado de um aumento
de tempo na resolucao do problema em aproximadamente 4 e 3.5 vezes nas situacoes de 10° e 25° ,
respectivamente. A simulacao foi efectuada num computador portatil com processador i7 2.30GHz e
8GB de RAM.
62
Capıtulo 6
Conclusoes e Trabalhos Futuros
O principal objectivo deste trabalho passava pela construcao de uma ferramenta computacional que
conseguisse calcular as caracteristicas aerodinamicas de um conjunto asa e flap.
A simulacao de um escoamento tridimensional invıscido, incompressivel e irrotacional em torno
de um corpo sustentador foi conseguido atraves da aplicacao do metodo dos paineis 3D, implemen-
tado em MATLAB. Este metodo revelou-se rapido no calculo da solucao final e preciso no que toca a
comparacoes com resultados experimentais.
Ao metodo dos paineis 3D foi introduzido o metodo das imagens para a simulacao de escoamen-
tos em torno de corpos sustentadores a operarem na proximidade de uma superficie plana. Com a
aplicacao deste metodo, demonstrou-se a capacidade de previsao de dois fenomenos caracteristicos
destes escoamentos: efeito solo e Venturi. Relativamente ao primeiro, verificou-se que a variacao das
caracteristicas aerodinanmicas sao predictas com uma boa precisao. Ja para o efeito Venturi, a previsao
da intensidade dos coeficientes de sustentacao e resistencia nao foi a mais correcta. Viu-se que essa
discrepancia nos resultados estava relacionada com o facto do metodo dos paineis nao contabilizar os
efitos viscosos. No entanto, as tendencias verificadas para o coeficiente de sustentacao e resistencia
demonstraram-se bastante satisfatorias quando comparados com resultados CFD.
Com a utilizacao deste metodo mostrou-se que uma asa na vertical com uma superficie na proximi-
dade de uma das extremidades, o coeficiente de sustentacao pode aumentar em 9.5% e o coeficiente
de resistencia induzida sofrer um decrescimo de 15%.
A formulacao do metodo dos paineis foi ligeiramente alterada com o objectivo de contabilizar a
variacao da velocidade no interior de uma camada limite atmosferica. Isto foi conseguido atraves da
adopcao de um modelo de variacao da velocidade segundo uma funcao ln.
Ao escoamento invıscido foram acomplados os efeitos viscosos atraves de uma velocidade de
transpiracao aplicada na superficie do corpo com o objectivo de simular o deslocamento das linhas
de corrente. Esta velocidade foi calculada a partir da variacao da espessura de deslocamento ao longo
da asa. A camada limite que se desenvolve na superficie do corpo sustentador admitiu-se bidimensio-
nal e portanto, a asa foi dividida em varias seccoes ao longo da envergadura onde foram tratadas cada
uma de forma independente. A analise a cada seccao foi realizada com recurso ao programa XFOIL.
63
Com a aplicacao deste acoplamento os resultados melhoraram significativamente, principalmente na
previsao do coeficiente de resitencia, em relacao aos invıscidos mas, ainda assim, sem a capacidade
de previsao da entrada em perda do corpo sustentador. A esta melhoria de resultados esta associado
um acrescimo no tempo de resolucao do problema em aproximadamente 4 vezes.
Concui-se assim que o codigo desenvolvido pode ser uma ferramenta bastante util numa fase preli-
minar de projecto de superficies sustentadoras, nomeadamente para velas rıgidas de catamaras, uma
vez que foi demonstrado neste trabalho boas aproximacoes a situacoes possıveis de operacao destas.
No que diz respeito a trabalhos futuros, este e um codigo que podera estar ligado a uma ferramenta
de optimizacao de geometria ou, de posicoes relativas enre as duas superficies sustentadoras imersas
no mesmo escoamento.
64
Bibliografia
[1] http://www.optimal.pt/Site/Home.html (05/12/2013).
[2] http://www.tonycastro.co.uk/ (05/12/2013).
[3] John L Hess and A.M.O. Smith. Calculation of non-lifting potential flow about arbitrary three-
dimensional bodies. Technical report, Douglas Aircraft Division, 1962.
[4] John L Hess. Calculation of potential flow about arbitrary three-dimensional lifting bodies. Technical
report, DTIC Document, 1972.
[5] Joseph Katz and Allen Plotkin. Low-Speed Aerodynamic: from wing theory to panel methods.
McGraw-Hill, 1991.
[6] R.L. Carmichael and L.L. Erickson. Pan air - a high order panel method for prediction subsonic
or supersonic linear potential flows about arbitrary configurations. In AIAA 14th Fluid and Plasma
Dynamics Conference, 1981.
[7] Brian Maskew. Prediction of subsonic aerodynamic characteristics: A case for low-order panel
methods. Journal of Aircraft, 19:157–163, 1982.
[8] Brian Maskew. Program vsaero theory document - a computer program for calculating nonlinear
aerodynamic characteristics of arbitrary configuration. Technical report, NASA, 1987.
[9] Dale L. Ashby. Potential flow theory and operation guide for the panel code pmarc. Technical report,
NASA, 1999.
[10] Mark Drela. Xfoil: An analysis and design system for low reynolds number airfoils. In Low Reynolds
number aerodynamics, pages 1–12. Springer, 1989.
[11] RJ Margason, SO Kjelgaard, WL Sellers III, CEK Morris Jr, KB Walkey, and EW Shields. Subsonic
panel methods- a comparison of several production codes. In AIAA 23 rd Aerospace Sciences
Meeting, 1985.
[12] William Michael Milewski. Three-dimensional viscous flow computations using the integral boun-
dary layer equations simultaneously coupled with a low order panel method. Master’s thesis, Mas-
sachusetts institute of technology, 1997.
65
[13] T. Cabeci and E. Besnard. An efficient and accurate approach for analysis and design of high lift
configurations. Canadian Aeronautics and Space Journal, 44:1–17, 1998.
[14] Arthur E.P. Veldman Edith G.M. Coenen and George Patrianakos. Viscous-inviscid interaction
method for wing calculations. In European Congress on Computational Methods in Applied Scien-
ces and Engineering, 2000.
[15] A. Veldman. New, quasi-simultaneous method to calculate interacting boundary layers. AIAA Jour-
nal, 19:79–85, 1981.
[16] Andre Deperrois. Analysis of foils and wings operating at low reynolds numbers.
[17] Anders Karlsson Kamran Mohseni Salvatore Lupo, Henrik Nyberg. Xwing - a 3d viscous design
tool for wings. In 46th AIAA Aerospace Sciences Meeting and Exhibit, 2008.
[18] Vasco de Brederode. Fundamentos de Aerodinamica Incompressıvel. Vasco de Brederode, 1997.
[19] Hanseong Lee. Modeling of Unsteady Wake Alignment and Developed Tip Vortex Cavitation. PhD
thesis, The University of Texas at Austin, 2002.
[20] Nilay Sezer-Uzol. High-accuracy wake and vortex simulations using a hybrid euler/discrete vortex
method. Master’s thesis, The Pennsylvania State University, 2001.
[21] Luis Vargas. Desenvolvimento e Implementacao de um Procedimento Numerico para Calculo de
Conjuntos Asa-Empenagens de Geometria Complexa em Regime de Voo Subsonico, Assimetrico
e Nao Linear. PhD thesis, Universidade Federal De Minas Gerais, 2006.
[22] John L. Hess and A.M.O. Smith. Calculation of potential flow about arbitrary bodies. Progress in
Aerospace Sciences, 8:1–138, 1966.
[23] Daniel Filkovic. Graduate work. Master’s thesis, Faculty of Mechanical Engineering and Naval
Architecture, 2008.
[24] Diogo Chaves. Metodo de painel para o calculo do escoamento potencial em torno de varios perfis
sustentadores. Master’s thesis, Instituto Superior Tecnico, 2012.
[25] T. Cabeci and J. Cousteix. Modeling and Computation of Boundary-Layer Flows. HORIZONS
PUBLISHING, 2005.
[26] Juhee Lee. Aerodynamic characteristics of rectangular wing with square laternal tip. In 51st
AIAA Aerospace Sciences Meeting Iincluding The New Horizons Forum And Aerospace Exposi-
tion, 2013.
[27] Th. Lutz S. Schmid and E. Kramer. Impact of modelling approaches on the prediction of ground
effect aerodynamics. 3:419–429, 2009.
[28] Nils Haack. C-class catamaran wing performance optimisation. Master’s thesis, University of Man-
chester, 2011.
66
[29] Nestor Ramos Garcıa. Unsteady Viscous-Inviscid Interaction Technique for Wind Turbine Airfoils.
PhD thesis, Technical University of Denmark, 2011.
[30] Edith Gerda Maria Coenen. Viscous-Inviscid Interaction with the Quasi-Simultaneous Method for
2D and 3D Aerodynamic Flow. PhD thesis, University of Groningen, 2001.
[31] Herve Martins-rivas David L. Rodriguez Peter Sturdza, Yoshifumi Suzuki. A quasi-simultaneous
interactive boundary-layer model for a cartesian euler solver. In 50th AIAA Aerospace Sciences
Meeting including the New Horizons Forum and Aerospace Exposition, 2012.
[32] H.A. Bijleveld and A.E.P. Veldman. Rotorflow: A quasi-simultaneous interaction method for the
prediction of three-dimensional aerodynamic flow over wind turbine blades.
[33] Mehmet Mersinligil. Airfoil boundary layer calculations using interactive method and e transition
prediction technique. Master’s thesis, The Graduate School of Natural and Applied Sciences Middle
East Technical University, 2006.
[34] Mark Drela. http://web.mit.edu/drela/public/web/xfoil/ (05/12/2013).
[35] Mark Drela and Michael B. Giles. Viscous-inviscid analysis of transonic and low reynolds number
airfoils. AIAA Journal, 25:1347–1355, 1987.
[36] Mark Drela. Xfoil 6.94 user guide.
[37] Gary S. Hufford. Viscous flow around marine propeller using boundary layer strip theory. Master’s
thesis, Massachusetts Institute Of Technology, 1992.
[38] Xiangming Yu. Three dimensional viscous/inviscid interactive method and its application to propeller
blades. Master’s thesis, The University of Texas at Austin, 2012.
[39] Gus Brown. http://www.mathworks.com/matlabcentral/fileexchange/30446-xfoil-interface
(05/12/2013).
[40] Long P. Yip and Gary L. Shubert. Pressure distributions on a 1 by 3 meter semispan wing at sweep
angles from 0 to 40 in subsonic flow. Technical report, NASA, 1976.
[41] James C. Sivells. Experimental and calculated characteristics of three wings of naca 64-210 and
65-210 airfoil sections. Technical report, NACA, 1947.
[42] D. A. Lovell. A wind-tunnel investigation of the effects of flap span and deflection angle, wing
planform and body on the high-lift performance of a 28º swept wing. Technical report, R.A.E. CP
1372, 1977.
67