Tratamento Numérico de Escoamentos Aerodinâmicos'breno.freeshell.org/Aer_num.pdf · 2009. 9. 20. · Breno Moura Castro Tratamento Numérico de Escoamentos Aerodinâmicos SÃO

Breno Moura Castro

Tratamento Numérico de EscoamentosAerodinâmicos

SÃO JOSÉ DOSCAMPOS - 2009

2

Sumário

1 Introdução 51.1 Objetivos gerais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .51.2 Objetivos específicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .51.3 Ementa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.4 Conteúdo programático . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .5

2 Noções básicas 72.1 Números de Mach e de Reynolds . . . . . . . . . . . . . . . . . . . . . . . . . . . .72.2 Classificação dos escoamentos segundo o Número de Mach . . . . . . . . . . . . . .82.3 Equações governantes do escoamento . . . . . . . . . . . . . . . . . . . . . . . . .82.4 Aproximações das equações governantes . . . . . . . . . . . . . . . . . . . . . . . .9

2.4.1 Escoamentos incompressíveis e não-viscosos . . . . . . . . . . . . . . . . .92.4.2 Pequenas perturbações . . . . . . . . . . . . . . . . . . . . . . . . . . . . .11

2.5 Coeficiente de Pressão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .132.6 Expansão em Série de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .142.7 Elementos de vibrações livres e forçadas . . . . . . . . . . . . . . . . . . . . . . . .14

2.7.1 Movimento harmônico simples . . . . . . . . . . . . . . . . . . . . . . . . .14

3 Métodos básicos de discretização 173.1 Métodos de Singularidades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .17

3.1.1 Escoamento potencial sobre uma placa plana . . . . . . . . . . . . . . . . .173.1.2 Escoamento potencial sobre uma asa finita . . . . . . . . . . . . . . . . . .20

3.2 Diferenças finitas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .233.3 Volumes finitos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .263.4 Estabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .27

4 Elementos de Aeroelasticidade 294.1 Conceito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .294.2 Aeroelasticidade estática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .31

4.2.1 Divergência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .314.2.2 Reversão de Controle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .32

5 Experimentos com estruturas 355.1 Extensometria elétrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .35

5.1.1 Conceitos básicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .35

A Programa em FORTRAN para cálculo de uma Placa Plana 39

B Programa em FORTRAN para o Método “Vortex-Lattice” 43

3

4 SUMÁRIO

Capítulo 1

Introdução

Este texto tem a finalidade de auxiliar o aluno da disciplina Tratamento Numérico de EscoamentosAerodinâmicos, na modalidade de ensino presencial, a passar pelo processo de aprendizagem domaterial proposto.

1.1 Objetivos gerais

Proporcionar ao aluno a utilização e compreensão de ferramentas de simulação numérica paraescoamentos de fluidos, compressíveis e incompressíveis,em malhas, permitindo que o aluno se con-centre nos conceitos e nas técnicas numéricas básicas.

1.2 Objetivos específicos

Estudar o tratamento numérico aplicado ao escoamento de fluidos, proporcionando a necessáriaexperiência para o tratamento de simulações elaboradas.

1.3 Ementa

Análise de aerofólios em escoamentos bi e tridimensionais incompressíveis e compressíveis. In-trodução à Aerodinâmica Não-Estacionária. Fenômeno de Flutter.

1.4 Conteúdo programático

Introdução à Dinâmica dos Fluidos Computacional, Equações de Diferenças Finitas, Técnicas deSoluções Numéricas, Discretização Uni e Multidimensional, Métodos Numéricos para Navier-Stokes,Modelamento em software para simulação.

5

6 CAPÍTULO 1. INTRODUÇÃO

Capítulo 2

Noções básicas

Este capítulo trata de alguns dos principais conceitos relacionados à Aerodinâmica. Dentre elesdestacam-se as definições dos Números de Mach e de Reynolds, parâmetros importantes para a com-partimentalização do estudo desta disciplina, e das principais equações para o estudo de escoamentosviscosos e não-viscosos.

2.1 Números de Mach e de Reynolds

O som é uma onda mecânica que se propaga num gás, líquido ou sólido. A velocidade com queesta onda se propaga depende de características elásticas e inerciais do meio. No ar, esta velocidadedepende fundamentalmente da temperatura do fluido, de acordo com a Eq. (2.1).

a =√

γRT (2.1)

em quea é a velocidade do som,γ é a razão de calores específicos do ar (γ = 1, 4), R é a constantedos gases para o ar (R = 287 J/kg.K) eT é a temperatura absoluta do ar.

Associada à velocidade do som, utiliza-se a definição de um número adimensional conhecidocomo Número de Mach, em homenagem a Ernest Mach. Este é um dos parâmetros mais relevantesno estudo da Aerodinâmica e um dos principais em Aeronáutica. A definição de Número de Mach éfeita de acordo com a Eq. (2.2).

M =V

a(2.2)

na qualM corresponde ao Número de Mach,V representa a velocidade do escoamento ea é avelocidade do som.

Outro parâmetro de especial relevância para a Aerodinâmica é o Número de Reynolds, que repre-senta a relação entre as forças de inécia e as forças viscosas no escoamento. A definição do Númerode Reynolds encontra-se representada na Eq. (2.3).

Re =ρ∞V∞lref

µ∞(2.3)

na qualρ∞ é a densidade,V∞ representa a velocidade eµ∞ a viscosidade, todas relativas ao escoa-mento não-perturbado, enquanto quelref simboliza um comprimento característico do corpo.

Estes dois números são referências para as definições das características de um escoamento e sãoutilizados, então, para auxiliar na classificação de seus tipos. O Número de Mach, por exemplo, servepara diferenciar escoamentos compressíveis daqueles que podem ser considerados incompressíveis.O Número de Reynolds indica aqueles escoamentos em que a viscosidade tem um papel importante,diferenciando escoamentos que podem ser tratados como viscosos ou não-viscosos.

7

8 CAPÍTULO 2. NOÇÕES BÁSICAS

2.2 Classificação dos escoamentos segundo o Número de Mach

Os escoamentos de ar em torno de um corpo podem ser classificados segundo a velocidade. Noentanto, deve-se levar em consideração que essa velocidade não é a mesma em todos os pontos doescoamento. A presença do corpo altera a direção e a intensidade da velocidade de algumas partículasdo fluido, fazendo com que ocorra uma variação quanto à velocidade das partículas que ainda estãodistantes do corpo.

Desta forma, para a classificação dos tipos de escoamento, toma-se por base a velocidade ditado escoamento não-perturbado, ou seja, a velocidade das partículas que estão longe o suficiente dainfluência do corpo sobre o escoamento. Esta velocidade é referenciada comoV∞.

O Número de Mach do escoamento não-perturbado é definido tomando-se como referência avelocidade do som na região afastada do corpo. A Eq. (2.4) traduz matematicamente este parâmetro.

M∞ =V∞a∞

(2.4)

Os regimes de vôo podem ser classificados de acordo com o Número de Mach associado aoescoamento não-perturbado,M∞. A divisão mais comum é apresentada a seguir, juntamente comalgumas das características mais importantes de cada regime.

• Subsônico:0 ≤ M∞ ≤ 0, 8 - ausência de ondas de choque; as perturbações geradas pelapresença do corpo são percebidas em todas as direções.

• Transônico:0, 8 ≤ M∞ ≤ 1, 2 - início da formação de ondas de choque; apresenta regiões emque o escoamento é subsônico e outras onde as velocidades são supersônicas.

• Supersônico:1, 2 ≤ M∞ ≤ 5, 0 - presença de ondas de choque oblíquas; as perturbaçõesgeradas em um ponto somente são percebidas a jusante desse ponto.

• Hipersônico:M∞ ≥ 5, 0 - ondas de choque oblíquas muito próximas à superfície do corpo;surgimento de dissociações das moléculas que constituem o fluido; transferência de calor in-tensa e presença de íons na região da onda de choque.

O Número de Mach do escoamento não-perturbado para o qual surge o primeiro ponto de veloci-dade local do escoamento igual a um (velocidade local sônica) recebe o nome deMach Crítico. Estenúmero de Mach depende da forma do corpo; em especial, depende da espessura do corpo medida nadireção perpendicular à do escoamento não-perturbado.

2.3 Equações governantes do escoamento

As equações fundamentais que descrevem a dinâmica do escoamento de um fluido, na notação deEinstein, são:

1. Equação da Continuidade:∂ρ

∂t+

∂(ρuj)

∂xj

= 0 (2.5)

2. Equação do Momentum:

∂

∂t(ρui) +

∂

∂xj

(Πij + ρuiuj) = 0 Πij = pδij + τij (2.6)

2.4. APROXIMAÇÕES DAS EQUAÇÕES GOVERNANTES 9

3. Equação da Energia:

∂e

∂t+

∂

∂xi

(eui + Πijuj + qi) = 0 e =1

2ρuiui + ρcvT qi = −k

∂T

∂xi

(2.7)

nas quaisρ é a densidade do fluido,p corresponde à pressão,ui representa as componentes do vetorvelocidade do fluido et o tempo.τij simboliza o tensor de tensões no fluido.

A Equação do Gás Perfeito é utilizada para equalizar o número de equações e incógnitas dosistema:

p = ρRT (2.8)

Em alguns casos, é necessário avaliarDui/Dt, que é a componente da aceleração de uma partículado fluido. Para expressar este parâmetro, é necessário conhecer as componentes da velocidade dofluido:

ui =dxi

dt= ui(x1, x2, x3, t)

que, como pode-se observar, são funções do ponto que a partícula de fluido ocupa no espaço e tambémdo tempo. A regra de diferenciação fornece:

Dui

Dt=

∂ui

∂t+

∂ui

∂xj

∂xj

∂t=

∂ui

∂t+ uj

∂ui

∂xj

(2.9)

Por sua vez, o tensor de tensões pode ser escrito como:

τij = µ

[(∂ui

∂xj

+∂uj

∂xi

)− 2

3δij

∂uk

∂xk

]i, j, k = 1, 2, 3 (2.10)

Substituindo as Eqs. (2.9) e (2.10) na Eq. (2.6):

∂ui

∂t+ uj

∂ui

∂xj

= Fi −1

ρ

∂p

∂xi

+1

ρ

∂τ ij

∂xj

(2.11)

2.4 Aproximações das equações governantes

2.4.1 Escoamentos incompressíveis e não-viscosos

Uma aproximação comum é obtida por meio da hipótese de que a densidade do fluido é constante,ou seja, não varia nem com o tempo e nem com a posição dentro do escoamento (ρ = constante).Esta hipótese é conhecida simplesmente como escoamento incompressível. Assim, a Equação daContinuidade (2.5) pode ser escrita na forma:

∂ui

∂xi

= 0 (2.12)

Da mesma forma, se for assumido que não existem forças de volume atuando no fluido (Fi = 0)e escoamento permanente (∂(.)/∂t = 0), a Equação da Conservação da Quantidade de Movimentotorna-se:

uj∂ui

∂xj

= −1

ρ

∂p

∂xi

(2.13)

A vorticidade do escoamento é definida como sendo o rotacional do vetor velocidade e é dada por:

∇× u =

(∂u3

∂x2

− ∂u2

∂x3

)e1 +

(∂u1

∂x3

− ∂u3

∂x1

)e2 +

(∂u2

∂x1

− ∂u1

∂x2

)e3 (2.14)


em quee1, e2 ee3 representam os versores unitários nas direçõesx1, x2 ex3, respectivamente.Se for admitido, ainda, que o escoamento é irrotacional (∇× u = 0), pode-se utilizar a seguinte

identidade matemática para simplificar as equações obtidas para escoamento incompressível:

∇×∇Φ ≡ 0 (2.15)

em queΦ representa uma função escalar das coordenadas (x1,x2,x3).Por inspeção, observa-se que, para escoamentos irrotacionais, o vetor velocidadeu pode ser es-

crito como:

u = ∇Φ =∂Φ

∂x1

e1 +∂Φ

∂x2

e2 +∂Φ

∂x3

e3 (2.16)

Também por inspeção, pode-se constatar que as componentes do vetor velocidadeu são dadaspor:

ui =∂Φ

∂xi

(2.17)

Substituindo a Eq. (2.17) em (2.12) obtém-se a Equação de Laplace:

∇2Φ = 0 ou∂2Φ

∂x21

+∂2Φ

∂x22

+∂2Φ

∂x23

= 0 (2.18)

Lembra-se que a Eq. (2.18) é válida somente quando:

1. o fluido é não-viscoso;

2. o escoamento é incompressível (ρ = constante); e

3. o escoamento é irrotacional (∇× u = 0).

Por sua vez, a Equação da Conservação da Quantidade de Movimento, dentro das hipóteses deescoamento permanente, incompressível e sem forças de volume para um fluido não-viscoso, podeser reduzida à seguinte equação:

p +1

2ρV 2 = constante (2.19)

na qualp é a pressão,ρ a densidade eV representa o módulo do vetor velocidadeu. Cabe salientarque a Eq. (2.19) é válida somente ao longo de uma linha de corrente do escoamento. No entanto, se oescoamento for irrotacional, a Eq. (2.19) é válida em todo o domínio.

V =√

u21 + u2

2 + u23

Sintetizando, a Eq. (2.19) é válida em todo o escoamento nas seguites condições:

1. o fluido é não-viscoso;

2. o escoamento é incompressível (ρ = constante);

3. o escoamento é permanente (∂(.)/∂t = 0);

4. não existem forças de volume (Fi = 0); e

5. o escoamento é irrotacional (∇× u = 0).

2.4. APROXIMAÇÕES DAS EQUAÇÕES GOVERNANTES 11

2.4.2 Pequenas perturbações

A hipótese de pequenas perturbações consiste de:

1. o escoamento é irrotacional (∇× u = 0);

2. o escoamento não-perturbado encontra-se na direçãox1; e

3. o corpo imserso no escoamento é esbelto, ou seja, sua maior dimensão nas direçõesx2 e x3 émuito menor que o seu maior comprimento na direçãox1.

Desta forma, pode-se admitir que:

u1 = U∞ + u1 u2 = u2 u3 = u3

em queu1, u2, u3 U∞.Da hipótese de escoamento irrotacional, da mesma forma como no caso de escoamentos incom-

pressíveis, pode-se escrever a velocidade como sendo o gradiente de um campo potencial:

u = ∇Φ

No caso de pequenas perturbações, o potencial de velocidadesΦ pode ser escrito da forma adiante:

Φ = φ + U∞x1 (2.20)

Assim, por comparação e utilizando a Eq. (2.16), pode-se escrever:

u1 =∂φ

∂x1

u2 =∂φ

∂x2

u3 =∂φ

∂x3

(2.21)

A Equação da Continuidade para escoamento compressível e permanente pode ser escrita como:

∂(ρui)

∂xi

= 0 (2.22)

Nas mesmas condições de escoamento e sem forças de volume, a Equação de Conservação daQuantidade de Movimento torna-se:

uj∂ui

∂xj

= −1

ρ

∂p

∂xi

(2.23)

A velocidade do som em um gás perfeito é definida como sendo a variação na pressão com relaçãoà variação na densidade do gás em um processo isentrópico (p/ργ = constante). Desta forma:

a2 =

(∂p

∂ρ

)(2.24)

Combinando as Eqs. (2.22), (2.23) e (2.24) e observando que o escoamento é irrotacional, obtém-se:(

1− u21

a2

)∂u1

∂x1

+

(1− u2

2

a2

)∂u2

∂x2

+

(1− u2

3

a2

)∂u3

∂x3

− 2u1u2

a2

∂u1

∂x2

− 2u2u3

a2

∂u2

∂x3

− 2u3u1

a2

∂u3

∂x1

= 0

(2.25)Invocando a hipótese de pequenas perturbações, a Eq. (2.25) torna-se:(

1− u21

a2

)∂u1

∂x1

+∂u2

∂x2

+∂u3

∂x3

= 0 (2.26)


A Equação (2.26) pode ser simplificada utilizando-se uma relação conhecida entre a velocidadedo som local (a) e a velocidade do som no escoamento não-perturbado (a∞) para escoamentos adia-báticos:

a2

a2∞

= 1− γ − 1

2

(u2

1 + u22 + u2

3

U2∞

− 1

)M2∞ (2.27)

Ainda dentro do contexto de pequenas perturbações, pode-se utilizar o Teorema do Binômio paragerar a seguinte aproximação:

a2∞

a2= 1 +

γ − 1

2M2∞

(2

u1

U∞+

u21 + u2

2 + u23

U2∞

)(2.28)

Para fazer simplificações adicionais na Eq. (2.26), é necessário escrever:

1− u21

a2= 1− (U∞ + u1)

2

a2

u1

U∞

U2∞

U2∞

a2∞

a2∞

= 1− (U2∞ + 2u1U∞ + u2

1)

U2∞

M2∞

a2∞

a2(2.29)

Substituindo Eq. (2.28) em (2.29) e abandonando termos de ordem superior:

1− u21

a2≈ 1−M2

∞

[1 +

2u1

U∞

(1 +

γ − 1

2M2∞

)](2.30)

Assim, a Eq. (2.26) pode ser reescrita como:

(1−M2∞)

∂u1

∂x1

+∂u2

∂x2

+∂u3

∂x3

= M2∞

(1 +

γ − 1

2M2∞

)2u1

U∞

∂u1

∂x1

(2.31)

Utilizando as componentes do vetor velocidade dentro da hipótese de pequenas perturbações:

(1−M2∞)

∂u1

∂x1

+∂u2

∂x2

+∂u3

∂x3

=2

U∞

(1 +

γ − 1

2M2∞

)M2∞u1

∂u1

∂x1

(2.32)

A Equação (2.32) representa a Equação de Pequenas Perturbações.Existem casos, excetuando-se os escoamentos na faixa do transônico, que a Eq. (2.32) pode ser

linearizada, ou seja, os termos à direita da equação podem ser negligenciados. Assim, obtém-se aEquação de Pequenas Perturbações Linearizada:

(1−M2∞)

∂u1

∂x1

+∂u2

∂x2

+∂u3

∂x3

= 0 (2.33)

Para escoamentos isentrópicos e, portanto, irrotacionais, pode-se introduzir a função potencialφ,nos moldes da definição apresentada na Eq. (2.21), e a Eq. (2.33) torna-se a Equação do PotencialLinearizada:

(1−M2∞)

∂2φ

∂x21

+∂2φ

∂x22

+∂2φ

∂x23

= 0 (2.34)

Para escoamentos transônicos deve-se reter os termos à direita da Eq. (2.32). Utilizando, tam-bém, o potencialφ para escoamentos irrotacionais, essa equação torna-se a Equação do PotencialTransônica:

(1−M2∞)

∂2φ

∂x21

+∂2φ

∂x22

+∂2φ

∂x23

=2

U∞

(1 +

γ − 1

2M2∞

)M2∞

∂φ

∂x1

∂2φ

∂x21

(2.35)

Novamente, é importante reforçar as condições em que as Equações do Potencial Linearizada eTransônica são válidas:

2.5. COEFICIENTE DE PRESSÃO 13

1. o fluido é um gás perfeito e não-viscoso;

2. o escoamento é isentrópico (p/ργ = constante);

3. o escoamento é permanente (∂(.)/∂t = 0);

4. não existem forças de volume (Fi = 0);

5. o escoamento é irrotacional (∇× u = 0); e

6. pequenas perturbações (u1, u2, u3 U∞).

2.5 Coeficiente de Pressão

O coeficiente de pressão é definido por meio da Eq. (2.36).

Cp ≡p− p∞12ρ∞V 2

∞(2.36)

Lembrando que o ar pode ser considerado como um gás perfeito e que neste caso:

p∞ = ρ∞RT∞ a2∞ = γRT∞ (2.37)

Fazendo a substituição das Eqs. (2.37) em (2.36):

Cp =2

γM2∞

(p

p∞− 1

)(2.38)

Sabe-se que, para escoamentos adiabáticos de um gás termica e caloricamente perfeito:

T +V 2

2cp

= T∞ +U2∞

2cp

(2.39)

na qualcp representa o calor específico a pressão constante do ar.Movimentando os termos e utilizando quecp = γR/(γ − 1):

T − T∞ =U2∞ − V 2

2γR/(γ − 1)(2.40)

Recordando quea∞ =√

γRT∞:

T

T∞− 1 =

γ − 1

2

U2∞ − V 2

γRT∞=

γ − 1

2

U2∞ − V 2

a2∞

(2.41)

Em termos de pequenas perturbações:

V 2 = (U∞ + u1)2 + u2

2 + u23 (2.42)

Substituindo (2.42) em (2.41) e fazendo as simplificações decorrentes:

T

T∞= 1− γ − 1

2a2∞

(2u1U∞ + u21 + u2

2 + u23) (2.43)


É conhecida, também, a relação entre pressão e temperatura para escoamentos isentrópicosp/p∞ =(T/T∞)γ/(γ−1) que, substituída em (2.43) gera:

p

p∞=

[1− γ − 1

2a2∞

(2u1U∞ + u21 + u2

2 + u23)

]γ/(γ−1)

(2.44)

Rearranjando os termos:

p

p∞=

[1− γ − 1

2M2∞

(2u1

U∞+

u21 + u2

2 + u23

U2∞

)]γ/(γ−1)

(2.45)

Aplicando novamente a expansão binomial:

p

p∞≈ 1− γ

2M2∞

(2u1

U∞+

u21 + u2

2 + u23

U2∞

)(2.46)

Substituindo (2.46) em (2.38)

Cp ≈2

γM2∞

[1− γ

2M2∞

(2u1

U∞+

u21 + u2

2 + u23

U2∞

)− 1

](2.47)

Cancelando termos e desprezando termos de ordem superior, chega-se à equação para o coeficientede pressão linearizado:

Cp = −2u1

U∞(2.48)

2.6 Expansão em Série de Taylor

A expansão utilizando Série de Taylor é uma técnica para aproximar o valor de uma função emtorno de um ponto dado. Supondo que se queira aproximar o valor def(x) em torno do pontox0 eque o pontox está a uma distância∆x dex0:

f(x) = f(x0 + ∆x) =∞∑

n=0

1

n!

dnf

dxn(x0)(∆x)n (2.49)

A expansão em Série de Taylor é uma das bases para a transformação de uma equação diferencialem uma equação de diferenças finitas.

2.7 Elementos de vibrações livres e forçadas

2.7.1 Movimento harmônico simples

O movimento harmônico simples caracteriza-se por uma oscilação que se repete no tempo deforma periódica. Uma propriedade importante deste movimento é a freqüência ou número de os-cilações que são completadas por segundo.

O movimento harmônico simples é definido matematicamente pela Eq. (2.50), que é obtida pormeio do modelo da Fig. 2.1.

x(t) = xm cos(ωt + φ) (2.50)

2.7. ELEMENTOS DE VIBRAÇÕES LIVRES E FORÇADAS 15

Figura 2.1: Modelo de um sistema massa-amortecedor-mola.

em quex(t) é a posição ocupada pelo corpo no instantet, xm corresponde à amplitude do movimento,ω simboliza a freqüência eφ é a fase do movimento.

ω =2π

T= 2πf (2.51)

na qualω é a freqüência em rad/s,T corresponde ao período do movimento ef representa a freqüênciaem Hz.

A velocidade do movimento é dada por:

v(t) =dx

dt=

d

dt[xm cos(ωt + φ)]

v(t) = −ωxm sin(ωt + φ) (2.52)

A aceleração, por sua vez, pode ser obtida de:

a(t) =dv

dt=

d

dt[−ωxm sin(ωt + φ)]

a(t) = −ω2xm cos(ωt + φ) = −ω2x (2.53)

A Segunda Lei de Newton estabelece que:

F = ma = −(mω2)x

Assim, se uma força elástica atua no corpo:

F = −kx

Donde se conclui que:k = mω2

n

ou

ωn =

√k

m(2.54)

em que a variávelωn é chamada de freqüência natural do sistema.Supondo o movimento de um conjunto massa+mola+amortecedor forçado por uma forçaF (t):

mx + Bx + kx = F (t)

em quex é a coordenada representando a posição ocupada pelo corpo no instantet, m é a massa docorpo,B é a constante de amortecimento viscoso,k é a constante da mola eF (t) é a função forçante.


Dividindo a equação anterior porm:

x +B

mx +

k

mx =

F (t)

m(2.55)

A Equação (2.55) pode ser reescrita na forma:

x + 2ζωnx + ω2nx = a(t) (2.56)

na qualB/m = 2ζωn, k/m = ω2n ea(t) = F (t)/m.

Uma oscilação em queF (t) = 0 é ditaOscilação Livre, enquanto que paraF (t) 6= 0 é conhecidacomoOscilação Forçada. Para o caso de oscilações livres e sem amortecimento:

x + ω2nx = 0 (2.57)

A solução para a Eq. (2.56) é composta de duas partes: a solução homogênea [a(t) = 0] e asolução particular. A solução da equação homogênea é dada por:

x(t) = xme−ζωnt cos(ωdt)

em queωd = ωn

√1− ζ2.

No caso de oscilação forçada comF (t) = Fm cos(ωt), a solução será:

x(t) = xm cos(ωt + φ) (2.58)

As amplitudes das oscilações dependerão dos valores deω eωn e serão máximas quandoω = ωn.Esta condição é denominadaRessonância.

Exercícios

1. Escreva as definições dos números de Mach e de Reynolds. Especifique, também, como podemser calculadas a velocidade do som e a viscosidade do fluido.

2. Comente sobre a utilidade do conhecimento dos números de Mach e de Reynolds do escoa-mento, ou seja, comente sobre o que eles indicam sobre o escoamento.

3. Classifique os escoamentos segundo o número de Mach, fornecendo a faixa de variação desteparâmetro em cada regime de velocidade.

4. Escreva a Equação de Laplace em coordenadas cartesianas e cite as condições em que estaequação pode ser utilizada.

5. Em termos práticos, ou seja, em relação a parâmetros relativos a características mensuráveis doescoamento (números de Mach e de Reynolds) e características geométricas do corpo imersono escoamento, quais são as limitações para utilização da Equação de Laplace?

6. Forneça a Equação do Potencial Linearizada e descreva as condições em que ela pode ser apli-cada.

7. Em termos práticos, considerando os números de Mach e Reynolds e a geometria do corpo,quais são as restrições para aplicação da Equação do Potencial Linearizada?

8. Repita as dois exercícios anteriores considerando a Equação do Potencial Transônica.

9. Escreva a definição do coeficiente de pressão e forneça a equação do coeficiente de pressãolinearizado. Em que condições esta última equação pode ser utilizada?

Capítulo 3

Métodos básicos de discretização

3.1 Métodos de Singularidades

Este texto não pretende cobrir os métodos baseados em singularidades com profundidade. Umexemplo deste tipo é o Método dos Painéis. Serão vistos, apenas, alguns exemplos de cálculo com autilização de vórtices bidimensionais. Estes exemplos, no entanto, apresentam a essência dos méto-dos de singularidades e servem como uma excelente introdução aos princípios que envolvem estesmétodos.

3.1.1 Escoamento potencial sobre uma placa plana

O estudo do escoamento sobre uma placa plana corresponde ao primeiro passo para o entendi-mento do escoamento em torno de uma asa. Imagine uma asa plana e sem espessura. Além disto,ela possui uma corda finita,c e uma envergadura infinita, de tal forma que o escoamento em tornoda mesma é bidimensional, ou seja, o escoamento em todas as estações ao longo da envergaduraapresentam as mesmas características.

A solução deste escoamento pode ser obtida através da utilização da superposição de soluçõessimples do tipo vórtice e escoamento uniforme. Seja o esquema apresentado na Fig. 3.1.

Um vórtice de intensidadeΓ encontra-se situado a uma distânciac/4 do bordo de ataque do perfil.Um ponto de controle posicionado a3c/4 é definido para imposição da condição de contorno detangência do escoamento à superfície do perfil. Note que precisamos de apenas um ponto de controlepois a única incógnita do problema é a intensidade do vórticeΓ.

A velocidade tangencial induzida por um vórtice de intensidadeΓ sobre um pontoP situado a

6y

-x

----

V∞

-Γ a V∞αVθ

XXXXXzXXXXXy

c

XXXz

XXXy c/2XXzXXyc/4

XXXXXXXXXX

Figura 3.1: Modelamento de uma placa plana.

17

18 CAPÍTULO 3. MÉTODOS BÁSICOS DE DISCRETIZAÇÃO

uma distânciar da origem do vórtice, é dada por:

Vθ =1

r

∂φ

∂θ= − Γ

2πr(3.1)

Existem duas componentes de velocidade atuando no ponto de controle. A primeira é a com-ponente da velocidade do escoamento não-perturbado perpendicular à superfície do perfil que podeser aproximada, para pequenos ângulos de ataque, porV∞α. A segunda corresponde à velocidadeinduzida pelo vórtice no ponto de controle. Esta velocidade também é perpendicular à superfíciedo perfil e depende apenas da distânciac/2 entre o ponto onde está situado o vórtice e o ponto decontrole. O seu valor, de acordo com a Eq. (3.1), é−Γ/(πc).

A condição de tangência à superfície do perfil impõe que a velocidade normal seja nula nestaregião. Assim, a soma das duas componentes de velocidade mencionadas acima deve ser nula. Estacondição é descrita pela Eq. (3.2).

V∞α− Γ

πc= 0 (3.2)

A intensidade do vórtice é obtida diretamente da Eq. (3.2). Isolando o termoΓ obtém-se:

Γ = παV∞c (3.3)

A sustentação gerada pelo perfil pode ser obtida por meio da aplicação do Teorema de Kutta-Joukowsky que estabelece a relação entre a força de sustentaçãoL e a circulaçãoΓ:

L = ρV∞Γ (3.4)

Substituindo o valor da circulação obtido em (3.3) na Eq. (3.4):

L = ρπαV 2∞c (3.5)

O coeficiente de sustentação da placa plana de envergadura unitária, para uma área de referênciaequivalente aSref = c.(1) = c, corresponde ao descrito pela Eq. (3.6).

CL =L

12ρV 2

∞Sref

=ρπαV 2

∞c12ρV 2

∞c= 2πα (3.6)

A derivada do coeficiente de sustentação com relação ao ângulo de ataque,CLα , para uma placaplana fica, então:

CLα = 2π (3.7)

A Equação (3.7) representa um resultado clássico para a Aerodinâmica. Apesar da simplicidadedo modelo utilizado, este resultado é exato para uma placa plana sujeita a um escoamento potencial.Isto se explica, no entanto, pela escolha conveniente do ponto de aplicação do vórtice e do ponto decontrole do perfil.

O mesmo valor deCLα pode ser obtido utilizando-se uma discretização mais refinada do perfil.Pode-se utilizar, por exemplo, uma divisão da corda do perfil emn painéis de comprimentos iguaisa ci = c/n. Em cada painel designado pelo índicei coloca-se, nos mesmos moldes do que foifeito anteriormente, um vórtice de intensidadeΓi no ponto a um quarto do comprimento do painel eestabelece-se um ponto de controle para cada painel situado a três quartos do seu comprimento. Asintensidades dos vórtices são determinadas aplicando-se a condição de tangência à superfície do perfilpara cada ponto de controle. Como cada painel possui uma intensidade de vórtice desconhecida e umponto de controle, o problema está determinado pois o número de equações corresponde ao númerode incógnitas.

3.1. MÉTODOS DE SINGULARIDADES 19

Seja, por exemplo, a velocidade normal induzida pelo vórtice posicionado a um quarto do com-primento do painelj no ponto de controle do paineli designada porvij. A velocidade normal totalinduzida por todos osn vórtices no paineli é, assim, dada pela Eq. (3.8).

vi =n∑

j=1

vij = −n∑

j=1

Γj

2π(xci − xv

j )(3.8)

ondexci corresponde à posição do ponto de controle do paineli e xv

j é a coordenada do vórtice dopainelj.

A velocidade normal total induzida no pontoi deve compensar a componente normal ao perfildo escoamento não-perturbado para que a condição de tangência do escoamento à superfície sejasatisfeita. Assim, para pequenos ângulos de ataque:

V∞α−n∑

j=1

Γj

2π(xci − xv

j )= 0 (3.9)

A aplicação da Eq. (3.9) para todos os pontos de controle leva a um sistema de equações linearesrepresentado na Eq. (3.10).

a11 a12 · · · a1(n−1) a1n

a21 a22 · · · a2(n−1) a2n...

......

......

an1 an2 · · · an(n−1) ann

Γ1

Γ2...

Γn

=

V∞αV∞α

...V∞α

(3.10)

Os termosaij da Eq. (3.10) dependem somente da geometria da discretização utilizada, uma vezque:

aij =1

2π(xci − xv

j )(3.11)

O sistema de equações representado pela Eq. (3.11) pode ser representado simbolicamente por:

[A]Γ = b (3.12)

O valor do vetorΓ pode ser obtido por intermédio da inversão da matrix[A], conforme mostradona Eq. (3.13), uma vez que o lado direito da Eq. (3.12),b, é conhecido.

Γ = [A]−1b (3.13)

O Apêndice A contém a listagem de um código computacional que traduz a teoria apresentadaacima. O referido programa foi executado para um ângulo de ataqueα = 1o e para uma discretizaçãohomogênea da corda equivalente an = 40. O resultado obtido para a distribuição deΓ ao longo dacorda é apresentada na Fig. 3.2. Nesta figura, o valorΓ0 é obtido da sustentação total do perfil e dadopela Eq. (3.14).

Γ0 =1

2V∞SrefcL (3.14)

Os resultados apresentados na Tab. 3.1 indicam que os valores deCLα fornecidos pelo códigocomputacional estão muito próximos do valor teórico deCLα = 2π. Além disso, a posição docentro de pressão também coincide com o valor teórico de1/4 da corda. Outro aspecto realçadoé a necessidade de se fazer um estudo paramétrico com a discretização para confirmar que a mesmanão está influenciando o resultado obtido.


Figura 3.2: Distribuição de sustentação ao longo da corda da placa plana.

Tabela 3.1: Resultados obtidos com o programa do Apêndice A.

n 10 20 30 40

CL 0,109662294 0,109662235 0,109662235 0,109662212Cm -0,0274155717 -0,0274155512 -0,0274155419 -0,0274155606CLα 6,28318691 6,2831831 6,2831831 6,28318214Cmα -1,57079661 -1,57079542 -1,57079482 -1,57079589

xCP /c 0,249999985 0,24999994 0,249999851 0,25000006

3.1.2 Escoamento potencial sobre uma asa finita

Para obter-se a solução do escoamento potencial sobre uma asa finita utiliza-se uma solução daEquação de Laplace conhecida como “Vórtice Ferradura”. Esta solução elementar é representadageometricamente na Fig. 3.3.

A idéia é dividir a superfície da asa em diversos quadriláteros, chamados comumente de painéis. Opainel é representado na Fig. 3.3 pelos pontosA, B, C eD. No interior de cada painel é posicionadoum vórtice em ferradura de intensidadeΓj e um ponto de controleP ′

c(x′, y′), conforme ilustrado.

A aplicação de um vórtice em ferradura necessita do conhecimento da velocidade normal induzidapor este tipo de vórtice em um ponto genérico da superfície da asaP ′(x′, y′). Esta velocidade érepresentada matematicamente pela Eq. (3.15). O sistema de coordenadasOx′y′ consiste de umsistema local, ou seja, posicionado no painel. A velocidade fornecida pela Eq. (3.15) é calculada emfunção das coordenadas locais do pontoP ′(x′, y′).

vz =−Γj

2πl′

1

ξ

[(1 + η)

RP+

(1− η)

RM

]+

1

(1− η)

[1 +

ξ

RM

]+

1

(1 + η)

[1 +

ξ

RP

](3.15)

3.1. MÉTODOS DE SINGULARIDADES 21

?x′,ξ

-y′,η

gP ′c(x

′, y′)

6

?

c′/4

6

?

c′/2

6

?

c′

- l′

- l′/2

?Γj

?Γj

-ΓjO

A B

CD

∞ ∞? ?

Figura 3.3: Representação de um vórtice ferradura.

ondeξ = 2x′/l′ eη = 2y′/l′. Os valores deRP eRM são fornecidos pela Eq. (3.16).

RP =√

ξ2 + (1 + η)2 RM =√

ξ2 + (1− η)2 (3.16)

O princípio utilizado para a solução do escoamento potencial sobre uma asa finita é o mesmoapresentado para a solução da placa plana. A velocidade normal induzida por todos os vórtices fer-radura solidários aos painéis que discretizam a superfície da asa deve anular a componente normaldo escoamento não-perturbado. A superfície da asa é dividida em painéis de forma a facilitar a im-plementação numérica do método. O mais comum é utilizar-se um arranjo conhecido como “VortexLattice”.

A velocidade normal induzida pelo vórtice ferradura do painelj no ponto de controle do paineli é representada, novamente, porvij. Note que a velocidade normal induzida pelo vórtice ferraduraé diretamente proporcional à sua intensidadeΓj e o fator de proporcionalidade depende apenas dageometria da discretização e pode ser representado poraij. O valor deaij é dado pela Eq. (3.15).Assim, a velocidade normal induzida por todos os vórtices no ponto de controle do paineli é dadapor:

vi =n∑

j=1

vij =n∑

j=1

aijΓj (3.17)

onden representa o número total de painéis da discretização da superfície da asa. Caso o número depainéis ao longo da corda e da envergadura sejam representados, respectivamente, pornc ens, entãon = ncns.

Para cada ponto de controle é aplicada a condição de tangência do escoamento à superfície da asa,


Figura 3.4: Distribuição de circulação ao longo da envergadura.

que equivale à Eq. (3.18).

V∞α +n∑

j=1

aijΓj = 0 (3.18)

Uma vez mais, como no caso da placa plana, a aplicação da condição de contorno em todos ospontos de controle leva a um sistema linear de equações que pode ser, novamente, representado pelaEq. (3.12). Neste caso, o vetorb é dado por−V∞α,−V∞α, ...,−V∞αT .

O método “Vortex-Lattice” descrito acima foi aplicado ao caso de uma asa retangular, de alonga-mentoA = 100, para comparação com os resultados obtidos com o código desenvolvido para aplaca plana. Como o alongamento da asa retangular é grande, espera-se que, longe das pontas, ascaracterísticas do escoamento sejam próximas daquelas para o escoamento bidimensional.

Foi utilizada uma discretização com 20 faixas ao longo da envergadura (ns = 20) e 10 painéisao longo da corda por faixa (nc = 10). A distribuição deΓ ao longo da envergadura é apresentadana Fig. 3.4. Nota-se que os valores deΓ não variam significativamente próximo do centro da asa,como era de se esperar. O efeito tridimensional causado pela presença de uma ponta na asa tambémé evidente.

A distribuição de circulação ao longo da corda é mostrada na Fig. 3.5. Nesta figura são apre-sentados os valores deΓ para cada faixa ao longo da envergadura calculados para a asa retangularde alongamentoA = 100. Além disso, os valores calculados pelo programa bidimensional parauma discretização equivalente, com 10 painéis (n = 10), são representados com símbolos diferentes.Observa-se que os resultados obtidos para o caso bidimensional concordam muito bem com os valorescalculados para as estações centrais da asa retangular.

3.2. DIFERENÇAS FINITAS 23

Figura 3.5: Distribuição de circulação ao longo da corda.

3.2 Diferenças finitas

Supondo que deseja-se resolver a Eq. (3.19 com condições iniciais e de contorno definidas.

∂u

∂t= α

∂2u

∂x2(3.19)

em queα é uma constante eu(x, t) é a função que se quer determinar.As condições iniciais e de contorno podem ser, por exemplo:

u(x, 0) = A sin(πx

L

)u(0, t) = 0 u(L, t) = 0

nas quaisA eL são constantes.Supondo, ainda, que o domínio de cálculo seja dividido em vários pontos separados pela mesma

distância∆x. Assim, tem-se o que é normalmente chamada de discretização do domínio de cálculo.O mesmo vale para o tempo. Pode-se considerar uma variávelt iniciando-se emt = 0 e crescendoem intervalos constantes∆t.

Desta forma, é possível avaliar a primeira derivada da funçãou num pontoxi utilizando a expan-são em Série de Taylor pelo lado direito deste ponto:

u(xi + ∆x) = u(xi) +∂u

∂x(xi)∆x +

1

2

∂2u

∂x2(xi)(∆x)2 +

1

6

∂3u

∂x3(xi)(∆x)3 + . . . (3.20)

na qual a variávelt foi omitida das expressões envolvendo a funçãou apenas para simplificar anotação.

Pode-se simplificar ainda mais a notação utilizandou(xi) = ui e u(xi + ∆x) = ui+1. Assim, aEq. (3.20) pode ser reescrita e rearranjada como:

ui+1 − ui =

(∂u

∂x

)i

∆x +1

2

(∂2u

∂x2

)i

(∆x)2 +1

6

(∂3u

∂x3

)i

(∆x)3 + . . . (3.21)


Dividindo toda a Eq. (3.21) por∆x:

ui+1 − ui

∆x=

(∂u

∂x

)i

+1

2

(∂2u

∂x2

)i

(∆x) +1

6

(∂3u

∂x3

)i

(∆x)2 + . . . (3.22)

Assim, a primeira derivada parcial deu no pontoxi pode ser estimada como:(∂u

∂x

)i

=ui+1 − ui

∆x− 1

2

(∂2u

∂x2

)i

(∆x)− 1

6

(∂3u

∂x3

)i

(∆x)2 − . . . (3.23)

Uma aproximação para a primeira derivada parcial deu no pontoxi pode, então, ser obtida comosendo: (

∂u

∂x

)i

≈ ui+1 − ui

∆x(3.24)

Nota-se, comparando as Eqs. (3.23) e (3.24), que a parte da Eq. (3.23) que foi desprezada começacom um termo contendo∆x. Por isso, a precisão da aproximação feita por meio da Eq. (3.24) é ditada ordem de∆x ou, ainda,O(∆x).

De forma semelhante, pode-se utilizar a expansão em Série de Taylor pelo lado esquerdo dexi:

u(xi −∆x) = u(xi)−∂u

∂x(xi)∆x +

1

2

∂2u

∂x2(xi)(∆x)2 − 1

6

∂3u

∂x3(xi)(∆x)3 + . . . (3.25)

Observa-se que os termos de ordem par apresentam o sinal positivo e os de ordem ímpar possuemsinal negativo. Utilizando a mesma simplificação na notação utilizada anteriormente, em queu(xi) =ui eu(xi −∆x) = ui−1:

ui−1 − ui = −(

∂u

∂x

)i

∆x +1

2

(∂2u

∂x2

)i

(∆x)2 − 1

6

(∂3u

∂x3

)i

(∆x)3 + . . . (3.26)

Dividindo toda a Eq. (3.26) por∆x:

ui−1 − ui

∆x= −

(∂u

∂x

)i

+1

2

(∂2u

∂x2

)i

(∆x)− 1

6

(∂3u

∂x3

)i

(∆x)2 + . . . (3.27)

Assim, a primeira derivada parcial deu no pontoxi também pode ser estimada como:(∂u

∂x

)i

=ui − ui−1

∆x+

1

2

(∂2u

∂x2

)i

(∆x)− 1

6

(∂3u

∂x3

)i

(∆x)2 + . . . (3.28)

Uma nova aproximação para a primeira derivada parcial deu no pontoxi pode, então, ser obtidada expressão: (

∂u

∂x

)i

≈ ui − ui−1

∆x(3.29)

Semelhantemente ao caso anterior, percebe-se que os termos negligenciados da Eq. (3.28) come-çam com um termo contendo∆x. Assim, esta aproximação para a primeira derivada parcial deu comrelação ax na pontoxi é da ordem de∆x ouO(∆x).

Uma expressão para avaliar a segunda derivada parcial deu com relação ax no pontoxi pode serobtida somando-se as Eqs. (3.22) e (3.27):

ui+1 − 2ui + ui−1 =

(∂2u

∂x2

)i

(∆x)2 +1

12

(∂3u

∂x3

)i

(∆x)4 + . . . (3.30)

3.2. DIFERENÇAS FINITAS 25

Dividindo a Eq. (3.30) por(∆x)2:

ui+1 − 2ui + ui−1

(∆x)2=

(∂2u

∂x2

)i

+1

12

(∂3u

∂x3

)i

(∆x)2 + . . . (3.31)

A segunda derivada fica:(∂2u

∂x2

)i

=ui+1 − 2ui + ui−1

(∆x)2− 1

12

(∂3u

∂x3

)i

(∆x)2 − . . . (3.32)

Consequentemente, a segunda derivada parcial deu com relação ax no pontoxi pode ser aproxi-mada como: (

∂2u

∂x2

)i

≈ ui+1 − 2ui + ui−1

(∆x)2(3.33)

Comparando-se as Eqs. (3.32) e (3.33), conclui-se que esta aproximação é da ordem de(∆x)2 ouO[(∆x)2].

Caso uma aproximação para a primeira derivada parcial deu com relação at no pontot(n) sejadesejada, pode-se utilizar o que já foi obtido parax com as devidas alterações para a variávelt. Aderivada, para uma expansão em Série de Taylor pela direita det(n), pode ser aproximada por:(

∂u

∂t

)(n)

≈ u(n+1) − u(n)

∆t(3.34)

em queu(n) = u(t(n)) eu(n+1) = u(t(n) + ∆t).Retornando à equação que se deseja resolver numericamente, a Eq. (3.19), e substituindo nela as

aproximações constantes das Eqs. (3.33) e (3.34):

u(n+1)i − u

(n)i

∆t+ O(∆t) = α

[u

(n)i+1 − 2u

(n)i + u

(n)i−1

(∆x)2

]+ O[(∆x)2] (3.35)

Simplificando a notação:

u(n+1)i − u

(n)i

∆t= α

[u

(n)i+1 − 2u

(n)i + u

(n)i−1

(∆x)2

]+ O[(∆t), (∆x)2] (3.36)

Assim, diz-se que este método numérico para resolução da Eq. (3.19) e primeira ordem de precisãono tempo [O(∆t)] e segunda ordem de precisão no espaço O[(∆x)2] . A Eq. (3.36) é chamada deEquação de Diferenças Finitas da Eq. (3.19).

A solução numérica baseia-se no fato de que a solução no instante inicial é conhecida em virtudeda condição inicial:

u(x, 0) = A sin(πx

L

)Assim, para o instante inicialt(0) para cadaxi da discretização no espaço, tem-se:

u(0)i = A sin

(πxi

L

)(3.37)

A solução no tempo seguinte pode ser obtida a partir da Eq. (3.36):

u(n+1)i ≈ u

(n)i +

α∆t

(∆x)2

[u

(n)i+1 − 2u

(n)i + u

(n)i−1

](3.38)


Outro detalhe importante é que a solução para os pontosx = 0 ex = L, para qualquer tempot, éconhecida pelas condições de contorno:

u(0, t) = 0 u(L, t) = 0

Desta forma, a solução pode ser obtida para uma discretização em quex0 = 0 exN = L aplicando,a cada passo no tempo, as condições de contorno:

u(n)0 = 0 u

(n)N = 0 (3.39)

em que, neste caso,N é o número de divisões deL e, consequentemente,∆x = L/N . O intervalo detempo∆t pode ser escolhido, assim como o número de divisõesN para a discretização do domíniode cálculo no espaço.

3.3 Volumes finitos

Admite-se, nesta seção, uma equação com duas dimensões no espaço e uma no tempo:

ρc∂T

∂t=

∂

∂x

(k∂T

∂x

)+

∂

∂y

(k∂T

∂y

)= ∇ · (k∇T ) (3.40)

Considerando que a Eq. (3.40) será integrada em um volume de controle igual aR:∫∫∫R

(ρc

∂T

∂t+∇ · q

)dR = 0 (3.41)

na qualq = −k∇T .Aplicando o Teorema do Divergente no segundo termo da Eq. (3.41):∫∫∫

R

ρc∂T

∂tdR +

∫∫S

q · n dS = 0 (3.42)

Considerando, também, quen = dx ı + dy eq = qx ı + qy , a Eq. (3.42) torna-se:∫∫∫R

ρc∂T

∂tdR +

∫∫S

(qx dy − qy dx) = 0 (3.43)

em que:

qx = −k∂T

∂xqy = −k

∂T

∂y

A Equação (3.43) é válida para qualquer volume de controleR. Na sua dedução não feita qualquerrestrição, por exemplo, quanto à forma do volume. Assim, para um paralelepípedo de lados∆x, ∆ye profundidade unitária, a primeira integral da Eq. (3.43) pode ser aproximado como:

ρcT

(n+1)i,j − T

(n)i,j

∆t∆x ∆y

A segunda integral da Eq. (3.43), baseando-se nas expressões deqx eqy e para uma área compostado perímetro do retângulo de lados∆x e∆y com profundidade unitária, é aproximada por:

−k ∆y

(∂T

∂x

)i+ 1

2,j

− k ∆x

(∂T

∂y

)i,j+ 1

2

+ k ∆y

(∂T

∂x

)i− 1

2,j

+ k ∆x

(∂T

∂y

)i,j− 1

2

3.4. ESTABILIDADE 27

Resta, neste momento, a avaliação dos termos correspondentes às derivadas parciais da variávelT com relação ax e ay. Admitindo uma malha com pontos igualmente espaçados em cada direçãoxey, estes termos podem ser aproximados por meio de diferenças centradas O[(∆x)2, (∆y)2]:(

∂T

∂x

)i+ 1

2,j

≈T

(n)i+1,j − T

(n)i,j

∆x

(∂T

∂x

)i− 1

2,j

≈T

(n)i,j − T

(n)i−1,j

∆x

sendo que, os valores das derivadas com relação ay são aproximados de forma similar, com a devidaalteração das variáveis.

Fazendo as devidas substituições:

ρcT

(n+1)i,j − T

(n)i,j

∆t∆x ∆y − k ∆y

T(n)i+1,j − T

(n)i,j

∆x− k ∆x

T(n)i,j+1 − T

(n)i,j

∆y

+k ∆yT

(n)i,j − T

(n)i−1,j

∆x+ k ∆x

T(n)i,j − T

(n)i,j−1

∆y= 0

(3.44)

Multiplicando toda a Eq. (3.44) por∆t/(∆x∆y) e agrupando os termos semelhantes:

ρc[T

(n+1)i,j − T

(n)i,j

]=

k∆t

∆x2

[T

(n)i+1,j − 2T

(n)i,j + T

(n)i−1,j

]+

k∆t

∆y2

[T

(n)i,j+1 − 2T

(n)i,j + T

(n)i,j−1

](3.45)

Nota-se que a aproximação por volumes finitos representada pela Eq. (3.45) corresponde a umaaproximação por diferenças finitas em que as derivadas parciais de primeira ordem são aproximadaspela Eq. (3.34) e as de segunda pela Eq. (3.33). Desta forma, a aproximação por volumes finitosconduziu a uma precisão deO[∆t, (∆x)2, (∆y)2].

O novo valor deTi,j para cada ponto da malha pode ser estimado, então, por:

T(n+1)i,j ≈ T

(n)i,j +

k∆t

ρc∆x2

[T

(n)i+1,j − 2T

(n)i,j + T

(n)i−1,j

]+

k∆t

ρc∆y2

[T

(n)i,j+1 − 2T

(n)i,j + T

(n)i,j−1

](3.46)

O problema pode ser resolvido pela utilização da Eq. (3.46) pois as condições iniciais e ascondições de contorno são conhecidas. Observa-se, também, a semelhança desta solução com a docaso unidimensional, representado pela Eq. (3.38), quando se substitui o parâmetroα pork/ρc.

3.4 Estabilidade


Exercícios

1. Faça a discretização de uma placa plana utilizando dois painéis correspondentes a cada metadeda corda da placa, colocando um vórtice a 1/4 da corda e um ponto de controle a 3/4 da corda decada painel. Compondo os escoamentos induzidos pelos dois vórtices com o escoamento uni-forme de velocidadeV∞ e sabendo que a placa plana encontra-se com um ângulo de ataqueα,descubra as equações do coeficiente de sustentação por unidade de envergadura, do coeficientede momento por unidade de envergadura em torno do bordo de ataque e a posição do centro depressão do perfil. Para a avaliação da posição do centro de pressão, deve-se utilizar a seguinterelação:

(xCP)BA

c= −(Cm′)BA

CL′

2. Obtenha uma aproximação para a primeira derivada de uma funçãou com relação ax, numcontexto de diferenças finitas, utilizando os pontos(xi, ui) e (xi+1, ui+1), tal quexi+1 − xi =∆x. Obtenha, também, a ordem do erro associado a esta aproximação.

3. Refaça o exercício anterior utilizando, desta vez, os pontos(xi−1, ui−1) e (xi, ui), de forma quexi − xi−1 = ∆x. Qual é a ordem do erro desta aproximação?

4. Obtenha uma aproximação para a primeira derivada de uma funçãou com relação ax, numcontexto de diferenças finitas, utilizando os pontos(xi−1, ui−1), (xi, ui) e (xi+1, ui+1), tal quexi − xi−1 = xi+1 − xi = ∆x. Obtenha, ainda, a ordem do erro associado a esta aproximação.

5. Forneça uma equação de diferenças finitas associada à equação diferencial parcial especificadaa seguir:

∂2u

∂t2= c2 ∂2u

∂x2

em quec é uma constante conhecida. Especifique a ordem da aproximação obtida.

6. Simule a solução numérica da equação diferencial parcial

∂u

∂t= α

∂2u

∂x2

em queα = 0, 2 e cuja equação de diferenças finitas associada pode ser

u(n+1)i − u

(n)i

∆t= α

[u

(n)i+1 − 2u

(n)i + u

(n)i−1

(∆x)2

]+ O[(∆t), (∆x)2]

Utilize a tabela a seguir para montar sua solução. A segunda linha corresponde aos valores dex e a segunda coluna aos det. Note que∆x = 0, 2 e que∆t = 0, 1. Utilize-a, também, como aespecificação das condições iniciais do problema. As condições de contorno são de temperaturaconstante nas extremidades do intervalo, ou seja, parax = 0 e x = 1. Pode ser utilizada umaplanilha tipo Excel para facilitar as contas. Neste caso, simule atét = 2, 0.

i 0 1 2 3 4 5n 0,0 0,2 0,4 0,6 0,8 1,0

0 0,0 10 20 30 30 20 101 0,12 0,23 0,34 0,45 0,5

Capítulo 4

Elementos de Aeroelasticidade

Aeroelasticidade é o estudo do efeito das forças aerodinâmicas em corpos elásticos.

Breve histórico:

Tacoma Narrows Bridge - Puget Sound, Washingtonem 7 de novembro de 1940, quatro mesesapós a ponte ter sido aberta ao tráfego, um vento de 42 mph induziu um movimento oscilatóriotorsional com um nó no centro do vão principal da ponte; havia, entretanto, um movimentoconsiderável no centro; a freqüência de oscilação mudou repentinamente de 37 para 14 ciclospor minuto; o movimento cresceu violentamente e o colapso da ponte ocorreu meia hora depois.

Fios de transmissão de energiaoscilações provocadas por ventos fortes; os cabos ou fios podemoscilar como um todo ou formando um ou mais nós ao longo da envergadura; quando as os-cilações tornam-se severas, os cabos movem-se irregularmente, mas de forma livre, atingindodistâncias verticais de 20 a 35 ft para envergaduras em torno de 500 ft.

Flutter de asas e superfícies de controle de aeronavesum dos primeiros casos documentados deflutter ocorreu com a empenagem horizontal do bombardeiro bimotor Handley Page O/400(biplano) no começo da I Guerra Mundial; problemas desta natureza com a asa começaram aaparecer quando os biplanos começaram a ser substituídos pelos monoplanos com uma rigideztorsional da asa menor (exemplo: Fokker D-8, da I Guerra Mundial).

A Figura 4.1 ilustra o movimento torsional e o colapso da Ponte Tacoma-Narrows.

4.1 Conceito

Aeroelasticidade é o ramo da Ciência que estuda a interação mútua entre forças aerodinâmicas eelásticas e a influência desta interação no projeto e na operação de aeronaves.

O diagrama apresentado na Fig. 4.2 mostra os diversos fenômenos decorrentes da interação dasforças aerodinâmicas (A), elásticas (E) e de inércia (I).

Alguns fenômenos são decorrentes da interação de apenas dois tipos de forças. Por exemplo, asvibrações mecânicas, como já visto anteriormente, são o produto da interação das forças elásticas ede inércia.

Existem fenômenos, no entanto, que são decorrentes da interação dos três tipos de força. Exem-plos são o flutter, o buffeting e a resposta dinâmica da aeronave.

Definições:

29

30 CAPÍTULO 4. ELEMENTOS DE AEROELASTICIDADE

Figura 4.1: Tacoma Narrows Bridge durante o flutter e o colapso.

Figura 4.2: Diagrama de interações das forças aerodinâmicas, elásticas e de inércia.

4.2. AEROELASTICIDADE ESTÁTICA 31

Flutter uma instabilidade dinâmica que ocorre em uma aeronave em vôo, a uma velocidade chamadade velocidade de flutter, em que a elasticidade da estrutura exerce um papel essencial na insta-bilidade.

Buffeting vibrações transientes de componentes estruturais da aeronave devido a impulsos aerodinâ-micos, produzidos por esteiras provenientes de asas, naceles, “pods´´ da fuselagem ou outroscomponentes da aeronave.

Resposta dinâmicaresposta transiente de componentes estruturais da aeronave produzida por car-gas aplicadas rapidamente tais como rajadas de vento, aterrissagem, movimentos abruptos docontrole, ondas de choque em movimento ou outras cargas dinâmicas.

Efeitos aeroelásticos na estabilidadeinfluência de deformações elásticas da estrutura na estabili-dade estática e dinâmica da aeronave.

Divergência uma instabilidade estática de uma superfície sustentadora de uma aeronave em vôo, auma velocidade chamada de velocidade de divergência, em que a elasticidade da superfíciesustentadora exerce um papel essencial na instabilidade.

Reversão de controleuma condição que ocorre em vôo, a uma velocidade chamada de velocidadede reversão de controle, na qual os efeitos esperados de deflexão de um dado componente dosistema de controle são completamente anulados por deformações elásticas da estrutura.

4.2 Aeroelasticidade estática

4.2.1 Divergência

Considerando a Figura 4.3, pode-se escrever, para a sustentação do aerofólio por unidade de en-vergadura:

Figura 4.3: Modelo esquemático de um aerofólio para análise de divergência.

L′ = qClc = qca(θ + α)

M ′0 = qCm0c

2

em que:

c → corda do perfilCl → coeficiente de sustentação bidimensional do perfila → derivada da curva de sustentação do perfil com o ângulo de ataqueq → pressão dinâmica =1/2ρV 2

α → ângulo de ataqueθ → ângulo de torsão


M ′a = M ′

0 + L′ec = qc2Cm0 + qec2a(θ + α)

M ′a = qec2aθ + qec2a

(α +

Cm0

ea

)

em que:ec → distância do Centro Aerodinâmico (CA) ao eixo de rotação (ou eixo elástico)e → razão entre a distância do CA ao eixo de rotação e a corda do aerofólio

Cm0 → coeficiente de momento para ângulo de ataque nulo

Fazendo−α0 = Cm0/ea eαe = α− α0:

M ′a = qec2a(θ + αe)

M ′e = Kαθ

Equilíbrio:M ′

a = M ′e ⇒ qec2a(θ + αe) = Kαθ

(qec2a−Kα)θ = −qec2aαe

Assim:

θ =qec2aαe

Kα − qec2a

A Divergência ocorre quando o denominador da equação anterior for nulo:

Kα − qdivec2a = 0 ⇒ qdiv =

Kα

ec2a

A Velocidade de Divergência pode ser obtida:

1

2ρV 2

div =Kα

ec2a

Finalmente:

Vdiv =

√2Kα

ρec2a(4.1)

4.2.2 Reversão de Controle

Considerando a Figura 4.4, pode-se escrever, para a sustentação do aerofólio por unidade de en-vergadura:

Figura 4.4: Modelo esquemático de um aerofólio para análise de reversão de controle.

4.2. AEROELASTICIDADE ESTÁTICA 33

CL =

(∂CL

∂ββ +

∂CL

∂αθ

)Cm =

∂Cm

∂ββ

T = qS

ec

(∂CL

∂ββ +

∂CL

∂αθ

)+ c

∂Cm

∂ββ

= Kαθ

Kα

qSθ = ec

(∂CL

∂ββ +

∂CL

∂αθ

)+ c

∂Cm

∂ββ

Kα

qS

θ

β= ec

(∂CL

∂β+

∂CL

∂α

θ

β

)+ c

∂Cm

∂β(Kα

qS− ec

∂CL

∂α

)θ

β= ec

∂CL

∂β+ c

∂Cm

∂β

θ

β=

ec∂CL

∂β+ c

∂Cm

∂βKα

qS− ec

∂CL

∂α

(4.2)

em que:

ec → distância do CA ao eixo de rotação (ou eixo elástico)e → razão entre a distância do CA ao eixo de rotação e a corda do aerofólioβ → deflexão da superfície de controleθ → torsão da asa em virtude da deflexão da superfície de controle

Reescrevendo a equação para oCL:

CL =

(∂CL

∂β+

∂CL

∂α

θ

β

)β (4.3)

Substituindo (4.2) em (4.3):

CL =

∂CL

∂β+

∂CL

∂α

ec∂CL

∂β+ c

∂Cm

∂βKα

qS− ec

∂CL

∂α

β

Após algum algebrismo, a equação anterior torna-se:

CL =

∂CL

∂β

Kα(∂CL

∂α

)qSec

+1

e

∂Cm

∂β

β

Kα(∂CL

∂α

)qSec

− 1

O aileron torna-se completamente inefetivo a partir deCL = 0. Neste momento, denomina-seqrev

a pressão dinâmica associada à Reversão de Controle. Assim:

∂CL

∂β

Kα(∂CL

∂α

)qrevSec

+1

e

∂Cm

∂β= 0


Novamente, após alguma álgebra:

qrev =

−Kα∂CL

∂β(∂CL

∂α

) (∂Cm

∂β

)Sc

=1

2ρV 2

rev

Isolando a velocidade de reversão de controleVrev:

Vrev =

√√√√√√√−Kα

∂CL

∂β(∂CL

∂α

) (∂Cm

∂β

)1

2ρSc

(4.4)

Como um resultado importante da derivação da Eq. (4.4), é importante observar que a velocidadede reversão de controle não depende da localização do eixo elástico (e).

Capítulo 5

Experimentos com estruturas

5.1 Extensometria elétrica

5.1.1 Conceitos básicos

Um fio, quando varia o seu comprimento, varia também a sua resistência à passagem de correnteelétrica. Esta característica foi aproveitada para a fabricação de dispositivos que medem a deformaçãoem estruturas. Sabe-se, também, que a deformação local de uma estrutura está associada à tensãonaquele local.

O dispositivo idealizado para a obtenção da deformação em um determinado local da estrutura éo Strain Gage, cuja forma mais comum encontra-se ilustrada na Fig. 5.1.

Figura 5.1: Modelo esquemático de umStrain Gage.

Na análise experimental de tensões, as deformações específicas podem ser determinadas por meioda variação da resistência elétrica dada pela Eq. (5.1).

∆R

R= kε (5.1)

na qual

k → constante característica do extensômetroR → resitência elétrica do fio [Ω]

∆R → variação da resitência elétrica do fio [Ω]ε → deformação específica do fio [m/m]

Exemplo:

35

36 CAPÍTULO 5. EXPERIMENTOS COM ESTRUTURAS

Para a investigação de deformações em peças metálicas, com deformações específicas deε =0, 001 m/m, empregando extensômetros elétricos de resistênciaR0 = 120Ω e constante característicak = 2, a variação da resistência será da ordem de:

∆R = 2× 0, 001× 120Ω

∆R = 0, 24Ω

A variação da resistência encontrada no exemplo anterior é típica de problemas nos quais são em-pregados extensômetros elétricos. Essa variação é muito baixa comparativamente ao valor da resistên-cia utilizada para a medida. Assim, para que seja possível medir deformações específicas menoresque 1 mm/m, torna-se necessária a utilização de um dispositivo que possa medir essa variação. Taldispositivo é chamado dePonte de Wheatstone.

Referências Bibliográficas

[1] FORTUNA, A. O. Técnicas Computacionais para Dinâmica dos Fluidos - Conceitos Básicos eAplicações, São Paulo: EdUSP, 2006.

[2] BISMARK-NASR, M. N. Structural Dynamics in Aeronautical Engineering, Reston, Virginia:AIAA Education Series, 1999.

[3] BISPLINGHOFF, R. L.; ASHLEY, H.; HALFMAN , R. L. Aeroelasticity, New York: Dover Pub-lications, 1996.

[4] CLARK , R. ET AL . A Modern Course in Aeroelasticity, 4th Edition, Kluwer Academic Publish-ers, 2004.

[5] DOWELL, E. H. ET AL . A Modern Course in Aeroelasticity, Rockville: Sijthoff & Noordhoff,1980.

[6] FUNG, Y. C. An Introduction to the Theory of Aeroelasticity, New York: Dover Publications,1993.

[7] SCHLICHTING, H.; TRUCKENBRODT, E. Aerodynamics of the airplane, New York: McGraw-Hill, 1979.

[8] ANDERSON, J. D. Fundamentals of Aerodynamics, 2nd Edition, New York:McGraw-Hill, 1991.

[9] BERTIN, J. J.; SMITH , M. L. Aerodynamics for Engineers, 3rd Edition, Upper SaddleRiver:Prentice-Hall, 1998.

[10] HALLIDAY , D.; RESNICK, R.; WALKER , J. Fundamentals of Physics, 4th Edition, New York:John Wiley & Sons, 1993.

[11] SCHMIDT, L. V. Introduction to Aircraft Flight Dynamics, AIAA Educational Series, 1998.

[12] ANDERSON, J. D. Introduction to Flight, McGraw-Hill, 199?.

[13] BOLTON, W. Engenharia de Controle, 1995.

37

38 REFERÊNCIAS BIBLIOGRÁFICAS

Apêndice A

Programa em FORTRAN para cálculo deuma Placa Plana

Este apêndice fornece um código computacional, escrito em FORTRAN, para o cálculo de umaplaca plana submetida a um escoamento potencial. O ângulo de ataque do perfil é de um grau. Acorda do aerofólio e a velocidade do escoamento encontram-se fixadas em uma unidade para cada umdestes parâmetros.

O número máximo de divisões da corda do perfil, para o cálculo das características aerodinâmicas,é definido pelo parâmetroNMAX que é fornecido no arquivo ’param.f’. Um exemplo desse arquivoé fornecido abaixo. Caso seja necessário aumentar o número máximo de divisões da corda, bastaalterar o valor estabelecido para o parâmetroNMAX.

Os programas utilizados neste texto são fornecidos para que o aluno tenha contato com ferramen-tas computacionais e para que possa familiarizar-se com este tipo de código. Além disto, o aluno éincentivado a modificar o código para que tenha acesso a informações que não foram disponibilizadasna versão fornecida aqui.

Alguns trabalhos do curso poderão ser elaborados a partir destes programas. Dessa forma, éimportante que o aluno tenha acesso ao código fonte e cuide para o seu constante aprimoramento.

integer nmaxparameter (nmax=50)

39

40 APÊNDICE A. PROGRAMA EM FORTRAN PARA CÁLCULO DE UMA PLACA PLANA

program perfilc--------------------------------------------------------------------c Propósito:c Resolver o escoamento potencial sobre uma placa plana utilizandoc vórtices bidimensionais. A corda do perfil é dividida em painéisc e cada um deles recebe um vórtice posicionado a 1/4 de seuc comprimento e um ponto de controle localizado a 3/4 desse mesmoc comprimento.cc Descrição das variáveis:c a - matriz de coeficientes de influênciac b - vetor de condições de contornoc gam - vetor contendo as intensidades dos vórticesc n - número de painéis utilizado na discretização da placa planac--------------------------------------------------------------------cc Declaração das variáveisc

implicit noneinclude ’param.f’integer n,i,jreal dx,pi,alpha,gamt,momt,cl,cm,cla,cma,xcpreal xc(nmax),xv(nmax),gam(nmax),b(nmax)real a(nmax,nmax)open(7,file=’cl_dist.dat’,status=’unknown’)

cc Inicialização das variáveisc

do i=1,nmaxdo j=1,nmax

a(i,j)=0.end doxc(i)=0.xv(i)=0.gam(i)=0.b(i)=0.

end docc Entrada da discretização do perfilc

print*, ’ ’print*, ’ Entre com o número de divisões da corda (N):’read(5,*) n

cc Teste de validade dos dados de entradac

if (n.gt.nmax) stop ’ N não deve ser maior que NMAX!’cc Construção dos pontos de posição dos vórtices e de controle

41

cpi=4.*atan(1.)dx=1./ndo i=1,n

xv(i)=(i-1)*dx+0.25*dxxc(i)=xv(i)+0.5*dx

end docc Construção do vetor de condições de contornoc

alpha=pi/180.do i=1,n

b(i)=alphaend do

cc Montagem da matriz dos coeficientes de influênciac

do i=1,ndo j=1,n

a(i,j)=1./(2.*pi*(xc(i)-xv(j)))end do

end docc Solução do sistema de equaçõesc

call invma(a,n)do i=1,n

do j=1,ngam(i)=gam(i)+a(i,j)*b(j)

end doend do

cc Cálculo da sustentação da placa planac

gamt=0.momt=0.do i=1,n

gamt=gamt+gam(i)write(7,*) xv(i),gam(i)momt=momt-xv(i)*gam(i)

end docl=2.*gamtcm=2.*momtcla=cl/alphacma=cm/alphaxcp=-cma/cla

cc Impressão dos resultadosc

42 APÊNDICE A. PROGRAMA EM FORTRAN PARA CÁLCULO DE UMA PLACA PLANA

write(6,*) ’ ’write(6,*) ’ Coeficiente de sustentação = ’,clwrite(6,*) ’ Coeficiente de momento = ’,cmwrite(6,*) ’ Derivada do coeficiente de sustentação = ’,clawrite(6,*) ’ Derivada do coeficiente de momento = ’,cmawrite(6,*) ’ Posição do centro de pressão = ’,xcpstopend

cc####################################################################c

subroutine invma(a,n)implicit noneinclude ’param.f’integer n,i,m,jreal a(nmax,nmax),p,q,p1

cdo 10 i=1,n

a(i,i)=a(i,i)+1.10 continue

do 20 m=1,np=a(m,m)-1.if (p.eq.0.) go to 500do 30 j=1,n

a(m,j)=a(m,j)/p30 continue

do 40 i=1,nif (i.eq.m) go to 40q=a(i,m)do 50 j=1,n

a(i,j)=a(i,j)-q*a(m,j)50 continue40 continue20 continue

do 60 i=1,na(i,i)=a(i,i)-1.

60 continue500 p1=abs(p)

if (p1.lt.10e-6) write(4,111) p111 format(1x,’PIVOT=’,f12.5)

returnend

Apêndice B

Programa em FORTRAN para o Método“Vortex-Lattice”

Este apêndice fornece um código computacional, escrito em FORTRAN, traduzindo uma imple-mentação do Método “Vortex-Lattice”. O ângulo de ataque é admitido como sendo fixo e igual aum grau. A geometria da planta da asa é fornecida por meio de um arquivo de entrada chamado’asa.dat’, que contém parâmetros tais como a corda na raiz, corda na ponta, enflechamento do bordode ataque e envergadura, conforme ilustrado na Fig. B.1. O número de Mach do escoamento e osnúmeros de faixas e de painéis por faixa também são fornecidos através deste arquivo. A velocidadedo escoamento encontra-se fixada em uma unidade.

6y

-x

@@

@@

@@

@@

AAAAAAAA

Λ

-ct

- cr

6

?

b

Figura B.1: Parâmetros geométricos da asa.

Este código computacional também possui um arquivo de parâmetros chamado ’param.f’ que ésintaticamente idêntico ao do programa apresentado no Ap. A. O parâmetroNMAX representa onúmero máximo de painéis. No entanto, para este caso, deve-se utilizar um número bem maior paraeste parâmetro, uma vez que a envergadura também está sendo discretizada. Um valor típico seriaNMAX = 500. Na verdade, o mesmo arquivo pode ser utilizado para ambos os códigos.

43

44 APÊNDICE B. PROGRAMA EM FORTRAN PARA O MÉTODO “VORTEX-LATTICE”

program vortexc--------------------------------------------------------------------c Propósito: Calcular derivada do coeficiente de sustentaçãoc com o ângulo de ataque e a posição do centro de pressão dec uma asa trapezoidal, utilizando o método Vortex-Lattice ec admitindo que o vórtice ferradura possui ângulos retos.c A área de referência é a área em planta da asa e oc comprimento de referência corresponde à corda na raiz. Ac correção de compressibilidade é feita através do método dec Prandtl-Glauert.cc Variáveis de entrada:c cr => corda na raiz da asac ct => corda na ponta da asac span => envergadura da asac sle => enflechamento do bordo de ataque em grausc nc => número de painéis ao longo da cordac ns => número de painéis ao longo da envergadurac mach => número de Mach do escoamento não-perturbadocc Variáveis de saída:c cla => derivada do coeficiente de sustentação com oc ângulo de ataquec cma => derivada do coeficiente de momento de arfagemc com o ângulo de ataquec xcp => posição do centro de pressão, a partir do bordoc de ataque na raiz da asa, em fração da corda nac raizc--------------------------------------------------------------------

implicit noneinclude ’param.f’integer i,j,nc,ns,n,ii,jj,ic,jcreal cr,ct,span,sle,ss,dx,vz,alpha,sref,machreal pi,cl,cla,cm,cma,xcp,gamt,momt,gaml,gams,xle,yc,ccreal y1,y2,dtheta,thet1,thet2,betareal xc(nmax,nmax),xv(nmax,nmax),y(nmax),dy(nmax)real a(nmax,nmax),b(nmax),gam(nmax)open(3,file=’asa.dat’,status=’unknown’)open(7,file=’gamma.dat’,status=’unknown’)open(8,file=’geo.dat’,status=’unknown’)open(9,file=’gam_s.dat’,status=’unknown’)

cc Inicializando as variáveisc

do i=1,nmaxdo j=1,nmax

a(i,j)=0.end dob(i)=0.

45

gam(i)=0.end do

cc Entrada de dados da geometria da asa e da discretizaçãoc

read(3,*) cr,ct,span,sle,nc,ns,machcc Verificação da integridade dos dados de entradac

n=nc*nsif (n.gt.nmax) stop ’ NC*NS não deve ser maior que NMAX!’if (mach.ge.1.) stop ’ O número de MACH deve ser menor que 1.0!’

cc Transformação da geometria em função do número de Machc (Prandtl-Glauert)c

pi=3.1415927sle=pi*sle/180.beta=sqrt(1.-mach*mach)cr=cr/betact=ct/betasle=atan(tan(sle)/beta)

cc Cálculo dos pontos dos vórtices e de controlec

alpha=pi/180.sref=(cr+ct)*span/2.ss=span/2.dtheta=pi/ns

cc Divisão da envergadura em ns partes, concentrando pontosc nas pontasc

do j=1,nsthet1=j*dthetathet2=(j-1)*dthetay1=ss*cos(pi-thet1)y2=ss*cos(pi-thet2)y(j)=0.5*(y1+y2)dy(j)=y1-y2yc=abs(y(j))xle=yc*tan(sle)cc=(cr-ct)/ss*(ss-yc)+ct

cc Divisão da corda em nc painéis de mesma cordac

dx=cc/ncdo i=1,nc

xv(i,j)=xle+(i-1)*dx+0.25*dx


xc(i,j)=xv(i,j)+0.5*dxwrite(8,*) xv(i,j),xc(i,j),y(j)

end dowrite(8,*) ’ ’

end docc Cálculo dos coeficientes da matriz de influênciac

do jj=1,nsdo ii=1,nc

ic=(jj-1)*nc+iido j=1,ns

do i=1,ncjc=(j-1)*nc+icall vind(xv(i,j),y(j),xc(ii,jj),y(jj),dy(j),vz)a(ic,jc)=vz

end doend dob(ic)=alpha

end doend do

cc Resolução do sistema linear de equaçõesc

call invma(a,n)do i=1,n

do j=1,ngam(i)=gam(i)+a(i,j)*b(j)

end doend do

cc Cálculo da sustentação da asac

gamt=0.momt=0.do j=1,ns

gams=0.do i=1,nc

ic=(j-1)*nc+igaml=gam(ic)*dy(j)gamt=gamt+gamlgams=gams+gam(ic)write(7,*) xv(i,j),gam(ic)momt=momt-xv(i,j)*gaml

end dowrite(7,*) ’ ’write(9,*) y(j),gams

end doc

47

c Cálculo dos coeficientes aerodinâmicos com a correção dec Prandtl-Glauertc

cl=2.*gamt/sref/betacm=2.*momt/sref/cr/betacla=cl/alphacma=cm/alphaxcp=-cma/cla

cc Impressão dos resultadosc

write(6,*) ’ ’write(6,*) ’ Número de Mach = ’,machwrite(6,*) ’ Derivada do coeficiente de sustentação = ’,clawrite(6,*) ’ Derivada do coeficiente de momento = ’,cmawrite(6,*) ’ Posição do centro de pressão/corda na raiz = ’,xcpstopend

cc################################################################c

subroutine vind(xv,yv,x,y,xl,vn)implicit nonereal xv,yv,x,y,xl,vnreal pi,xi,eta,rp,rm,vi1,vi2,vi3xi=2.*(x-xv)/xleta=2.*(y-yv)/xlrp=sqrt(xi*xi+(1.+eta)*(1.+eta))rm=sqrt(xi*xi+(1.-eta)*(1.-eta))vi1=((1.+eta)/rp+(1.-eta)/rm)/xivi2=(1.+xi/rm)/(1.-eta)vi3=(1.+xi/rp)/(1.+eta)vn=(vi1+vi2+vi3)/(6.2831854*xl)returnend

cc################################################################c

subroutine invma(a,n)implicit noneinclude ’param.f’integer n,i,m,jreal a(nmax,nmax),p,q,p1

cdo 10 i=1,n

a(i,i)=a(i,i)+1.10 continue

do 20 m=1,np=a(m,m)-1.


if (p.eq.0.) go to 500do 30 j=1,n

a(m,j)=a(m,j)/p30 continue

do 40 i=1,nif (i.eq.m) go to 40q=a(i,m)do 50 j=1,n

a(i,j)=a(i,j)-q*a(m,j)50 continue40 continue20 continue

do 60 i=1,na(i,i)=a(i,i)-1.

60 continue500 p1=abs(p)

if (p1.lt.10e-6) write(4,111) p111 format(1x,’PIVOT=’,f12.5)

returnend

Documents

Tratamento Numérico de Escoamentos Aerodinâmicos'breno.freeshell.org/Aer_num.pdf · 2009. 9. 20. · Breno Moura Castro Tratamento Numérico de Escoamentos Aerodinâmicos SÃO