Anais do 12O Encontro de Iniciação Científica e Pós-Graduação do ITA – XII ENCITA / 2006
Instituto Tecnológico de Aeronáutica, São José dos Campos, SP, Brasil, Outubro, 16 a 19, 2006
EXTENSÃO DO MÉTODO DOS PAINÉIS PARA O REGIME
COMPRESSÍVEL ATRAVÉS DA DISTRIBUIÇÃO DE SINGULARIDADES
UNICAMENTE NA FRONTEIRA DE PERFIS AERODINÂMICOS SEM
SUSTENTAÇÃO
Aline Sousa da Silveira
Instituto Tecnológico de Aeronáutica – ITA
Bolsista PIBIC-CNPq
[email protected]
André Valdetaro Gomes Cavalieri
Instituto Tecnológico de Aeronáutica – ITA
[email protected]
Paulo Afonso de Oliveira Soviero
Instituto Tecnológico de Aeronáutica – ITA
[email protected]
Resumo. Este trabalho procura estender para o regime transônico os métodos de análise aerodinâmica que distribuem soluções
elementares sobre a fronteira de perfis. Nesse regime, a equação mais simples possui um termo não-linear e, para resolvê-la para
um perfil, os métodos convencionais distribuem painéis bidimensionais sobre o campo de velocidades e fazem integrais duplas.
Usando o Teorema de Green, os dois métodos apresentados neste artigo substituem essas integrais duplas por integrais de linha
sobre a fronteira. O primeiro deles é totalmente uma formulação apenas sobre a fronteira, entretanto, não resolveu o problema
aerodinâmico, apesar de poder ser aplicado em outros problemas de equações similares, como em transferência de calor. O
segundo, que apresentou bons resultados no escoamento sobre um perfil parabólico, ainda necessita de valores no campo de
velocidades, mas não é necessário integrar esses valores, eles são usados apenas em cálculos algébricos. Esse último método
acabou se mostrando mais rápido e eficiente que os métodos convencionais, com resultados muito próximos.
Palavras chave: aerodinâmica transônica, método dos painéis, método da reciprocidade dual
1. Introdução
O método dos painéis, que procura resolver um escoamento unicamente por meio de integrais sobre a fronteira, é
hoje muito utilizado para o regime incompressível, cuja solução é extensível para o regime baixo subsônico por meio da
transformação de Prandtl-Glauert ou outras. Entretanto, essa extensão falha para números de Mach próximos de 1,
devido ao aparecimento de efeitos ligados à transição para o regime supersônico, e a equação mais simples que
descreve um escoamento nessas condições já possui uma não-linearidade. Para resolver esses problemas, a classe de
métodos mais comum são os métodos de campo, que fazem a correção para o regime transônico utilizando, no caso de
perfis, integrais numéricas duplas sobre o domínio do escoamento. Isso, além de demandar muito tempo, também traz a
inconveniência de se ter que determinar uma região finita sobre o domínio para se fazer as integrações e, assim,
desprezar o escoamento fora dela.
O objetivo deste trabalho foi encontrar uma maneira de substituir as integrais de área sobre o campo por integrais
de linha sobre a fronteira e, dessa forma, resolver a equação do regime transônico. A primeira tentativa foi um método
que faz essa substituição pela integral sobre a fronteira de uma série de funções meta-harmônicas cujos coeficientes são
obtidos pela derivação sucessiva do campo potencial também sobre a fronteira. Ele funcionou bem para problemas de
transferência de calor com geração, normalmente mais bem-comportados, mas falhou no problema aerodinâmico
devido à presença de singularidades nos pontos de estagnação decorrentes da aproximação para pequenas perturbações,
que falha nesses pontos. Elas produziram uma divergência do método porque as derivadas sucessivas iam-se tornando
cada vez maiores com o aumento do número de termos da série, cujos coeficientes, perto dessas singularidades,
cresciam mais rápido que o decaimento dos termos independentes, o que fez com que a própria série fosse
matematicamente divergente.
Então, na segunda parte deste trabalho, foi estudado outro método, uma adaptação do chamado Método da
Reciprocidade Dual, que foi bem-sucedido por não necessitar de derivações sucessivas, mas necessita de valores do
potencial em pontos fora da fronteira até uma altura de cinco vezes a espessura no extradorso e no intradorso. A
vantagem desse método em relação aos métodos de campo é que nele, não são necessárias integrais sobre o campo,
apenas o valor do termo não-linear da equação em pontos isolados. Ele procura aproximar o campo potencial por uma
combinação linear de funções radiais pré-definidas que satisfaçam a equação que se deseja resolver, cada uma centrada
em um dos pontos de uma malha que envolve toda a extensão da corda do perfil e a já mencionada distância de cinco
vezes a espessura nas duas faces do perfil.
Anais do XII ENCITA 2006, ITA, Outubro, 16-19, 2006
,
O problema resolvido pela implementação do método é o escoamento sem sustentação sobre um perfil parabólico
simétrico. Entretanto, é possível, utilizando a forma completa das equações, resolver problemas com sustentação. Já a
representação numérica de ondas de choque, freqüentes em escoamentos transônicos, é mais complexa e não foi feita
neste trabalho. Quando o perfil e o escoamento têm condições que causam esse fenômeno, o método diverge.
2. Equação do regime transônico e transformação afim
Com as hipóteses de escoamento não-viscoso e irrotacional, pode-se definir um potencial de velocidade φ tal que:
r
V = ∇φ
(1)
Além disso, define-se a velocidade de perturbação, que, bidimensionalmente, é dada por:
r r
Vˆ = (uˆ , vˆ ) = V − V∞
(2)
Conseqüentemente, o potencial de perturbação é definido como:
∇φˆ = Vˆ
(3)
A equação que descreve o potencial de perturbação no regime compressível bidimensional, se a direção do
escoamento for considerada como sendo o eixo x, é:
γ + 1 ˆ γ + 1 ˆ 2 γ −1 ˆ 2  ˆ
φ x + 2 φ x + 2 φ y φ xx +
+ φˆ yy = M ∞2 
2V ∞
2V ∞
 V∞

ˆ
 φ y  φˆ 
 γ −1 ˆ γ + 1 ˆ 2 γ −1 ˆ 2  ˆ
M ∞2 
φ x + 2 φ y + 2 φ x φ yy + M ∞2  1 + x  φˆxy + φˆyx
2V∞
2V∞
 V∞  V∞ 
 V∞

(1 − M )φˆ
2
∞
xx
(



(4)
)
onde M∞ é o número de Mach do escoamento não-perturbado e γ, a razão de calores específicos, CP/CV.
A Equação (4) é válida para qualquer escoamento não-viscoso e irrotacional. Entretanto, ela é extremamente nãolinear e, por isso, de difícil solução. Se for feita uma aproximação para pequenas perturbações, para números de Mach
fora das faixas transônica e hipersônica, essa equação se transforma em:
(1 − M )φˆ
2
∞
xx
+ φˆ yy = 0
(5)
A solução da Eq. (5) para números de Mach abaixo da faixa transônica pode ser encontrada a partir da equação de
Laplace, do regime incompressível:
∇ 2φ = 0
(6)
que também é válida para o potencial de perturbação, por meio de uma transformação afim, denominada transformação
de Prandtl-Glauert. No regime supersônico, a Equação (5) muda de elíptica para hiperbólica e sua solução é encontrada
usando uma transformação afim em:
φ yy − φ xx = 0
(7)
Devido a essa mudança de tipo, a Equação (5) não descreve adequadamente o regime transônico porque nele há
uma coexistência de regiões subsônica e supersônica sobre o perfil. Como essa transição não se dá de forma abrupta, o
primeiro dos termos não-lineares da Eq. (4) se torna importante quando o número de Mach vai se aproximando de 1, e a
equação mais simples para a faixa transônica é:
(1 − M )φˆ
2
∞
xx
M 2 (γ + 1) ˆ ˆ
+ φˆ yy = ∞
φ x φ xx
V∞
(8)
Anais do XII ENCITA 2006, ITA, Outubro, 16-19, 2006
,
Mediante a mudança de variáveis:
β 2 = 1 − M ∞2
K=
(9)
2
∞
M (γ + 1)
V∞
(10)
~
x=x
~
y = βy
~
φ =
(11)
(12)
Kφˆ
(13)
β2
A Equação (8) se torna:
~
~
~ ~
(14)
φ ~x ~x + φ ~y~y = φ ~x φ ~x ~x
que é a forma canônica da equação do regime transônico, para a qual serão descritos os dois métodos de solução
mencionados na introdução.
Nas implementações dos métodos, é resolvido o escoamento sobre um perfil parabólico e são pedidos ao usuário a
espessura relativa do perfil e o número de Mach. A saída é a distribuição de coeficientes de pressão (CP) sobre o perfil,
que, no caso sem sustentação, é igual no extradorso e no intradorso. O CP linearizado é dado por:
C P = −2
uˆ
V∞
(15)
Observando a equação (14), nota-se que a mudança de tipo de equação (ou seja, número de Mach igual a 1) ocorre
quando a velocidade horizontal de perturbação no plano transformado atinge o valor 1, que equivale a β2/K no plano
físico. Portanto, substituindo essa expressão na Eq. (15), tem-se:
C Pcrit = −2
β2
(16)
KV∞
3. Método da série meta-harmônica
3.1. Dedução da fórmula de integração sobre a fronteira
Seja uma equação qualquer da forma:
∇ 2φ + R(φ , P) = 0
(17)
O Teorema de Green, usado para transformar integrais de domínio em integrais de fronteira, é dado, em sua forma
bidimensional, por:
∫∫ (φ1∇ φ2 − φ2∇ φ1 )dS = ∫ (φ1∇φ2 − φ2∇φ1 ) ⋅ nˆ ds
2
2
S (C )
(18)
C(S )
onde φ1 e φ2 são campos escalares quaisquer. Já que:
∇φ ⋅ nˆ =
∂φ
∂n
(19)
Tem-se:
∫∫ (φ1∇ φ2 − φ2∇ φ1 )dS = ∫
2
S (C )
2
C(S )
((φ1
∂φ2
∂φ
− φ2 1 )ds
∂n
∂n
(20)
Anais do XII ENCITA 2006, ITA, Outubro, 16-19, 2006
,
Define-se:
∇ 2 Gi ( P, Q) = Gi −1 ( P, Q)
(21)
∇ 2 G 0 ( P, Q ) = δ ( P , Q )
onde δ(P, Q) é a função Delta de Dirac definida em P e centrada em Q.
[
Multiplicando a Eq. (17) por G0(P,Q), obtém-se:
G 0 ( P, Q)∇ 2φ ( P) + R (φ , P) G 0 ( P, Q) = 0
(22)
E ainda, a integração da Eq. (22) na variável P para Q constante resulta em (em duas dimensões):
2
∫∫ G0 ( P, Q)∇ φ ( P) dS + ∫∫ R(φ , P) G0 ( P, Q) dS
S
(23)
S
Para aplicar o Teorema de Green ao primeiro termo da equação (23), precisa-se calcular o valor de
2
∫∫φ ( P)∇ G0 ( P, Q) dS .
S
∇ 2 G0 ( P, Q ) só é diferente de zero em Q:
2
2
∫∫φ ( P)∇ G0 ( P, Q) dS = φ (Q) ∫∫ ∇ G0 ( P, Q) dS
S
(24)
S
∇ 2G0 ( P, Q) = δ ( P, Q ) → ∫∫ ∇ 2G0 ( P, Q ) dS = C
(25)
S
onde C é o valor principal de Cauchy, que depende do ângulo sólido em torno do ponto. C vale 1 nos pontos do interior
do domínio, ½ nos pontos de fronteira suave e, nos pontos onde a fronteira tiver quinas ou vértices, tem o valor da razão
entre o ângulo sólido em torno do ponto e o ângulo sólido máximo (4π).
2
∫∫φ ( P)∇ G0 ( P, Q) dS = Cφ (Q)
(26)
S
Subtraindo (23) de (26), tem-se:
2
2
∫∫ (φ ∇ G0 − G0∇ φ ) dS − ∫∫ ( R G0 ) dS = Cφ (Q)
S
(27)
S
∂G ( P, Q )
∂φ ( P) 

Cφ (Q ) = ∫  φ ( P ) 0
− G0 ( P, Q)
 ds − ∫∫ R (φ , P )G0 ( P, Q ) dS
∂n
∂n 
S
S
(28)
Definindo:
f i = ∇ 2 f i −1
(29)
f 1 = R(φ , P )
Tem-se:
2
∫∫ R(φ , P)G0 ( P, Q) dS = ∫∫ f1 (φ , P)G0 ( P, Q) dS = ∫∫ f1∇ G1 dS =
S
S
S
∂f1 
 ∂G1
2
2
2
∫∫ ( f1∇ G1 − G1∇ f1 ) dS + ∫∫ G1∇ f1 dS = ∫  f1 ∂n − G1 ∂n  dS + ∫∫ f 2G1 dS
S
S
C
S
(30)
Anais do XII ENCITA 2006, ITA, Outubro, 16-19, 2006
,
Para um i arbitrário, obtém-se:
2
2
∫∫ f iGi −1 dS = ∫∫ f i∇ Gi dS = ∫∫ ( fi∇ Gi − Gi∇
S
S
2
f i ) dS + ∫∫ Gi ∇ 2 f i dS =
S
S
∂f i 
 ∂Gi
∫  f i ∂n − Gi ∂n  dS + ∫∫ f i +1Gi dS
C
S
(31)
Utilizando os resultados (30) e (31), por indução, tem-se:

∞
∫∫ R(φ , P)G 0 ( P, Q) dS = ∑
∫  f i (φ , P)
i =1 
S
C
∂Gi ( P, Q)
∂f (φ , P) 
− Gi ( P, Q ) i
 ds
∂n
∂n 
(32)
Substituindo (32) em (28), chega-se à expressão final:
∞
∂G ( P, Q )
∂φ ( P ) 
∂G ( P, Q )
∂f (φ , P) 


− G0 ( P, Q )
− Gi ( P, Q ) i
Cφ (Q ) = ∫  φ ( P) 0
 ds − ∑ ∫  f i (φ , P) i
ds
∂n
∂n 
∂n
∂n 
i =1 S 
C
(33)
A Equação (33) fornece um resultado importante: é possível a resolução de uma equação da forma (17) se forem
conhecidos os valores do campo escalar φ e de todas as suas derivadas apenas sobre a fronteira. A necessidade de se
conhecer esses valores leva a um método iterativo, exceto no caso linear. Porém, esse método já economiza bastante
esforço computacional devido ao fato de não ser preciso fazer uma discretização de todo o campo, mas somente da
fronteira.
3.2. Singularidades bidimensionais
Em duas dimensões, a função G0 que satisfaz a definição (21) é:
G0 (P, Q ) =
ln r
2π
(34)
onde r é a distância entre P e Q. Essa é a solução fonte da equação de Laplace, (6), e sua derivada em qualquer direção é
uma solução dipolo orientada nessa direção. Portanto, o primeiro termo da Eq. (33) é a solução do escoamento
incompressível, e a série infinita seguinte é a correção para o regime transônico, desde que, na Eq. (17), se assuma
R = - φxφxx.
De acordo com a definição (21), cada função Gi é solução da equação ∇ 2 (i +1)φ = 0 , sendo, portanto, uma função
meta-harmônica, que foi denominada metafonte de ordem i e sua derivada, metadipolo de ordem i.
No método das singularidades usado no regime incompressível, que consiste em distribuir soluções elementares
sobre a corda do perfil, as fontes são responsáveis pela espessura do perfil e os dipolos, pela vorticidade, que é a causa
da sustentação. Na correção para o regime transônico, as metafontes complementam o papel das fontes e os metadipolos
complementam o dos dipolos. Como o problema que se tentou resolver por esse método (escoamento sobre um perfil
parabólico simétrico sem ângulo de ataque) não possui sustentação, foram utilizadas apenas fontes e metafontes.
Para encontrar uma fórmula geral de recorrência para as metafontes, será usado o laplaciano em coordenadas
bidimensionais polares:
∇ 2Gi (r ,θ ) =
1 ∂  ∂Gi  1  ∂ 2Gi
r
+ 
r ∂r  ∂r  r 2  ∂θ 2




do qual será utilizado somente o primeiro termo porque as metafontes são apenas funções de r:
1 ∂  ∂G1
r
r ∂r  ∂r

 = ln r

∂  ∂G1 
r
 = r ln r
∂r  ∂r 
(35)
Anais do XII ENCITA 2006, ITA, Outubro, 16-19, 2006
,
r
∂G1 r 2
r2
= ln r −
∂r
2
4
∂G1 r
r
= ln r −
∂r
2
4
2
r
r2 r2 r2
G1 =
ln r −
−
=
(ln r − 1)
4
8
8
4
1 ∂  ∂G i +1
r
r ∂r  ∂r
(36)

 = r 2i (C i ln r − Di )

∂  ∂Gi +1 
r
 = r 2i +1 (C i ln r − Di )


∂r  ∂r 
r
∂Gi +1
r 2(i +1)
r 2(i +1)
r 2(i +1)
=
C i ln r −
C
−
Di
i
2
∂r
2(i + 1)
2(i + 1)
4(i + 1)
∂G i +1
r 2i +1
r 2i +1
r 2i +1
=
C i ln r −
C
−
Di
i
2
∂r
2(i + 1)
2(i + 1)
4(i + 1)
Gi +1 =
Gi +1 =
2
C i  r 2 (i +1)
r 4 (i +1)  r 2(i +1)
r 2(i +1)
ln r −
−
C
−
Di


i
2
3
2
2(i + 1)  2(i + 1)
4(i + 1)  8(i + 1)
4(i + 1)
Ci
4(i + 1)
2
r 2(i +1) ln r −
r 2(i +1)
4(i + 1)
3
Ci −
r 2(i +1)
4(i + 1)
2
Di
 Ci
Ci
Di 
Gi +1 = r 2(i +1) 
ln r −
−
2
3
2 
4(i + 1)
4(i + 1) 
 4(i + 1)
(37)
Usando os resultados (36) e (37), por indução, chega-se à forma recursiva geral:
1 2i
r (C i ln r − Di )
2π
C0 = 1
Gi =
D0 = 0
C i +1 =
Di +1 =
(38)
Ci
4(i + 1)
Ci
2
3
4(i + 1)
+
Di
4(i + 1)
2
3.3. Implementação numérica
Em cada iteração, o programa vai calculando a próxima distribuição de potencial sobre a corda do perfil utilizando
a Eq. (33), tendo como estimativa inicial a solução do regime incompressível pelo método das singularidades. Quando
há a convergência do potencial, há como, usando a definição (29) (para R = - φxφxx), encontrar os coeficientes da série
infinita e, usando a Eq. (33), encontrar o potencial em qualquer ponto do campo. Na realidade, a variável iterada é a
velocidade em x, φx, e as derivadas são encontradas pelo método das diferenças finitas centradas na maioria dos painéis
e adiantada e atrasada no primeiro e no último painel. Além disso, são feitos cálculos também nos pontos
imediatamente acima e abaixo da corda, para que os laplacianos seja calculados numericamente usando uma função do
software MATLAB, que foi usado na implementação.
Como já mencionado, o método teve bons resultados em problemas de transferência de calor, mas divergiu no
problema aerodinâmico devido às singularidades nos pontos de estagnação decorrentes da falha nesses pontos da
hipótese de pequenas perturbações. Analisando a forma recursiva (38), nota-se, em cada metafonte, os termos que
multiplicam a forma geral têm um decaimento da ordem de (2i)! em relação à primeira. Se as derivadas sucessivas
crescerem mais rápido que esse decaimento, que é o que ocorre perto de uma singularidade, a série diverge, o que
Anais do XII ENCITA 2006, ITA, Outubro, 16-19, 2006
,
ocorre também com o método iterativo. Devido a isso, usou-se um método alternativo, uma variação do Método da
Reciprocidade Dual, que não elimina a singularidade, mas converge por não serem necessárias derivações sucessivas.
4. Método da Reciprocidade Dual
4.1. Fórmula de integração sobre a fronteira
Seja uma equação da forma:
∇ 2φ = σ (φ , P )
(39)
que é idêntica à forma (17), mas com o termo não-linear no segundo membro da equação por conveniência do método.
Aplicando o Teorema de Green apenas ao laplaciano, da mesma forma que no item 2.1, a forma integral da equação
(39) é:
∂φ
∂G
G dΓ + ∫ φ dΓ + ∫∫ σG dΩ
∂
n
Γ
Γ ∂n
Ω
cφ ( P ) = − ∫
(40)
onde Γ é a fronteira do domínio, Ω é o domínio e G é a função G0 do outro método, a solução fonte da equação de
Laplace. É essa a integral utilizada nos métodos convencionais, onde a última integral é uma integral de área, que será
eliminada neste método.
O método da Reciprocidade Dual utiliza um conceito análogo ao do método de solução de equações diferenciais
ordinárias lineares não-homogêneas, que resolve a equação homogênea e soma ao resultado uma solução particular da
equação não homogênea:
(41)
φ = φh + φ p
que satisfazem:
∇ 2φ h = 0
(42)
2
(43)
∇ φp =σ
A equação homogênea é a equação de Laplace, que descreve o regime incompressível. Sua solução pode ser
encontrada utilizando qualquer método de resolução desse regime. Já a solução particular é considerada como sendo
composta de uma série de funções φj:
∇ 2φ p = σ = ∑ α j ∇ 2φˆ j
(44)
j
O termo de compressibilidade, σ, pode ser escrito como uma série de funções pré-definidas fj:
σ = ∑α j f j
(45)
j
E os termos da solução particular são dados por:
∇ 2φˆj = f j
(46)
Com isso, fazendo uma expansão pelo Teorema de Green análoga à da equação original, tem-se:
cφˆ j (x i ) = − ∫
∂φˆ j
∂G
GdΓ + ∫ φˆ j
dΓ + ∫∫ f j GdΩ
∂n
∂n
(47)
Anais do XII ENCITA 2006, ITA, Outubro, 16-19, 2006
,
Logo:

α j cφˆ j ( xi ) = α j  − ∫


∂φˆ j

∂G
GdΓ + ∫ φˆ j
dΓ  + ∫∫ α j f j GdΩ

∂n
∂n

(48)
E, como:
σ = ∑α j f j
(49)
j
Tem-se:

∂φˆ j

∂G
GdΓ + ∫ φˆ j
dΓ  + ∫∫ σGdΩ

∂n
∂n
j
j




∂φˆ j
∂G
→ ∫∫ σGdΩ = ∑ α j  cφˆ j (Pi ) + ∫
GdΓ − ∫ φˆ j
dΓ 


∂n
∂n
j


∑ α j cφˆ j (xi ) = ∑ α j  − ∫
(50)
Portanto, a equação para φ em um ponto i arbitrário se torna:

∂φˆj
∂φ
∂G
∂G 
G dΓ + ∫ φ dΓ + ∑ α j  cφˆij + ∫
GdΓ − ∫ φˆ j
dΓ


∂n
j
Γ ∂n
Γ ∂n
Γ ∂n
Γ


cφi = − ∫
(51)
As duas primeiras integrais constituem a solução da equação homogênea, que é a equação de Laplace. Já o termo
seguinte é a correção de compressibilidade, devido ao termo σ, que, na equação do regime transônico, corresponde ao
termo não-linear φxφxx. Além disso, como já dito no outro método, os termos em G correspondem a distribuições de
fontes sobre a fronteira, que fornecem a espessura do perfil, e os termos em ∂G/∂n correspondem a dipolos,
responsáveis pela vorticidade e, conseqüentemente, pela sustentação. Como só foi simulado o caso sem sustentação, os
termos em ∂G/∂n foram descartados.
A solução do regime incompressível foi encontrada distribuindo fontes sobre a corda. Após a linearização das
condições de contorno, a densidade linear de fonte q(x) é dada por:
q( x ) = 2V∞
dzt
dx
(52)
onde zt é a linha da semi-espessura do perfil.
Utilizando esse tipo de solução, com distribuição de fontes sobre a corda (e não sobre a fronteira), as condições de
contorno para o Método da Reciprocidade Dual passam a ser aplicadas sobre a corda. Com isso, no caso sem
sustentação, a integral de linha desaparece porque extradorso e intradorso coincidem. Assim, a equação (51) se
transforma em:
∂φ
∂G
G dΓ + ∫ φ dΓ + ∑ α j cφˆij
j
Γ ∂n
Γ ∂n
cφi = − ∫
(53)
4.2. Determinação do termo de compressibilidade
No Método da Reciprocidade Dual, as funções fj dependem de r (P, Pj), que é a distância ||P - Pj|| do ponto onde se
quer calcular a função a um ponto Pj escolhido para cada fj. Nesta implementação, foram escolhidas funções fj = 1 + r.
Com isso, aplicando inversamente o laplaciano em coordenadas polares bidimensionais (Eq. (35)), tem-se:
φˆ j =
r2 r3
+
4
9
(54)
Anais do XII ENCITA 2006, ITA, Outubro, 16-19, 2006
,
No método, é necessário que se escolham pontos Pj também fora da fronteira. Nesta implementação, foram
escolhidas 10 camadas de pontos igualmente espaçados na direção x, e o número de painéis na corda é determinado
pelo usuário. Na direção y, é aproveitada a simetria (distribuição de camadas apenas no extradorso), e feito um
refinamento quadrático nas camadas próximas à corda, com a última camada a uma distância da corda cinco vezes
maior que a espessura relativa. Na direção x, são tomados pontos de uma distância de dois terços da corda à frente do
bordo de ataque até a mesma distância atrás do bordo de fuga.
Para encontrar os coeficientes αj, é resolvido o seguinte sistema linear:
{σ } = [F ]{α }
(55)
onde Fij = fj (Pi) e σi = φxφxx calculado no ponto Pi.
O método foi implementado utilizando o software Matlab®. A variável iterada é a velocidade induzida em x, u =
φx. A derivada φxx é encontrada usando métodos de diferenças finitas adiantadas, centradas e atrasadas de segunda
ordem nos dois pontos próximos a cada extremo horizontal (que estão fora da corda) e um método de diferenças finitas
centradas de quarta ordem nos pontos centrais.
4.3. Algoritmo numérico
1)
2)
3)
4)
5)
Resolver o campo da equação linear, usando a distribuição de fontes da equação (52)
Calcular em todos os pontos da malha o termo de compressibilidade, φxφxx
Encontrar os coeficientes αj, resolvendo o sistema (55)
Calcular o novo campo de velocidades usando a equação (53)
Retornar ao passo 2 e repetir os passos seguintes até atingir a convergência
4.4. Resultados
As simulações foram feitas com 20 painéis na corda, e a malha, com espaçamento na direção x constante e igual a
c/20. A teoria linear sofreu a correção de Prandtl-Glauert.
Figura 1. Distribuições de coeficiente de pressão para M = 0.8 e t/c = 6%.
Anais do XII ENCITA 2006, ITA, Outubro, 16-19, 2006
,
Figura 2. Distribuições de coeficiente de pressão para M = 0.83 e t/c = 6%.
Figura 3. Distribuições de coeficiente de pressão para M = 0.84 e t/c = 6%.
Anais do XII ENCITA 2006, ITA, Outubro, 16-19, 2006
,
De acordo com o presente método, o número de Mach crítico do perfil parabólico simétrico de espessura 6% com
ângulo de ataque nulo está entre 0.83 e 0.84. Como já mencionado, a representação numérica de ondas de choque é
complexa e não foi feita nesta implementação, na qual houve divergência para números de Mach 0.85 e mais altos.
Observando os gráficos, nota-se que já para um número de Mach igual a 0.8 há uma diferença significativa entre as
resoluções da equação potencial completamente linearizada (Eq. (5)) e da equação do potencial do regime transônico
(Eq. (8)). Essa última equação pode ser escrita como:
(1 − β
2
)
− Kφ x φ xx + φ yy = 0
(56)
Essa forma mostra que, no regime transônico, há um acréscimo de compressibilidade em relação à equação linear,
o que ocorreu nos resultados. Devido à disparidade entre os resultados, conclui-se que, para cálculos de parâmetros
aerodinâmicos ligados à faixa transônica, como, por exemplo, o número de Mach crítico, é necessário considerar a nãolinearidade.
Os resultados do presente método estão muito próximos dos obtidos a partir de integrais duplas de campo, que
conferem com a experiência. Além disso, o presente método leva vantagem por demandar menor esforço computacional
devido à eliminação das integrais duplas.
5. Conclusão
Foram aqui apresentados dois métodos que procuram transformar integrais de área em integrais de linha. O
primeiro deles, se houvesse sido bem-sucedido para aplicações aerodinâmicas, seria o de menor esforço computacional
por ser totalmente uma formulação apenas sobre a fronteira. Já o segundo método substitui as integrais de campo por
cálculo algébrico por meio da aproximação do campo por funções pré-definidas das quais só seria necessário encontrar
os coeficientes. Esse cálculo necessita de valores fora da fronteira, logo, é um método intermediário entre o primeiro
mostrado e os métodos convencionais que calculam integrais duplas. Por ser mais rápido que os métodos convencionais
e ter produzido bons resultados (ao contrário do primeiro), o Método da Reciprocidade Dual se mostrou a melhor entre
as soluções estudadas para resolver o problema aerodinâmico citado.
6. Agradecimentos
Agradecemos ao CNPq pelo apoio e financiamento deste projeto de iniciação científica
7. Referências
Anderson Jr., J.D., 1984, “Fundamentals of Aerodynamics”, McGraw-Hill
Hess, J. L., Smith, A. M. O., “Calculation of Potential Flow About Arbitrary Bodies”, Progress in Aeronautical
Sciences, 1967.
Hildebrand, F.B., 1948, “Advanced Calculus for Applications”, Prentice-Hall
Partridge, P.W., Brebbia, C.A. and Wrobel, L. C., 1992 “The Dual Reciprocity Boundary Element Method”,
Computational Mechanics Publications.
Ribeiro, R. S. and Soviero, P. A. O., 1987, “Análise do Escoamento Transônico Através do Método das
Singularidades”, Master Thesis, ITA.
Ribeiro, R. S. and Soviero, P. A. O., 1987, “Calculation of Transonic Flow About Airfoils by a Field Panel Method”,
Boundary Element Techniques: Applications in Fluid Flow and Computational Aspects, Computational Mechanics
Publications.
Schilichting, H., Truckenbrodt, E., 1979, “Aerodynamics of the airplane”, McGraw-Hill, New York
Shapiro, A.H., 1953, “The dynamics and termodynamics of compressible flow”, The Ronald Press, New York
Silva, B.G.O., 2001, “Estudo de perfis aerodinâmicos em alta velocidade”, Trabalho de Graduação, ITA
Spreiter, J. R. and Alksne, A., 1955, “Theoretical prediction of pressure distributions on nonlifting airfoils at high
subsonic speeds”, NACA Report 1217.
Uhl, B., Ostertag, J., Guidati, G. and Wagner, S., 1999, “Application of the Dual-Reciprocity Method to ThreeDimensional Compressible Flows Governed by the Full-Potential Equation”, 37th AIAA Aerospace Sciences
Meeting and Exhibit, January 11-14, 1999, Reno, NV.
Zuosheng, Y., 1988, “A complete boundary integral formulation for steady compressible inviscid flows governed by
non-linear equations”, Int. Journal for Numerical Methods in Fluids, Vol. 16, pp. 231-237