Download PDF
ads:
UNIVERSIDADE FEDERAL DE ITAJUBÁ
INSTITUTO DE ENGENHARIA MECÂNICA
PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA MECÂNICA
DISSERTAÇÃO DE MESTRADO
Análise do Comportamento Aerodinâmico de
um Corpo na Presença de uma Superfície
Plana Móvel
Autor: Alex Mendonça Bimbato
Orientador: Prof. Dr. Luiz Antonio Alcântara Pereira
Co-Orientador: Prof. Ph.D. Miguel Hiroo Hirata
Itajubá, 11 de Julho de 2008
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
UNIVERSIDADE FEDERAL DE ITAJUBÁ
INSTITUTO DE ENGENHARIA MECÂNICA
PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA MECÂNICA
DISSERTAÇÃO DE MESTRADO
Análise do Comportamento Aerodinâmico de
um Corpo na Presença de uma Superfície
Plana Móvel
Autor: Alex Mendonça Bimbato
Orientador: Prof. Dr. Luiz Antonio Alcântara Pereira
Co-Orientador: Prof. Ph.D. Miguel Hiroo Hirata
Curso: Mestrado em Engenharia Mecânica
Área de Concentração: Dinâmica dos Fluidos e Máquinas de Fluxo
Dissertação submetida ao Programa de Pós-Graduação em Engenharia Mecânica como
parte dos requisitos para obtenção do Título de Mestre em Engenharia Mecânica.
Itajubá, 11 de Julho de 2008
M.G. - Brasil
ads:
UNIVERSIDADE FEDERAL DE ITAJUBÁ
INSTITUTO DE ENGENHARIA MECÂNICA
PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA MECÂNICA
DISSERTAÇÃO DE MESTRADO
Análise do Comportamento Aerodinâmico de
um Corpo na Presença de uma Superfície
Plana Móvel
Autor: Alex Mendonça Bimbato
Orientador: Prof. Dr. Luiz Antonio Alcântara Pereira
Co-Orientador: Prof. Ph.D. Miguel Hiroo Hirata
Composição da Banca Examinadora:
Prof. Ph.D. Gustavo César Rachid Bodstein - COPPE/UFRJ
Profª. Drª. Ana Lúcia Fernandes de Lima e Silva - UNIFEI
Prof. Dr. Waldir de Oliveira - UNIFEI
Prof. Ph.D. Miguel Hiroo Hirata (Co-Orientador) - FAT/UERJ
Prof. Dr. Luiz Antonio Alcântara Pereira (Orientador) - UNIFEI
Dedicatória
Dedico este trabalho ao meu Deus e à minha querida família.
Agradecimentos
Agradeço aos meus pais, Jair e Nancy, pelo total apoio recebido durante todo o tempo
de elaboração deste trabalho.
Ao Professor Luiz Antonio Alcântara Pereira, por ter me apresentado o Método de
Vórtices, pela orientação clara, segura e objetiva, e pela cobrança para que este trabalho
pudesse ser concluído com êxito.
Ao Professor Miguel Hiroo Hirata, pela sua indispensável ajuda na co-orientação deste
trabalho.
À Takafumi Nishino, da Universidade de Southampton na Inglaterra, pelos valiosos
esclarecimentos prestados em relação aos resultados dos experimentos apresentados na sua
Tese de Doutorado.
Meu muito obrigado aos amigos Carlos Adriano, Jan e Washington, pelo auxílio no uso
de ferramentas da informática.
Ao Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), pelo
apoio financeiro através da concessão de uma bolsa de mestrado.
“No campo daqueles que procuram a verdade, não existe nenhuma autoridade humana. Todo
aquele que se fizer de magistrado encontrará imediatamente a risada dos deuses”
Albert Einstein.
Resumo
BIMBATO, A. M. (2008), Análise do Comportamento Aerodinâmico de um Corpo na
Presença de uma Superfície Plana Móvel, Itajubá, 131p. Dissertação (Mestrado em Dinâmica
dos Fluidos e Máquinas de Fluxo) – Instituto de Engenharia Mecânica, Universidade Federal
de Itajubá.
As características do escoamento ao redor de um corpo que se encontra nas
proximidades de uma superfície plana são governadas não só pelo número de Reynolds, Re,
mas também pela relação que define a distância do corpo até a superfície plana, dh . Os
detalhes dos fenômenos associados à presença do efeito solo ainda estão longe de serem
completamente entendidos, devido a uma variedade de outros fatores influentes, em
particular, a interferência entre a camada limite que se forma sobre o solo e a esteira formada
a jusante do corpo. O propósito de fazer com que a superfície plana se mova em relação ao
corpo está no fato de se eliminar a vorticidade gerada no solo, que consiste no principal
mecanismo que dificulta a compreensão do efeito solo. Neste trabalho apresenta-se um
algoritmo do Método de Vórtices para a análise do escoamento bidimensional e em regime
não-permanente de um fluido viscoso ao redor de um cilindro de seção circular na presença de
uma superfície plana móvel, onde a geração de vórtices discretos se dá apenas sobre a
superfície do corpo. As cargas aerodinâmicas atuantes sobre o corpo são calculadas a partir de
uma formulação integral derivada de uma equação de Poisson para a pressão e os resultados
numéricos obtidos são comparados com resultados experimentais disponíveis na literatura.
Palavras-chave
Método de Vórtices, Método de Painéis, Superfície Plana Móvel, Supressão de Camada
Limite, Cargas Aerodinâmicas, Descrição Lagrangiana.
Abstract
BIMBATO, A. M. (2008), Analysis of Moving Ground Effects on Aerodynamics Loads of a
Body, MSc. Dissertation – Instituto de Engenharia Mecânica, Universidade Federal de Itajubá,
131p.
The characteristics of the flow around a body placed near a plane boundary are
governed not only by the Reynolds number, Re, but also by the “gap ratio”, i.e., the ratio of
the gap distance between the body and the plane boundary, h, to the body diameter, d. The
details of the effects of dh , or “ground effect”, are still far from being fully understood due
the variety of other influencing factors, in particular the interference between the boundary
layer formed on the ground and the wake formed behind the body. The purpose of using the
moving ground is to eliminate the discussed interference, which is probably one of the most
crucial but that makes the problem understanding difficult. The present work deals with the
Vortex Method to simulate numerically the two dimensional, unsteady and viscous flow
around a circular cylinder in the vicinity of a moving ground. The vorticity which is generated
from the body surface is discretized and represented by a cloud of Lamb discrete vortices.
Each vortex in the cloud is followed in a purely Lagrangian description. In this methodology
the vorticity is not generated on the ground. The aerodynamics loads are computed using an
integral formulation derived from the Poisson equation for the pressure. The numerical results
for a circular cylinder in a moving ground are compared with experimental results available in
the literature.
Keywords
Vortex Method, Panels Method, Moving Ground, Boundary Layer Suppression,
Aerodynamics Loads, Lagrangian Description.
i
Sumário
LISTA DE FIGURAS
iv
LISTA DE TABELAS
viii
SIMBOLOGIA
ix
LETRAS LATINAS
ix
LETRAS GREGAS
xi
ABREVIATURAS
xii
SIGLAS
xiii
CAPÍTULO 1: INTRODUÇÃO
1
CAPÍTULO 2: REVISÃO BIBLIOGRÁFICA
7
2.1 – INTRODUÇÃO 7
2.2 – O MÉTODO DE VÓRTICES DISCRETOS 7
2.3 – O CILINDRO CIRCULAR ISOLADO 13
2.4 – O EFEITO SOLO 17
CAPÍTULO 3: FORMULAÇÃO GERAL DO PROBLEMA
23
3.1 – INTRODUÇÃO 23
3.2 – GEOMETRIA E DEFINIÇÕES 23
3.3 – HIPÓTESES SIMPLIFICADORAS 25
3.4 – EQUAÇÕES GOVERNANTES E CONDIÇÕES DE CONTORNO 26
3.4.1 – Equações Governantes 26
3.4.2 – Condições de Contorno 27
3.5 – ADIMENSIONALIZAÇÃO DO PROBLEMA 28
3.6 – EQUAÇÃO DO TRANSPORTE DA VORTICIDADE 31
CAPÍTULO 4: MÉTODO DE SOLUÇÃO: O MÉTODO DE VÓRTICES
33
4.1 – INTRODUÇÃO 33
4.2 – CARACTERÍSTICAS MARCANTES DO MÉTODO DE VÓRTICES 33
ii
4.3 – A CONVECÇÃO DA VORTICIDADE 35
4.3.1 – Contribuição do Escoamento Incidente 35
4.3.2 – Contribuição do Corpo: O Método dos Painéis 36
4.3.3 – Geração de Vorticidade 39
4.3.4 – Contribuição da Nuvem de Vórtices 44
4.3.5 – Esquemas de Avanço Convectivo 46
4.4 – A DIFUSÃO DA VORTICIDADE 47
4.5 – CARGAS AERODINÂMICAS 48
CAPÍTULO 5: IMPLEMENTAÇÃO DO ALGORITMO DO MÉTODO DE
VÓRTICES
52
5.1 – INTRODUÇÃO 52
5.2 – DESCRIÇÃO DO FUNCIONAMENTO DAS ROTINAS 52
5.2.1 – Rotina INPUT.FOR 53
5.2.2 – Rotina DATAPRGP.FOR 55
5.2.3 – Rotina DATAPRSB.FOR 55
5.2.4 – Rotina COUPS.FOR 56
5.2.5 – Rotina COUPV.FOR 57
5.2.6 – Rotina RHSS.FOR 57
5.2.7 – Rotina RHSV.FOR 58
5.2.8 – Rotina GAUSSPIV.FOR 58
5.2.9 – Rotina MODCOUP.FOR 58
5.2.10 – Rotina GENERAT.FOR 59
5.2.11 – Rotina COMPUMVM.FOR 59
5.2.12 – Rotina COMPUCVC.FOR 59
5.2.13 – Rotina PRESSURE.FOR 60
5.2.14 – Rotina DIFFUS.FOR 60
5.2.15 – Rotina CONVEC.FOR 60
5.2.16 – Rotina REFLECT.FOR 60
5.2.17 – Rotina REFLECS.FOR 60
5.2.18 – Rotina CONSERV.FOR 61
5.2.19 – Rotina MODM.FOR 61
5.2.20 – Rotina PRINT.FOR 61
5.2.21 – Rotina AVERAGE.FOR 61
CAPÍTULO 6: ANÁLISE DE RESULTADOS
63
6.1 – INTRODUÇÃO 63
iii
6.2 – PARÂMETROS UTILIZADOS NA SIMULAÇÃO NUMÉRICA 64
6.2.1 – Parâmetros Relacionados com o Método Numérico 64
6.2.2 – Parâmetros Relacionados com o Fenômeno Físico 66
6.3 – ESCOAMENTO AO REDOR DE UM CILINDRO CIRCULAR
AAAAIIISOLADO
68
6.4 – ESCOAMENTO AO REDOR DE UM CILINDRO CIRCULAR
AAAAIISUBMETIDO AO EFEITO SOLO
77
CAPÍTULO 7: CONCLUSÕES E SUGESTÕES
88
7.1 – CONCLUSÕES 88
7.2 - SUGESTÕES 91
REFERÊNCIAS BIBLIOGRÁFICAS
93
APÊNDICE A: DISTRIBUIÇÃO DA VORTICIDADE E DA
VELOCIDADE INDUZIDA
103
A.1 – O VÓRTICE POTENCIAL 103
A.2 – O VÓRTICE DE LAMB 105
APÊNDICE B: RESULTADOS ADICIONAIS DA SIMULAÇÃO
NUMÉRICA
108
iv
Lista de Figuras
Figura 2.1 Visualização do escoamento para uma ampla faixa de números
de Reynolds. (Figuras retiradas de Batchelor, 1967; Schlichting,
1979; van Dyke, 1982; Tritton, 1988).
15
Figura 2.2 Análise das cargas aerodinâmicas atuantes sobre um cilindro
circular na presença do efeito solo. (Figura retirada do trabalho
de Nishino (2007)).
21
Figura 3.1 Definição da geometria e da região fluida do problema com o
sistema de coordenadas fixo ao corpo.
24
Figura 3.2 Especificação das condições de contorno sobre as fronteiras que
delimitam o domínio fluido.
27
Figura 3.3 Representação do problema adimensionalizado. 31
Figura 4.1 Condições de Contorno. 36
Figura 4.2 Distribuição constante de fontes sobre o eixo x. 37
Figura 4.3 Velocidade induzida no ponto P(x,y) por uma distribuição de
fontes com densidade constante,
(
)
xσ , ao longo de uma
superfície de comprimento
(
)
12
xx
.
38
Figura 4.4 Geração de vorticidade: um enfoque fenomenológico. 40
Figura 4.5 Geração de vórtices discretos de Lamb sobre a superfície
discretizada de um corpo.
41
Figura 4.6 Comportamento da velocidade tangencial induzida. 45
Figura 5.1 Estrutura do programa computacional
“MOVINGROUND.FOR”.
53
Figura 5.2 Representação esquemática do corpo e da superfície plana
móvel.
54
Figura 5.3 Representação dos pontos extremos dos painéis, dos pontos de 56
v
controle dos painéis e dos pontos de desprendimento de vórtices
discretos.
Figura 6.1 Cilindro circular localizado próximo a uma superfície plana
móvel.
64
Figura 6.2 Distribuição média de pressão ao longo da superfície
discretizada do cilindro circular isolado (Euler, 300mb
=
,
0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re =
).
70
Figura 6.3 Evolução das cargas aerodinâmicas integradas ao longo do
tempo para o cilindro circular isolado (Euler,
300mb
=
,
0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re = ).
72
Figura 6.4 Distribuição instantânea de pressão ao longo da superfície
discretizada do cilindro circular isolado (Euler,
300mb
=
,
0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re = ).
72
Figura 6.5 Detalhes do desprendimento de estruturas vorticosas em um
cilindro circular isolado no instante 39,4t
=
(Euler,
300mb =
,
0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re = ).
73
Figura 6.6 Detalhes do desprendimento de estruturas vorticosas em um
cilindro circular isolado no instante 6,40t
=
(Euler, 300mb = ,
0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re =
).
74
Figura 6.7 Detalhes do desprendimento de estruturas vorticosas em um
cilindro circular isolado no instante
6,41t
=
(Euler, 300mb = ,
0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re = ).
75
Figura 6.8 Detalhes do desprendimento de estruturas vorticosas em um
cilindro circular isolado no instante 9,42t
=
(Euler,
300mb =
,
0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re = ).
75
Figura 6.9 Detalhes do desprendimento de estruturas vorticosas em um
cilindro circular isolado no instante 0,44t
=
(Euler, 300mb = ,
0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re =
).
76
Figura 6.10 Posição dos vórtices na esteira ao final da simulação numérica
do cilindro circular isolado (Euler,
300mb
=
, 0,05t
=
,
0,001epsσ
0
=
=
,
5
10Re = ).
77
Figura 6.11 Comparação entre os resultados numéricos e experimentais 78
vi
obtidos para o coeficiente de arrasto com relação aos efeitos
tridimensionais.
Figura 6.12 Aparato experimental usado nos testes em túnel de vento por
Nishino (2007).
79
Figura 6.13 Comparação entre resultados numéricos e experimentais obtidos
para o coeficiente de sustentação envolvendo o cilindro circular
submetido ao efeito solo.
80
Figura 6.14 Comparação da evolução das cargas aerodinâmicas para os casos
do cilindro circular (a) submetido ao efeito solo e (b) isolado
(Euler,
300mb = ,
0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ).
81
Figura 6.15 Distribuição instantânea de pressão ao longo da superfície
discretizada do cilindro circular submetido ao efeito solo (Euler,
300mb
=
, 0,05t = , 0,001epsσ
0
=
=
,
5
10Re = , 0,95dh = ).
82
Figura 6.16 Detalhes do desprendimento de estruturas vorticosas em um
cilindro circular submetido ao efeito solo no instante
4,46t =
(Euler,
300mb = , 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ,
0,95dh = ).
83
Figura 6.17 Detalhes do desprendimento de estruturas vorticosas em um
cilindro circular submetido ao efeito solo no instante
6,47t =
(Euler,
300mb = , 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ,
0,95dh = ).
84
Figura 6.18 Detalhes do desprendimento de estruturas vorticosas em um
cilindro circular submetido ao efeito solo no instante
7,48t =
(Euler,
300mb = , 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ,
0,95dh = ).
85
Figura 6.19 Detalhes do desprendimento de estruturas vorticosas em um
cilindro circular submetido ao efeito solo no instante
1,50t =
(Euler,
300mb = , 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ,
0,95dh = ).
85
Figura 6.20 Detalhes do desprendimento de estruturas vorticosas em um
cilindro circular submetido ao efeito solo no instante
51,2t =
(Euler,
300mb = , 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ,
86
vii
0,95dh = ).
Figura 6.21 Posição dos vórtices na esteira ao final da simulação numérica
do cilindro circular submetido ao efeito solo (Euler,
300mb = ,
0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re = , 0,95dh
=
).
87
Figura A.1 Velocidade induzida (retirada de Alcântara Pereira, 1999). 104
Figura A.2 Velocidade induzida pelo vórtice potencial (retirada de
Alcântara Pereira, 1999).
105
Figura A.3 Velocidade induzida pelo vórtice de Lamb (retirada de Alcântara
Pereira, 1999).
107
Figura B.1 Evolução das cargas aerodinâmicas ao longo do tempo para um
cilindro circular estacionado a pequenas distâncias da superfície
plana móvel (Euler,
300mb
=
, 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ).
109
Figura B.2 Evolução das cargas aerodinâmicas ao longo do tempo para um
cilindro circular estacionado a distâncias intermediárias da
superfície plana móvel (Euler,
300mb
=
, 0,05t
=
,
0,001epsσ
0
=
=
,
5
10Re = ).
110
Figura B.3 Posição dos vórtices na esteira ao final da simulação numérica
do cilindro circular submetido ao efeito solo para alguns valores
de
dh
(Euler, 300mb
=
, 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ).
110
viii
Lista de Tabelas
Tabela 6.1 Valores médios dos coeficientes de arrasto e de sustentação e
número de Strouhal para um cilindro circular isolado.
69
Tabela 6.2 Valor médio do coeficiente de arrasto e número de Strouhal para
o cilindro circular submetido ao efeito solo para
0,95dh
=
.
80
Tabela B.1 Valores médios dos coeficientes de arrasto e de sustentação e
número de Strouhal para um cilindro circular submetido ao
efeito solo para vários valores da relação
dh .
108
ix
Simbologia
Letras Latinas
Ad Matriz lado direito de pressão
D
C
Coeficiente de arrasto
'
D
C
Valor rms do coeficiente de arrasto
L
C
Coeficiente de sustentação
'
L
C
Valor rms do coeficiente de sustentação
co Ponto de controle de um painel plano
core
Raio do núcleo do vórtice de Lamb (
0
σcore
=
)
P
C
Coeficiente de pressão
m
P
C
Coeficiente de pressão em um painel genérico m
D Força de arrasto
d Diâmetro do cilindro circular
delt Incremento de tempo
dS
Elemento orientado de um contorno
eps Distância para a geração de vórtices discretos a partir do co
f Freqüência de desprendimento de pares contra-rotativos de vórtices
h Distância entre a superfície plana e o centro do cilindro circular
()
c
dh
Valor adimensional e crítico da distância do cilindro circular à
superfície plana móvel
L Força de sustentação
lm Comprimento de cada módulo que compõe a superfície plana
M Número total de painéis planos utilizados para discretizar a
x
superfície do corpo e a superfície plana móvel
m Painel arbitrário
Ma Número de Mach
mb Número de painéis utilizados na discretização do cilindro circular
mg Número de painéis utilizados na discretização da superfície plana
N Número de vórtices presentes na nuvem
nm Número de módulos que compõem a superfície plana
np Número de painéis utilizados para discretizar cada módulo da
superfície plana
P Número randômico entre 0 e 1
p Campo de pressões
j
p
Pressão no ponto de controle do painel j
Pr Número de Prandtl
pshed Ponto de desprendimento de vórtice
p
Pressão de referência
Q Número randômico entre 0 e 1
r Distância do ponto extremo de um painel até um ponto arbitrário do
domínio fluido
Re Número de Reynolds
kj
r
Distância entre o vórtice k e o vórtice j
m
k
r
Distância entre o vórtice k e o painel m
S Fronteira que delimita o domínio fluido semi-infinito
St Número de Strouhal
start Instante de início do cálculo das cargas aerodinâmicas
step Passo em que se encontra a simulação numérica
stop Número total (final) de incrementos de tempo da simulação
numérica
t Instante de tempo da simulação numérica
U Velocidade do escoamento incidente
u
Campo de velocidades
u Componente na direção x da velocidade total induzida pelos painéis
que depende da distribuição de fontes
(
)
xσ
N
k
u
Componente na direção x da velocidade total induzida no vórtice k
pela nuvem de vórtices
xi
mk
u
Componente na direção x da velocidade induzida por um vórtice
arbitrário k sobre o ponto de controle de um painel genérico m
k
t
u
Componente na direção x da velocidade total induzida no vórtice k
jk,
V
U
Velocidade induzida no vórtice arbitrário k pelo vórtice arbitrário j
na direção do eixo dos x, se este último possuir intensidade unitária
u
Componente na direção x da velocidade do escoamento incidente
v Componente na direção y da velocidade total induzida pelos painéis
que depende da distribuição de fontes
(
)
xσ
N
k
v
Componente na direção y da velocidade total induzida no vórtice k
pela nuvem de vórtices
mk
v
Componente na direção y da velocidade induzida por um vórtice
arbitrário k sobre o ponto de controle de um painel genérico m
k
t
v
Componente na direção y da velocidade total induzida no vórtice k
jk,
V
V
Velocidade induzida no vórtice arbitrário k pelo vórtice arbitrário j
na direção do eixo dos y, se este último possuir intensidade unitária
v
Componente na direção y da velocidade do escoamento incidente
x Coordenada cartesiana no plano complexo z
k
x
Componente na direção x da posição de um vórtice k
m
XM
Abscissa do ponto de controle do painel m
Y Trabalho específico total
y Coordenada cartesiana no plano complexo z
dy
e
Distância entre o cilindro circular e a extremidade inferior das placas
utilizadas por Nishino (2007) em seus experimentos
k
y
Componente na direção y da posição de um vórtice k
Letras Gregas
α
Ângulo de orientação do escoamento incidente
m
β
Ângulo de orientação de um painel genérico m
Γ
Intensidade de vórtice
Γ
Circulação
γ
Densidade de vórtice
xii
∆θ
Avanço de um vórtice no intervalo (
2π0
) para o cálculo da
difusão da vorticidade
Avanço de um vórtice na direção radial para o cálculo da difusão
da vorticidade
S
Comprimento de um painel plano
t
Incremento de tempo
CONVECÇÃO
x
Deslocamento de um vórtice na direção x pelo processo de
convecção
DIFUSÃO
x
Deslocamento de um vórtice na direção x pelo processo de difusão
CONVECÇÃO
y
Deslocamento de um vórtice na direção y pelo processo de
convecção
DIFUSÃO
y
Deslocamento de um vórtice na direção y pelo processo de difusão
dδ
B
Espessura da camada limite formada junto ao solo
ε
Valor extra obtido na resolução da equação matricial de vórtices
devido à imposição da conservação da circulação
ς
Constante utilizada no cálculo das cargas aerodinâmicas
θ
Ângulo de orientação de r
µ
Coeficiente de viscosidade dinâmica
ν
Coeficiente de viscosidade cinemática
π
3,14159265359
ρ
Massa específica
Σ
Representa um somatório
0
σ Raio do núcleo do vórtice de Lamb ( coreσ
0
=
)
()
xσ
Distribuição constante de fontes sobre um painel de comprimento
()
xS
Φ
Função potencial de velocidades
Domínio fluido
ω
Campo de vorticidades
Abreviaturas
C.F.D. Dinâmica dos Fluidos Computacional
COUPP Matriz de influência de pressão
xiii
COUPS Matriz de influência - fontes
COUPV Matriz de influência - vórtices
N-S Refere-se às Equações de Navier-Stokes
P.C.M. Princípio de Conservação da Massa
P.C.Q.M.L. Princípio de Conservação da Quantidade de Movimento Linear
RHSP Vetor coluna lado direito de pressão
RHSS Vetor coluna lado direito - fontes
RHSV Vetor coluna lado direito - vórtices
Siglas
UNIFEI Universidade Federal de Itajubá
1
Capítulo 1
INTRODUÇÃO
Todos os corpos presentes em nosso dia-a-dia estão expostos à correntes de ar, de água
ou de um outro tipo de fluido qualquer, sendo que estes corpos são conhecidos na
aerodinâmica como corpos de forma rombuda ou como corpos de forma esbelta. O
escoamento ao redor de corpos rombudos provoca uma grande quantidade de fenômenos
fluido-dinâmicos, como a separação, o desprendimento alternado de pares contra-rotativos de
vórtices e a transição para a turbulência; tais fenômenos despertam grande interesse científico
e têm grande impacto nas aplicações da engenharia. O mecanismo de desprendimento de
vórtices e a sua interrupção tem sido particularmente foco de vários estudos, uma vez que
afeta várias propriedades fluido-dinâmicas como, por exemplo, as forças de arrasto e de
sustentação, vibração, ruído e a eficiência nas transferências de calor e de massa. O
escoamento ao redor de corpos rombudos geralmente é bastante complicado, sendo
necessários esforços experimentais e computacionais para entendê-lo completamente.
Na tentativa de entender fenômenos tão complexos é razoável que se estude o
escoamento ao redor de corpos que possuam uma geometria simples. Dentre eles, os cilindros
de seção circular se apresentam como a melhor alternativa, na medida em que restringem a
complexidade do problema e permitem que se observe as características fundamentais do
escoamento. Na verdade, os cilindros de seção circular têm sido objetos de muitos estudos
não só para se observar as características fundamentais do escoamento ao passar por eles, mas
também pelo fato de terem grande importância em muitas aplicações práticas como, por
exemplo, nos escoamentos ao redor de grandes edifícios, pilares de pontes, chaminés, tubos
2
de trocadores de calor, cabos de transmissão de energia elétrica, entre outros. É interessante
mencionar, ainda, que para os escoamentos ao redor de cilindros de comprimento finito, ou
seja, nos casos práticos, as condições de ponta do cilindro têm influência direta sobre a
configuração final do escoamento.
Toda vez que um fluido viscoso incide sobre um corpo de forma qualquer e conhecida,
verifica-se a condição de aderência sobre a superfície do corpo e a presença de um gradiente
não-nulo de velocidades, o que propicia a formação da camada limite. Em havendo o efeito de
rotação das partículas fluidas, ocorre o fenômeno da geração de vorticidade. Esta vorticidade
se concentra praticamente na região de camada limite e é transportada para jusante, se
concentrando, também, na esteira viscosa do corpo. Especialmente no cilindro circular a
esteira possui um caráter oscilatório quando há um desprendimento alternado de pares contra-
rotativos de vórtices a partir dos pontos de separação; esta esteira é conhecida como esteira de
von Kármán. Entretanto, há muitas situações onde a freqüência de desprendimento dos pares
contra-rotativos de vórtices é reduzida, e uma delas ocorre quando o cilindro circular está
próximo a uma superfície plana (ou solo). À medida que o cilindro de seção circular se
aproxima de uma superfície plana que se move com a mesma velocidade do escoamento
incidente, verifica-se uma diminuição no desprendimento dos pares contra-rotativos de
vórtices gerados a partir dos pontos de separação do corpo. Observa-se, ainda, que à medida
que o cilindro se aproxima do solo, a vorticidade concentrada na esteira do corpo sofre um
efeito de bloqueio, distribuindo-se de forma alongada no domínio fluido, e não mais de forma
aproximadamente circular, como no caso do cilindro circular isolado.
O grande objetivo deste trabalho consiste numa investigação numérica do
comportamento das cargas aerodinâmicas atuantes sobre um cilindro circular que se encontra
na presença do efeito de uma parede plana. Particularmente esta parede plana se move com a
mesma velocidade do escoamento incidente, de forma que não seja criada uma camada limite
a partir da parede plana (Nishino, 2007). Sabe-se que as características aerodinâmicas do
corpo dependem significativamente da distância entre esse corpo e o solo. No entanto, a física
do efeito solo ainda está longe de ser completamente entendida, devido a uma série de fatores
que não podem ser negligenciados. Um destes fatores é a camada limite anteriormente
mencionada, que se forma sobre a superfície plana e se entrelaça com a esteira oriunda do
corpo, dificultando em muito a interpretação física do efeito solo. Muitos testes experimentais
envolvendo um cilindro circular estacionado nas vizinhanças de uma superfície plana fixa
foram feitos, mas devido à dificuldade acima comentada, ainda não existe um consenso sobre
o assunto.
3
O artifício de se utilizar movimento relativo entre o corpo e a superfície plana tem sido
recentemente utilizado no estudo do comportamento aerodinâmico de veículos de alto
desempenho, como os carros de Fórmula 1, mas raramente foi aplicado para o caso do
cilindro circular. Experimentalmente, é difícil usar este artifício, uma vez que são necessários
bons equipamentos para garantir a condição acima citada, além do conhecido alto custo dos
testes em túneis de vento.
A simulação numérica de escoamentos ao redor de corpos com geometrias complexas
(aqui se entendem como geometrias complexas, por exemplo, uma asa se movendo nas
proximidades de uma superfície plana móvel e o movimento das pás no interior de uma voluta
de máquina de fluxo) levando-se em conta os efeitos da viscosidade tem sido basicamente
realizada com a utilização de dois enfoques: descrição lagrangiana (que tem a característica de
acompanhar cada partícula durante todo o tempo de simulação) e descrição euleriana
(caracterizada pela necessidade de se criar uma malha de discretização do domínio fluido a
ser analisado).
Existe uma classe geral de métodos numéricos, os denominados Métodos de Partículas,
que utilizam uma descrição do tipo lagrangiana. Neste tipo de descrição, diferentemente dos
métodos numéricos clássicos (Método de Elementos Finitos, Método de Volumes Finitos,
etc.), preocupa-se com a discretização de uma propriedade dos fluidos. O Método de Vórtices,
certamente o representante mais conhecido dos Métodos de Partículas, discretiza a vorticidade
presente no escoamento representando-a por uma nuvem de vórtices discretos. Para isso, os
vórtices da nuvem são criados a partir das fronteiras sólidas (no interior da camada limite) e
são submetidos aos processos de convecção e difusão da vorticidade, sendo transportados para
o interior da massa fluida. É neste contexto que surge a estratégia numérica do presente
trabalho: uma vez que o fato da velocidade do chão ser igual à velocidade do escoamento
incidente implica na supressão da camada limite junto ao chão, para estudar a física testada
em túnel de vento por Nishino (2007) basta que não sejam gerados vórtices discretos na
superfície plana.
Dessa forma, vê-se que o Método de Vórtices Discretos pode ser empregado em
praticamente todas as aplicações da engenharia que envolvem o movimento de um fluido
(Kamemoto, 2004), especialmente na análise do escoamento ao redor de corpos rombudos ou
esbeltos providos de um razoável ângulo de incidência, onde se observam o fenômeno da
separação e a formação de uma generosa esteira a jusante do corpo. Estes corpos podem estar
4
ainda, animados de um movimento oscilatório ou apresentarem outras fronteiras móveis nas
suas vizinhanças.
Na UNIFEI existe um grupo de professores e de alunos que vem desenvolvendo e
aplicando um Método de Vórtices desde 1998. Os problemas que vêm sendo atacados
envolvem interferência entre corpos com e sem movimento relativo entre eles, efeitos de
oscilação de corpos, fenômenos de transferência de calor e aspectos de turbulência. Quando se
faz referência aos aspectos de transferência de calor, não mais se fala em Método de Vórtices,
mas sim em Métodos de Partículas de Calor (Alcântara Pereira & Hirata, 2003).
Um dos grandes desafios para o desenvolvimento do Método de Vórtices está
relacionado ao elevado tempo de simulação nurica, uma vez que a solução de um problema
utilizando uma nuvem de vórtices discretos envolve o cálculo do campo de velocidades, o
qual é composto por três parcelas: a contribuição do escoamento incidente, a contribuição das
fronteiras sólidas (que pode ser calculada através da utilização do Método de Painéis, o qual
possui uma metodologia baseada na escolha do tipo de singularidade e na escolha do tipo de
condição de contorno) e a contribuição da nuvem de vórtices discretos. Nesta última etapa,
cada vórtice discreto presente no domínio fluido induz velocidade sobre todos os outros
vórtices, sendo a parcela mais onerosa em termos de tempo de CPU, uma vez que o número
de operações realizadas por um processador é proporcional ao quadrado do número de
vórtices presentes na nuvem, se a Lei de Biot-Savart for utilizada. Uma solução alternativa
para o desafio da interação vórtice-vórtice consiste no emprego do Método de Expansão em
Multipolos associado à técnica da computação paralela, o que torna possível a utilização
simultânea de vários processadores durante os cálculos. Devido à situação atual da infra-
estrutura computacional do Grupo de Método de Vórtices da UNIFEI (a paralelização dos
computadores se encontra em estágio inicial de desenvolvimento), este trabalho discretizará a
superfície do cilindro circular em 300 painéis planos, de forma a se obter uma esteira
contendo 300000 vórtices ao final das simulações numéricas utilizando um processador
PENTIUM IV de 1700 MHz.
O momento é oportuno para lembrar que a inclusão dos efeitos da viscosidade no
Método de Vórtices exige, também, cálculos envolvendo métodos refinados para a solução da
difusão da vorticidade. O Método do Crescimento do Raio do Núcleo do Vórtice Modificado
(Rossi, 2006) se apresenta como a alternativa mais atraente para os propósitos do
desenvolvimento do Método de Vórtices na UNIFEI, uma vez que permite a simulação de
escoamentos com baixos valores de número de Reynolds, bem como a repetibilidade dos
5
resultados obtidos. Neste trabalho ainda não se faz uso do referido método, porque ele se
encontra em estágio de implementação.
Os programas computacionais que atualmente vêm simulando os problemas de
Mecânica dos Fluidos estudados no Grupo de Método de Vórtices da UNIFEI são todos
oriundos de uma mesma biblioteca inicial de rotinas, que sofre adaptações em função do
problema analisado. Este trabalho parte de um programa computacional que simulava o
escoamento ao redor de um cilindro circular oscilante e estacionado nas proximidades de uma
superfície plana fixa (Moura, 2007), e exclui a geração de vórtices discretos na superfície
plana, bem como os efeitos harmônicos de oscilação do corpo, de modo a representar a
situação em que um cilindro circular se move em relação a uma superfície plana.
É importante ressaltar que a análise do escoamento ao redor de um cilindro circular que
apresenta movimento em relação a uma superfície plana pode ser extremamente útil para
entender as características essenciais de escoamentos mais complexos, como por exemplo, o
escoamento ao redor de caminhões e de ônibus quando os mesmos trafegam em uma rodovia.
Para que os propósitos deste estudo fossem atingidos, as seguintes etapas foram
efetuadas:
- aferição do novo código computacional implementado, simulando a situação do escoamento
em torno de um cilindro circular isolado;
- cálculo dos coeficientes de arrasto e sustentação médios e análise da evolução temporal
destes coeficientes, bem como o cálculo do número de Strouhal para diferentes distâncias
entre o cilindro e o solo (relação
dh
) quando a superfície plana se move com a mesma
velocidade do escoamento incidente;
- investigação do comportamento das distribuições de pressão média e instantânea em torno
do corpo para diferentes relações
dh ;
- análise da esteira que se forma a jusante do corpo através de animações feitas com a
utilização do software TECPLOT, com destaque para o que acontece nas proximidades do
corpo;
- comparação dos resultados numéricos obtidos com os resultados experimentais disponíveis
na literatura (Blevins, 1984; Roshko
et al., 1975; Nishino, 2007).
6
O Capítulo 2 apresenta uma revisão bibliográfica específica sobre o fenômeno estudado
(incluindo os estudos experimentais), e sobre o Método de Vórtices, contendo importantes
trabalhos realizados com o uso desta ferramenta.
No Capítulo 3 é mostrada a Formulação Geral do Modelo Matemático utilizado e os
principais aspectos que definem a geometria do problema, assim como as hipóteses
simplificadoras e a adimensionalização das equações que governam o escoamento de um
fluido viscoso que incide sobre um corpo na presença do efeito solo (parede plana móvel).
No Capítulo 4 são apresentadas as características marcantes do Método de Vórtices e
efetua-se a solução numérica do problema proposto no Capítulo 3 com o uso do referido
método lagrangiano.
O Capítulo 5 traz um fluxograma que explica a estrutura do programa computacional
desenvolvido em linguagem FORTRAN chamado de “MOVINGROUND.FOR”, bem como a
função que cada rotina executa no presente código.
No Capítulo 6 é feita toda a análise dos resultados numéricos obtidos.
O Capítulo 7 apresenta as conclusões mais importantes tiradas da análise dos resultados
obtidos, e as sugestões para melhorar a implementação numérica em trabalhos futuros.
Na seqüência, encontram-se relacionadas as referências bibliográficas de todos os
trabalhos citados neste texto.
No Apêndice A são discutidos os comportamentos da velocidade tangencial induzida e
da distribuição de vorticidade para o modelo do vórtice potencial e para o modelo do vórtice
de Lamb.
No Apêndice B são apresentados resultados adicionais obtidos a partir das simulações
numéricas realizadas. No Capítulo 6 são apresentadas apenas as análises do efeito solo
(parede plana móvel) para uma relação de distância entre o centro do cilindro circular e a
superfície plana móvel (
dh
).
7
Capítulo 2
REVISÃO BIBLIOGRÁFICA
2.1 – INTRODUÇÃO
O presente capítulo compromete-se em fazer uma revisão atualizada a respeito do
Método de Vórtices para escoamentos incompressíveis, na intenção de situar o leitor no
contexto da ferramenta numérica a ser utilizada no trabalho, bem como uma revisão acerca do
fenômeno do efeito solo, a qual incluí trabalhos experimentais.
2.2 – O MÉTODO DE VÓRTICES DISCRETOS
O Método de Vórtices utilizado na simulação de escoamentos de fluidos viscosos
corresponde a uma abordagem numérica que possui três características fundamentais.
Primeiramente, as equações de Navier-Stokes são formuladas em termos do campo de
vorticidades e não do campo de velocidades. Em segundo lugar, faz-se uso de um dos
teoremas de Helmholtz (Batchelor, 1967; Sherman, 1990), o qual mostra a correspondência
entre os elementos de vorticidade (elementos computacionais ou vórtices discretos) e as
partículas materiais de fluido; a partir desta característica, é possível submeter os vórtices
discretos a um processo convectivo com a mesma velocidade das partículas fluidas, o que
confere uma característica puramente lagrangiana ao método. Por último, para se obter a
8
velocidade do fluido, faz-se uso do fato de que a vorticidade é definida por
uω ×= ; desta
forma, integrando-se o campo de vorticidades, determina-se o campo de velocidades
u . Esta
é a lei de Biot-Savart, que permite descrever completamente o escoamento através do
acompanhamento dos vórtices discretos (Batchelor, 1967).
As características fundamentais descritas acima são as responsáveis tanto pelos aspectos
desejáveis do Método de Vórtices quanto pelas maiores dificuldades encontradas para a sua
implementação numérica. Estudar o escoamento viscoso concentrando-se no campo de
vorticidades é desejável devido ao fato de se propiciar uma melhor visualização do que ocorre
no problema, especialmente em escoamentos vorticosos e em regime não-permanente; além
disso, a equação que rege o transporte da vorticidade na forma bidimensional é escalar. Outra
vantagem é o fato do termo de pressão desaparecer da equação da quantidade de movimento
(equações de Navier-Stokes); posteriormente pode-se, como alternativa, tomar o divergente
das equações de Navier-Stokes para se recuperar, via solução integral oriunda de uma
equação de Poisson para a pressão, o cálculo das cargas aerodinâmicas. Além disso, trabalhar
com o campo de vorticidades sendo discretizado de uma forma lagrangiana faz com que as
condições de contorno no infinito sejam satisfeitas automaticamente, devido à característica
que os vórtices discretos têm de marchar no tempo e de simular a dinâmica da vorticidade. É
importante salientar que o avanço temporal dos vórtices discretos é bastante simples, uma vez
que estas partículas não têm massa. Ao contrário, satisfazer estas condições nos métodos que
utilizam malha pode ser uma tarefa delicada, já que o domínio fluido é limitado. Ademais, os
vórtices discretos são convectados sem que haja dissipação numérica, uma vez que o termo
não-linear das equações de Navier-Stokes é tratado por um conjunto de equações diferenciais
ordinárias que descrevem as trajetórias das partículas; os métodos de malha inevitavelmente
sofrem com essa dissipação numérica, especialmente nos casos de elevados números de
Reynolds. Finalmente, os métodos que não utilizam malha são, na sua essência, mais
vantajosos, já que a geração de malha é uma das etapas mais custosas da dinâmica dos fluidos
computacional (CFD). Por outro lado, o Método de Vórtices Discretos também possui suas
desvantagens, e essas serão discutidas após um breve relato histórico do referido método.
O surgimento do Método de Vórtices é quase tão antigo quanto o Método de Diferenças
Finitas. Entretanto, sua implementação computacional é uma das mais recentes. Na verdade, o
trabalho de Präger (1928) com distribuição de vórtices potenciais sobre uma superfície foi a
origem do Método de Painéis utilizado neste trabalho (Katz & Plotkin, 1991), e largamente
usado na indústria aeronáutica atual.
9
Rosenhead (1931) apresentou o cálculo de folhas de vorticidade utilizando vórtices
discretos potenciais; este trabalho teve tanta importância que é citado até os dias atuais, sendo
um estudo clássico da área. Interessantemente, o Método de Vórtices Discretos Potenciais de
Rosenhead (1931) foi parcialmente desacreditado por volta de 1960 pela observação de que
faltavam provas de que o método convergia (Birkhoff, 1962), e pelos cálculos
computacionais, os quais exibiam um movimento aparentemente caótico das partículas
(Birkhoff & Fisher, 1959). Este último problema foi atribuído à característica singular da
velocidade induzida próximo à origem do vórtice potencial, e diferentes abordagens foram
propostas para solucioná-lo. Uma delas é a utilizada neste trabalho, onde cada vórtice discreto
possui um núcleo viscoso, sendo a vorticidade distribuída de forma gaussiana no interior deste
núcleo (Chorin, 1973), enquanto outras resolvem o problema analiticamente removendo a
singularidade da expressão de Biot-Savart (Krasny, 1983, 1986a, 1986b).
O Método de Vórtices moderno nasceu nos anos de 1970, e os estudiosos proeminentes
foram Alexander Chorin e Anthony Leonard, nos Estados Unidos, e Conrad Rehbach, na
França. O grande interesse do Método de Vórtices durante o começo dos anos de 1980 estava
focado nos aspectos matemáticos, como as propriedades de convergência (Hald & Mauceri,
1978; Hald, 1979; Beale & Majda, 1981, 1982a, 1982b, 1985). Nos últimos anos o
desenvolvimento do método foi muito rico, principalmente com relação às questões da
inclusão correta dos efeitos viscosos, do tratamento das condições de contorno sobre as
fronteiras sólidas e da redução eficiente dos custos computacionais, no intuito de tornar o
método apropriado para realizar simulações de alto nível de escoamentos viscosos em regime
não-permanente e de números de Reynolds elevados. Revisões abrangentes com relação ao
desenvolvimento do Método de Vórtices e suas aplicações podem ser encontradas em
(Leonard, 1980, 1985; Spalart, 1988; Sarpkaya, 1989; Sethian, 1991; Puckett, 1993).
Anderson & Greengard (1991) têm uma coleção de artigos que revelam a conjuntura do
assunto no início dos anos de 1990. Recentemente Cottet & Koumoutsakos (2000) publicaram
um livro dedicado ao assunto, contendo muitas considerações práticas com relação à
implementação do Método de Vórtices. Stock (2007) apresentou uma revisão bastante crítica
sobre a ferramenta numérica utilizada neste trabalho.
Apesar do desenvolvimento contemporâneo com o Método de Diferenças Finitas, o
Método de Vórtices ainda não está entre as ferramentas de CFD mais utilizadas. Ele às vezes
sofre principalmente pela reputação de ser grosseiro na tentativa de modelar escoamentos
viscosos não-permanentes e de alta complexidade, os quais são difíceis de serem tratados
pelos CFD’s tradicionais. Os três principais obstáculos que se apresentam para que o Método
10
de Vórtices seja aceito como uma ferramenta padrão da Dinâmica dos Fluidos Computacional
são: (i) a complexidade numérica para o cálculo da velocidade usando a lei de Biot-Savart,
que requer
2
N operações do processador para N vórtices discretos presentes na nuvem; (ii) o
inconveniente de incluir os efeitos viscosos em uma formulação puramente lagrangiana, uma
vez que a difusão é bem mais facilmente computada nos métodos de malha; e, (iii) o efeito da
evolução lagrangiana no tempo, pois os resultados serão tão precisos quanto menores forem
os incrementos de tempo utilizados; entretanto, quanto menores os incrementos de tempo,
mais onerosas serão as simulações, no que diz respeito ao tempo de CPU. Cabe reforçar
ainda, que o presente trabalho lança mão de uma simulação numérica bidimensional, onde a
evolução temporal do termo de deformação (esticamento) dos tubos de vorticidade é
desprezada; caso este termo fosse considerado, ou seja, caso a simulação numérica fosse
tridimensional, o tempo de CPU seria ainda mais crítico.
A primeira dificuldade mencionada tem sido resolvida com sucesso através do Método
de Expansão em Multipolos (Greengard & Rokhlin, 1987) para o cálculo da interação vórtice-
vórtice. Neste aspecto, o Método de Expansão em Multipolos é uma alternativa à lei de Biot-
Savart no que tange à diminuição dos esforços computacionais. Tal método possibilita que o
processador faça NlogN operações (Interação Vórtice-Caixa) ou até mesmo N operações
para N vórtices da nuvem (Interação Caixa-Caixa), conforme discutido na Tese de Doutorado
de Koumoutsakos (1993). A aplicação do referido método tem sido muito estudada e
implementações eficientes têm sido desenvolvidas (Salmon
et al., 1994), mas devem-se
ressaltar os grandes esforços despendidos na programação do método. Nesta mesma linha de
atuação, alguns pesquisadores têm resolvido o problema usando um sistema de descrição
híbrido (euleriano-lagrangiano), o chamado Método de Vórtices em Célula (Christiansen,
1973; Couet
et al., 1981; Rottman et al., 1987; Cottet, 1987; Smith & Stansby, 1988; Brecht
& Ferrante, 1990; Ebiana & Bartholomew, 1996; Liu & Doorly, 2000; Liu, 2001; Cottet &
Poncet, 2004), o qual tem o ônus de adicionar erros de interpolação e possuir um tempo de
CPU maior do que o Método de Expansão em Multipolos ( MlogMN
+
, sendo M o número
de pontos da malha).
Com relação à inclusão dos efeitos da difusão viscosa, pode ser bastante difícil simulá-
la em um método lagrangiano. Nas últimas três décadas, uma grande quantidade de pesquisas
sobre esse assunto produziu pelo menos seis esquemas diferentes para adicionar a difusão
viscosa nos cálculos do Método de Vórtices, são eles: o Método de Avanço Randômico
(Chorin, 1973; Lewis, 1991), o Método do Crescimento do Raio do Núcleo do Vórtice
(Kamemoto
et al., 1990; Rossi, 1996, 1997, 2005, 2006), o Método da Mudança da
11
Intensidade da Partícula (Degond & Mas-Gallic, 1989), o Método da Redistribuição da
Vorticidade (Shankar & van Dommelen, 1996), o Método da Velocidade de Difusão (Ogami
& Akamatsu, 1991) e o Método de Fishelov (Fishelov, 1990).
O Método de Avanço Randômico, utilizado neste trabalho, foi a primeira técnica
aplicada ao Método de Vórtices para tratar a difusão viscosa (Chorin, 1973); este método é de
simples implementação e de rápida execução, porém possui uma taxa de convergência baixa
(1/
N ). O Método do Crescimento do Raio do Núcleo do Vórtice foi proposto por Leonard
(1980) e utilizado por Kamemoto
et al. (1990); entretanto, tal método foi abandonado quando
Greengard (1985) provou que esta metodologia não convergia para as equações de Navier-
Stokes. Diante disso, Rossi (1996) corrigiu o Método do Crescimento do Raio do Núcleo
fazendo com que os vórtices discretos crescessem até um valor máximo; a partir disto os
vórtices são divididos (partição) dando origem a novos vórtices, cujos raios podem
novamente se expandir. Posteriormente, Rossi (2005, 2006) promoveu modificações para
melhorar o Método do Crescimento do Raio do Núcleo do Vórtice através do uso de Funções
de Base Radial. O Método da Mudança da Intensidade da Partícula foi introduzido por
Degond & Mas-Gallic (1989) e se caracteriza por distribuir a intensidade de um vórtice entre
seus vizinhos através de uma malha. O Método da Redistribuição da Vorticidade
desenvolvido por Shankar & van Dommenlen (1996) é similar ao Método da Mudança da
Intensidade da Partícula, mas dispensa o uso de malha. O Método da Redistribuição da
Vorticidade é resolvido através de um sistema de equações algébricas que otimiza o processo
de redistribuição da intensidade de um vórtice entre seus vizinhos. O Método da Velocidade
de Difusão desenvolvido por Ogami & Akamatsu (1991) simula a difusão da vorticidade
através da inserção de uma velocidade extra no processo convectivo devido exclusivamente à
difusão; tal velocidade é relacionada ao coeficiente de viscosidade cinemática do fluido, ao
campo de vorticidades e ao gradiente do campo de vorticidades. Este método necessitou de
correção do raio do núcleo dos vórtices para que convergisse conforme foi discutido por
Kempka & Strickland (1993). O Método de Fishelov (1990) adiciona derivadas com funções
suaves na equação da vorticidade e transfere tais derivadas para estas funções.
O efeito da evolução lagrangiana no tempo tem sido resolvido pela aplicação de
esquemas que calculam a vorticidade dos vórtices discretos através de uma malha em cada
instante de tempo; estes esquemas utilizam “kernels” de interpolação de alta ordem em uma
formulação cartesiana de produtos tensoriais. Tais esquemas têm sido utilizados há bastante
tempo, e possibilitam o cálculo preciso de escoamentos complexos; no entanto, eles têm
causado controvérsias pelo fato de adicionar uma malha em um método caracterizado por não
12
utilizar malhas. Além disso, eles introduzem alguns erros de interpolação, mas que são
geralmente toleráveis, a menos que se queira simular escoamentos com números de Reynolds
mais altos, quando tais erros podem tornar-se uma limitação. Buscando solucionar tais
inconvenientes, Barba
et al. (2004) utilizaram o Método do Crescimento do Raio do Núcleo
do Vórtice em conjunto com Funções de Base Radial para eliminar a necessidade de
utilização de malha; o controle do crescimento do raio do núcleo dos vórtices discretos é feito
automaticamente e esta técnica se mostrou eficiente quanto à precisão dos resultados obtidos.
No que diz respeito ao cálculo das cargas aerodinâmicas, Lewis (1991) determina o
coeficiente de pressão pela integração do termo de pressão das equações de Navier-Stokes. He
& Su (1999) apresentaram uma nova formulação para a determinação da distribuição de
pressões isolando o termo de pressão das equações de N-S e acrescentando o termo de
aceleração convectiva, não considerado por Lewis (1991). Um estudo de extrema relevância
para a precisão dos resultados é a formulação proposta por Shintani & Akamatsu (1994),
oriunda dos trabalhos de Uhlman (1992) e Kamemoto (1993), que leva em consideração a
contribuição de todos os vórtices presentes na esteira, diferentemente de Lewis (1991) e He &
Su (1999), que consideram apenas a influência dos vórtices nascentes em um dado incremento
de tempo. Buscando obter resultados mais precisos, o presente trabalho utiliza a formulação
proposta por Shintani & Akamatsu (1994).
No que se refere à maior precisão das simulações é essencial destacar a necessidade de
desprendimento de um grande número de vórtices discretos, o que torna o tempo de CPU
bastante oneroso. Portanto, a tendência atual parece ser a procura por algoritmos que tornem
as simulações menos demoradas, como o apresentado por Alcântara Pereira (1999) (Lei de
Biot-Savart Modificada), além da utilização do Método de Expansão em Multipolos
(Greengard & Rokhlin, 1987) para o cálculo da interação vórtice-vórtice, e da expectativa de
desenvolvimento de técnicas como a computação paralela, a qual consiste em conectar
diversas CPU’s entre si para que vários processadores possam trabalhar em conjunto na
simulação (Takeda
et al., 1999).
Na busca por resultados cada vez mais realísticos, Kamemoto
et al. (2000) fazem uma
revisão do Método de Vórtices descrevendo a importância do desenvolvimento de modelos de
turbulência para os métodos lagrangianos. O trabalho de Alcântara Pereira
et al. (2002) foi
preparado visando à realização de simulações numéricas mais refinadas envolvendo os
aspectos de turbulência, tendo como contribuições principais: uma modelagem sub-malha de
turbulência utilizando-se um modelo de Função Estrutura de Velocidade de Segunda Ordem
13
adaptada ao Método de Vórtices Discretos e o desenvolvimento e implementação de um
algoritmo numérico, para incluir, no contexto do Método de Vórtices, a modelagem de
turbulência. Recentemente, Mustto (2004) utilizou com sucesso a formulação proposta por
Alcântara Pereira
et al. (2002) na sua Tese de Doutorado. Neste trabalho não será adotada
modelagem de turbulência, ficando como sugestão para um trabalho futuro a incorporação do
modelo da Função Estrutura de Velocidade de Segunda Ordem ao algoritmo do programa
“MOVINGROUND.FOR”.
Os pesquisadores que desenvolvem as diferentes versões do Método de Vórtices vêm se
reunindo a cada dois anos em uma Conferência Internacional idealizada para discutir os
avanços e as aplicações do método. Esta Conferência já foi realizada em quatro ocasiões:
Kobe, no Japão, em 1999, Istambul, na Turquia, em 2001, Yokohama, no Japão, em 2005 e a
última em Deajeon na Coréia do Sul em abril de 2008.
2.3 – O CILINDRO CIRCULAR ISOLADO
O escoamento ao redor de um cilindro pode ser influenciado pela sua forma e pelo
número de corpos que se apresentam num arranjo no meio fluido. No trabalho de Zdravkovich
(1987) são mostrados vários arranjos de cilindros circulares na presença de um escoamento
incidente que podem ser estudados. Nota-se, portanto, que a forma circular da seção
transversal do cilindro tem um grande número de trabalhos experimentais disponíveis na
literatura, uma vez que o escoamento ao redor de um cilindro circular consiste em um
problema clássico da Mecânica dos Fluidos. A forma circular permite uma combinação ideal
de geometria simples com a configuração complexa do escoamento ao redor de um corpo
rombudo. Essa combinação é extremamente atrativa, uma vez que permite analisar de maneira
simplificada escoamentos ao redor de grandes edifícios, “risers” de plataformas de petróleo e
trocadores de calor, por exemplo.
Conforme mencionado anteriormente, a configuração do escoamento ao redor de um
cilindro é influenciada por uma grande variedade de parâmetros. Para um cilindro de seção
circular longo, submetido a um escoamento uniforme, um parâmetro governante é o número
de Reynolds (
νUdRe = , onde U, d e ν são respectivamente, a velocidade do escoamento
incidente, o diâmetro do cilindro e o coeficiente de viscosidade cinemática). O escoamento ao
14
redor de um cilindro circular pode ser classificado em função do número de Reynolds em
quatro regimes: subcrítico, crítico, supercrítico e transcrítico (Nishino, 2007).
Quando o número de Reynolds é baixo (
1Re
) o escoamento se comporta sem que se
verifiquem fenômenos como separação e predomínio do arrasto de forma; observa-se ainda
uma simetria quase perfeita do escoamento (Figura 2.1a). À medida que o número de
Reynolds aumenta o coeficiente de arrasto sofre uma queda. Para
30Re2 ocorre a
separação da camada limite, o escoamento se apresenta assimétrico e, apesar da presença dos
pares contra-rotativos de vórtices, o regime do escoamento ainda é permanente; a esteira
formada possui um comprimento limitado (Figura 2.1b).
Para
70Re40
começa-se a observar uma oscilação da esteira; a Figura 2.1c mostra
uma situação aproximadamente limite para o regime permanente. Quando o número de
Reynolds se encontra por volta de 90 os pontos de separação não são mais fixos e observa-se
um desprendimento alternado de pares contra-rotativos de vórtices, o que determina o caráter
oscilatório da esteira – a esteira de von Kármán – (Figura 2.1d); neste momento, o arrasto de
pressão (arrasto de forma) é responsável por cerca de 90% do arrasto total. Assim, verifica-se
que existe uma camada limite laminar que evolui para uma camada limite turbulenta. O
escoamento laminar é muito vulnerável ao gradiente adverso de pressão na traseira do cilindro
e a separação ocorre por volta de 82° (Figura 2.1e); a ampla esteira e a pressão muito baixa na
região de separação laminar causam um grande arrasto. Para
54
10Re10 nota-se a
existência de uma generosa esteira pulsante a jusante do cilindro de seção circular. No
entanto, em muitos casos práticos pode haver outros fatores que irão afetar consideravelmente
o escoamento, tais como: a incidência de uma velocidade turbulenta, a rugosidade da
superfície do corpo, a razão de aspecto do cilindro e a superposição de um movimento
oscilatório.
Para o nível de desenvolvimento em que se encontra este trabalho, dentre os principais
estudos que utilizaram o Método de Vórtices para simular o escoamento ao redor de um
cilindro de seção circular está o trabalho de Lewis (1986), que para a época tratava-se de uma
simulação de alto nível, mas que no contexto atual deixa a desejar. Lewis (1986) simulava a
evolução da vorticidade desprendendo dois vórtices discretos por instante de tempo a partir de
pontos de separação pré-determinados sobre o cilindro circular, o que não representava a
física realmente envolvida no fenômeno.
Mustto (1998) estudou o escoamento em torno de um cilindro circular utilizando o
Método de Vórtices associado ao Teorema do Círculo (Milne-Tompson, 1955), o que garantia
15
que a condição de impenetrabilidade fosse satisfeita em todos os pontos da superfície real do
corpo, mas tinha o inconveniente de ser limitado a corpos de forma circular ou
aproximadamente circular (no caso de se lançar mão de uma transformação conforme). A
condição de escorregamento nulo era satisfeita desprendendo-se vórtices sobre pontos pré-
determinados sobre a superfície do cilindro (foram utilizados 16 e 32 pontos de
desprendimento).
(a) 16,0Re
=
(b) 26Re =
(c) 41Re
=
(d) 140Re =
(e)
2000Re
=
(f) 10000Re =
Figura 2.1 – Visualização do escoamento para uma ampla faixa de números de Reynolds.
(Figuras retiradas de Batchelor, 1967; Schlichting, 1979; van Dyke, 1982; Tritton, 1988).
16
Alcântara Pereira (1999) discretizou a superfície de um cilindro circular utilizando o
Método de Painéis (Hess & Smith, 1967), que tem como vantagem discretizar corpos de
forma qualquer e conhecida, mas com o ônus de satisfazer as condições de contorno
(impenetrabilidade e escorregamento nulo) apenas sobre alguns poucos pontos da superfície
discretizada do corpo, os chamados pontos de controle. Além disso, a referida simulação
utilizava um número reduzido de painéis planos (
32M
=
), parâmetro que tem forte influência
na precisão dos resultados. Comparando os resultados de Mustto (1998) e de Alcântara
Pereira (1999) com o resultado experimental de Blevins (1984), verifica-se que ambos os
trabalhos numéricos apresentaram resultados satisfatórios. Os resultados de Mustto (1998)
apresentam uma distribuição de pressão próxima (para
°
<
50θ
) e bem diferente (para
°> 50θ ) da experimental devido ao emprego da formulação de Lewis (1991) para o cálculo
da pressão. A definição do ângulo
θ pode ser encontrada no Capítulo 6.
Dois outros parâmetros que afetam significativamente os resultados são a posição de
geração dos vórtices discretos e o raio do núcleo destes (Koumoutsakos, 1993). Atualmente, a
maneira como esses parâmetros são manipulados no presente trabalho faz com que sejam
inseridos erros nos resultados numéricos, uma vez que o raio do núcleo dos vórtices e as
posições em que estes são gerados são mantidos fixos. Esta metodologia implica em uma
representação errônea dos efeitos viscosos ao longo da superfície do corpo, quando se sabe
que, na verdade, não é dessa forma que a física do problema se processa (Yang & Huang,
1999). Para o regime de Reynolds simulado neste trabalho (
5
10Re = ) existe uma camada
limite laminar que contorna o corpo até certo ponto e que, após isto, ela evolui para uma
camada turbulenta, ficando clara a diferença entre estes dois casos no que tange à intensidade
dos efeitos viscosos presentes no escoamento.
Muitas aplicações práticas de engenharia envolvendo o cilindro circular foram
simuladas utilizando o Método de Vórtices. Dentre as mais recentes está o trabalho de Recicar
(2007), que analisou o escoamento bidimensional, incompressível e em regime não-
permanente ao redor de um cilindro circular submetido a oscilações de grandes amplitudes.
Diferentemente dos trabalhos anteriores encontrados na literatura, Recicar (2007) analisou o
escoamento considerando número de Reynolds alto (
5
10Re = ). Três tipos de regime foram
identificados durante um aumento na freqüência de oscilação do corpo. O primeiro tipo é
observado para baixas freqüências de oscilação do corpo; nesta situação o número de Strouhal
permanece quase constante, correspondendo ao número de Strouhal de um corpo sem
oscilação. O segundo tipo corresponde a um regime de transição, onde aparentemente a
freqüência de emissão de vórtices não se correlaciona com a freqüência de oscilação do corpo.
17
Finalmente, para altas freqüências de oscilação do corpo, a freqüência de emissão de vórtices
coincide com a freqüência de oscilação do corpo, denominada freqüência de “lock-in”.
2.4 – O EFEITO SOLO
As características do escoamento que incide sobre um cilindro circular situado nas
vizinhanças de uma superfície plana (solo) são influenciadas não só pelo número de
Reynolds, mas também pela distância do cilindro ao solo (h). Entretanto, o fenômeno do
efeito solo ainda não está completamente esclarecido, uma vez que existe uma grande
quantidade de fatores que afetam o problema, em especial a camada limite que se forma no
solo e dificulta o entendimento da física envolvida no fenômeno.
Um dos primeiros trabalhos que investigou a influência da distância do cilindro ao solo
foi feito por Taneda (1965). Visualizou-se o escoamento da água ao redor de um cilindro
circular de forma que a água e o solo se moviam em relação ao cilindro, e com a mesma
velocidade (
soloágua
UU
=
). Tais testes foram feitos com número de Reynolds baixo
(
170Re = ). Nestas condições, verificou-se o desprendimento de vórtices do tipo von Kármán
a partir do cilindro para
1,1dh = , enquanto que para 0,6dh
=
, apenas uma modesta fileira
de vórtices foi gerada. Mais detalhes dessa diminuição de geração de vorticidade à medida
que a relação
dh diminui são apresentados nos trabalhos de Zdravkovich (1985a) e Lin et al.
(2005), mas apenas para números de Reynolds baixos, 3350 e 780, respectivamente. Mais
tarde, Roshko
et al. (1975) estudaram experimentalmente o escoamento ao redor de um
cilindro circular estacionado nas proximidades de uma superfície plana fixa em um túnel de
vento, para
4
102,0Re = . Foram observados os comportamentos dos coeficientes de arrasto e
de sustentação, constatando-se que o arrasto diminuía rapidamente, ao passo que a
sustentação aumentava à medida que o cilindro se aproximava do solo.
Bearman & Zdravkovich (1978) investigaram experimentalmente a distribuição do
coeficiente de pressão sobre um cilindro circular estacionado nas vizinhanças de uma
superfície plana fixa para
4
104,8Re = . Constataram que a diminuição do coeficiente de
arrasto para pequenos valores da relação
dh era acompanhada por uma pressão de base alta,
e que o coeficiente de sustentação aumentava quando
dh
decrescia para valores menores do
que 0,9, devido a uma distribuição de pressão assimétrica ao redor do corpo. À medida que o
18
cilindro era afastado da superfície plana fixa (
0,9dh ), a distribuição de pressão se tornava
simétrica e a pressão de base decrescia até
1,5dh
=
. Também se mediu a velocidade nas
regiões próximas ao corpo na tentativa de investigar a freqüência com que vórtices eram
desprendidos na situação de efeito solo. Verificou-se que o número de Strouhal (parâmetro
adimensional utilizado para se medir a freqüência de desprendimento de pares de estruturas
vorticosas contra-rotativas) permanecia aproximadamente constante ( 0,2St
) para qualquer
0,8dh < , e que o pico da medição de velocidade desaparecia para esta faixa de dh .
Portanto, esta distância do cilindro ao solo (
0,8dh
=
) pode ser uma distância crítica, apesar
de que não é sempre que se pode obter uma relação direta entre o desaparecimento do pico da
medição de velocidade e a interrupção do desprendimento de vórtices (Price
et al., 2002).
Buresti & Lanciotti (1979) mediram o número de Strouhal em dois casos distintos: um
cilindro circular rugoso, e outro sem rugosidade, na situação de efeito solo sem movimento
relativo entre o corpo e o chão. Para o cilindro circular não-rugoso imerso em um escoamento
cujo número de Reynolds era de até
5
101,9 , observou-se que a distância crítica entre o
cilindro e o chão era de 0,9, e que o valor do número de Strouhal ficava por volta de 0,2 para
qualquer relação de
dh maior do que 0,9. Assim, vê-se que a distância crítica e o valor do
número de Strouhal dependem do regime do escoamento, sendo impossível definir valores
precisos para tais variáveis; o único consenso parece ser o de que para escoamentos de altos
Re, o número de Strouhal diminui à medida que
dh decresce, mas mesmo assim, este tipo de
comportamento obedece certos limites.
Apesar dos efeitos fundamentais causados pela relação
dh terem sido estudados com
sucesso nos trabalhos acima citados, a influência da camada limite formada junto à superfície
plana fixa é muito complicada e ainda não está totalmente esclarecida.
Zdravkovich (1985b) observou que nos escoamentos com Re variando de
4
104,8 a
5
103,0 o arrasto decrescia rapidamente à medida que a distância do cilindro ao solo fixo
diminuía a um valor menor do que a espessura da camada limite formada no solo (
dδ
B
).
Com isso, concluiu-se que a variação do coeficiente de arrasto era dominada pela relação
B
δh , e não mais pela relação convencional dh . Notou-se ainda, que mesmo o coeficiente
de sustentação poderia ser afetado pelas condições de camada limite, apesar de ser insensível
à espessura desta.
19
Mais tarde, Hiwada
et al. (1986) analisaram o comportamento das forças aerodinâmicas
para
4
102,0Re = e para uma ampla faixa de valores de dδ
B
(de 0,23 a 2,82). Concluiu-se
que o coeficiente de arrasto sofria uma queda à medida que
dh decrescia e que dδ
B
aumentava, o que concordava com os resultados apresentados por Zdravkovich (1985b).
Entretanto, para o caso de menor camada limite (
0,23dδ
B
=
), a diminuição do coeficiente
de arrasto iniciou-se quando a distância do cilindro até o solo era igual a 1,0, caso onde o
cilindro ainda estava fora da camada limite formada no solo, sugerindo que a queda do arrasto
poderia ser causada não só pela interferência direta da camada limite, como havia se pensado.
Com relação ao coeficiente de sustentação, resultados similares aos de Roshko
et al. (1975)
foram obtidos para todos os valores de
dδ
B
investigados, ou seja, a sustentação não sofre
influência da camada limite formada junto ao solo. Hiwada
et al. (1986) ainda mediram a
velocidade na região próxima ao corpo e encontraram
0,8dh
=
como sendo o valor crítico
da distância do cilindro ao chão, para todas as relações
dδ
B
testadas. O número de Strouhal,
no entanto, diminuía à medida que
dh decrescia e que dδ
B
aumentava, contrariando os
resultados obtidos por Bearman & Zdravkovich (1978), para
4
104,8Re = .
Outro estudo relatando os efeitos da relação
dδ
B
foi feito por Taniguchi & Miyakoshi
(1990) para
4
109,4Re = . Ao contrário dos resultados obtidos por Hiwada et al. (1986),
verificou-se que a distância crítica
(
)
c
dh tornava-se mais ampla (de 0,8 para 1,4) à medida
que
dδ
B
aumentava (de 0,34 para 1,05). A relação entre
(
)
c
dh e dδ
B
não era linear, mas
exponencial. Neste trabalho, os autores obtiveram os valores rms dos coeficientes de arrasto e
de sustentação (
'
D
C e
'
L
C ), e identificaram que a relação
(
)
c
dh correspondia a um valor de
dh onde os valores de
'
D
C e
'
L
C eram mínimos. A relação entre o valor de
'
L
C e o início ou
interrupção do desprendimento de vórtices foi confirmada por Lei
et al. (1999). Neste último
estudo, entretanto, a relação
()
c
dh diminuía suavemente (de 0,9 para 0,8) à medida que
dδ
B
aumentava (de 0,14 para 2,89), contrariando as observações de Taniguchi & Miyakoshi
(1990).
Zdravkovich (2003) estudou o comportamento do coeficiente de arrasto em um cilindro
circular localizado próximo a uma superfície plana móvel; a superfície se movimentava com a
mesma velocidade do escoamento incidente. O número de Reynolds do escoamento era alto
(
5
102,5Re = ) e, nas condições descritas acima, praticamente não houve formação de
camada limite junto à superfície plana móvel. Ao contrário de todos os estudos comentados
20
anteriormente, o autor não verificou a queda do arrasto quando
dh diminuía. Entretanto, não
ficou claro se o fato observado ocorreu devido ao alto número de Reynolds ou à inexistência
da camada limite junto ao solo, ou ainda, devido a algum outro fator influente.
Recentemente, Nishino (2007) estudou o comportamento dos coeficientes de arrasto e
de sustentação atuantes sobre um cilindro circular localizado próximo a uma esteira rolante
em um túnel de vento. A esteira rolante se movimentava com a mesma velocidade do
escoamento incidente, ou seja, não houve formação de camada limite junto ao solo. Foi
investigado o comportamento aerodinâmico do corpo para
4
104,0Re = e
5
101,0Re = . Os
efeitos de ponta do cilindro também foram estudados, mediante o uso de placas nas suas
extremidades. Para o cilindro sem placas (caso essencialmente tridimensional) verificou-se
que o coeficiente de arrasto aumentava à medida que
dh diminuía, e que isto ocorria devido
à inexistência de camada limite no chão, esclarecendo a dúvida deixada no trabalho de
Zdravkovich (2003). Com relação às placas, à medida que a distância do cilindro até as
extremidades das placas aumentava, retirava-se parte da tridimensionalidade do problema;
para o caso mais bidimensional estudado, o valor obtido para o coeficiente de arrasto foi o de
1,3, quando
2,50dh = , o que é bem próximo do valor obtido para o caso de um cilindro
circular isolado para
5
101,0Re = (Blevins, 1984). Quando o cilindro era posicionado de
forma a tangenciar as extremidades das placas ((
0dy
e
=
), caso menos bidimensional
estudado), verificou-se que para
0,85dh
<
as placas não exerciam influência alguma sobre o
coeficiente de arrasto. Assim, obteve-se um valor de
D
C aproximadamente constante e igual a
0,95. No que diz respeito ao coeficiente de sustentação, a situação de movimento relativo com
a velocidade do escoamento incidente igual à velocidade da esteira rolante não exerceu
modificação sobre este parâmetro, em comparação com o caso de corpo e esteira rolante
parados. Sendo assim, sempre que
dh diminuiu observou-se que o
L
C aumentou, tal como
observado por Roshko et al. (1975). A Figura 2.2 ilustra algumas das análises anteriormente
discutidas, onde a expressão “Present study” se refere ao trabalho experimental de Nishino
(2007). Ainda nesta figura, a relação
dy
e
corresponde à distância entre o cilindro de seção
circular e as extremidades inferiores das placas.
É importante ressaltar que para efeito de comparação com os resultados obtidos neste
trabalho, os valores da relação
dh
testados pelos autores anteriormente referidos estão todos
acrescidos do valor 0,5, ou seja, a relação
dh corresponde à distância entre o centro do
cilindro de seção circular e o solo.
21
(a) Coeficiente de Arrasto (b) Coeficiente de Sustentação
Figura 2.2 – Análise das cargas aerodinâmicas atuantes sobre um cilindro circular na presença
do efeito solo. (Figura retirada do trabalho de Nishino (2007)).
Com relação à ferramenta numérica em desenvolvimento que foi utilizada neste
trabalho, o Método de Vórtices, alguns trabalhos envolvendo o fenômeno do efeito solo (sem
movimento da superfície plana) merecem destaque.
Chacaltana
et al. (1994) consideraram os efeitos da presença de um vórtice livre
movendo-se nas vizinhanças de um aerofólio, ao mesmo tempo em que consideraram o efeito
do solo em um mesmo escoamento. O efeito solo, para aerofólios delgados, foi bem
analisado, porém sem considerar a convecção do vórtice. A modelagem do problema era
potencial, restrita a corpos delgados. A formação da esteira era levada em conta com a difusão
de vorticidade a partir do bordo de fuga do aerofólio. Não foi prevista a geração de vórtices
junto à superfície do solo.
Chacaltana
et al. (1995) adotaram o mesmo modelo do trabalho anterior, porém
permitindo ao vórtice se deslocar por convecção, utilizando um esquema lagrangiano de
avanço no tempo. O problema era formulado sob o ponto de vista potencial. Este trabalho foi
uma evolução ao trabalho anterior, cabendo as mesmas observações. Nota-se que,
provavelmente, em função da ausência dos efeitos viscosos junto ao solo, os resultados destes
trabalhos indicam uma diminuição da sustentação para uma aproximação à superfície plana.
Ricci (2002) utilizou o Método de Vórtices para analisar o escoamento de um fluido
Newtoniano e homogêneo em torno de um corpo de forma arbitrária, quando disposto nas
proximidades de uma superfície plana fixa. A superfície do solo foi simulada com a utilização
do Método de Imagens, no qual a condição de impenetrabilidade era satisfeita
automaticamente em todos os pontos da superfície, porém tinha o ônus de demandar maior
22
tempo de simulação, uma vez que os vórtices-imagens também deviam ser considerados no
cálculo da velocidade total induzida nos vórtices discretos. Vale lembrar que a velocidade
nunca deve ser calculada sobre os vórtices-imagens.
Silva de Oliveira
et al. (2005) simularam a interferência da esteira de um perfil NACA
0012 com a esteira formada a partir de uma superfície plana fixa e ainda com a esteira oriunda
de uma geração contínua de vórtices numa fila que se movia com velocidade constante e
ficava posicionada a montante do corpo. Neste trabalho foi investigado o comportamento das
cargas aerodinâmicas e do campo de velocidades nas proximidades do bordo de ataque do
perfil. Diferentemente de Ricci (2002), Silva de Oliveira
et al. (2005) representaram a
superfície do solo mediante a utilização do Método de Painéis (Hess & Smith, 1967).
Recentemente, Moura (2007) utilizou a mesma metodologia de Silva de Oliveira
et al.
(2005) e analisou o escoamento bidimensional, incompressível e em regime não-permanente
de um fluido Newtoniano com propriedades constantes que incidia sobre um cilindro circular
na presença do efeito solo sem a fila móvel composta por vórtices discretos. O corpo
apresentava um movimento de oscilação linear e de pequena amplitude na direção
perpendicular ao escoamento incidente. Estudou-se a influência que o movimento oscilatório
de pequena amplitude exercia sobre o coeficiente de sustentação na presença do efeito solo.
O presente trabalho vem contribuir no contexto do Método de Vórtices, com o
aperfeiçoamento de uma ferramenta lagrangiana em desenvolvimento. Foi investigado o
fenômeno do efeito solo em uma situação de grande aplicação prática da engenharia: a
situação de movimento relativo. A simulação numérica do escoamento bidimensional,
incompressível e em regime não-permanente de um fluido Newtoniano com propriedades
constantes que incide sobre um cilindro circular estacionado próximo a uma parede plana, foi
feita gerando-se vórtices discretos de Lamb apenas sobre a superfície do corpo, a qual foi
discretizada em painéis planos, sobre os quais foram distribuídas singularidades do tipo fontes
com densidade constante. Como não há geração de vórtices na superfície plana, tal situação
equivale àquelas testadas em túneis de vento por Zdravkovich (2003) e Nishino (2007).
Detalhes da simulação numérica realizada neste trabalho podem ser vistos no Capítulo 5. Ao
contrário dos trabalhos de Ricci (2002) e Moura (2007), o presente estudo não se limitou
apenas à análise do comportamento do coeficiente de sustentação; o coeficiente de arrasto
também foi investigado e comparado com os resultados obtidos experimentalmente por
Nishino (2007), para o caso mais bidimensional estudado pelo citado autor.
23
Capítulo 3
FORMULAÇÃO GERAL DO PROBLEMA
3.1 - INTRODUÇÃO
Este capítulo apresenta a formulação matemática necessária para a realização da
simulação numérica do escoamento bidimensional, incompressível e em regime não-
permanente de um fluido Newtoniano com propriedades constantes, que incide sobre um
corpo de forma qualquer e conhecida utilizando o Método de Vórtices Discretos. Este corpo
se encontra na presença de uma superfície plana, a qual se move com a mesma velocidade do
escoamento incidente; sendo assim, não há formação de camada limite junto ao solo. Para
cumprir esta proposta devem-se especificar as hipóteses simplificadoras de forma a garantir
que as equações que governam o problema traduzam a física do caso em estudo.
3.2 - GEOMETRIA E DEFINIÇÕES
Uma representação esquemática do escoamento descrito anteriormente pode ser vista na
Figura 3.1, a qual ilustra um cilindro circular imerso numa região fluida bidimensional e
semi-infinita. Ainda, pode-se identificar a presença de uma superfície plana se movimentando
em relação ao corpo com a mesma velocidade do escoamento incidente.
24
Figura 3.1 – Definição da geometria e da região fluida do problema com o sistema de
coordenadas fixo ao corpo.
Na figura acima são definidos:
=
321
SSSS fronteira que delimita o domínio fluido semi-infinito;
1
S
contorno da superfície do cilindro circular;
2
S contorno da superfície plana móvel;
3
S contorno da superfície definida a grandes distâncias do corpo, definida para
r
;
d
diâmetro do cilindro de seção circular (comprimento característico do problema);
U velocidade do escoamento incidente (velocidade característica do problema);
α ângulo de orientação do escoamento incidente;
h distância do centro do cilindro circular até a superfície plana móvel.
25
3.3 - HIPÓTESES SIMPLIFICADORAS
De maneira geral, os problemas de interesse prático da engenharia, como o analisado
neste trabalho, são de difícil solução. O problema em particular analisado aqui depara-se com
a presença do campo de vorticidades superposto ao campo de velocidades; a vorticidade
presente no domínio fluido não permite que se encontre uma solução analítica. Assim, na
solução de tais problemas deve-se, também, fazer uma análise cuidadosa do caso em estudo
para se lançar mão de hipóteses que simplifiquem as equações governantes sem que haja uma
perda de generalidade do problema. No presente trabalho, as hipóteses que sustentam a
formulação matemática relacionam-se à geometria do corpo, às propriedades dos fluidos e às
propriedades do escoamento, sendo mais bem explicadas abaixo:
a) Hipótese relativa à geometria
1
H região fluida semi-infinita e bidimensional: já que a superfície
3
S foi definida a
grandes distâncias do corpo, conforme visto na Figura 3.1.
b) Hipótese relativa às propriedades dos fluidos
2
H fluido Newtoniano: caracterizado pela existência de uma relação linear entre a tensão
tangencial aplicada e a respectiva taxa de deformação sofrida; além disso, admite-se que as
propriedades do fluido sejam constantes, ou seja, a massa específica e o coeficiente de
viscosidade mantêm-se invariáveis.
Esta última hipótese relativa às propriedades do fluido faz com que a equação do
movimento (Princípio da Conservação da Quantidade de Movimento Linear) seja
representada pelas equações de Navier-Stokes; veja a seguir.
c) Hipóteses relativas às propriedades do escoamento
3
H Escoamento incompressível: os efeitos de compressibilidade do escoamento são
desprezados. Normalmente, adota-se 0,3Ma
<
;
4
H Escoamento laminar: embora o número de Reynolds seja elevado, esta hipótese é
adotada devido à ausência de um modelo de turbulência. O modelo de turbulência da Função
Estrutura de Velocidade de Ordem 2 (Alcântara Pereira
et al., 2002) será futuramente
incorporado nesta formulação.
26
A presença dos efeitos da viscosidade dentro da camada limite faz com que o perfil do
campo de velocidades seja não retangular, provendo as partículas fluidas de um movimento
de rotação; conseqüentemente, próximo de fronteiras sólidas tem-se o fenômeno da geração
de vorticidade (veja mais detalhes no Capítulo 4). A vorticidade presente no escoamento, no
contexto deste trabalho, fica representada por uma nuvem de vórtices discretos de Lamb
(Panton, 1984) e o escoamento é rotacional na camada limite criada ao redor da superfície do
corpo, nas camadas cisalhantes surgidas dos pontos de separação e na esteira viscosa formada
a jusante do corpo.
As hipóteses anteriormente apresentadas, seguidas das equações governantes, bem
como das condições de contorno que definem o problema, permitem a formulação matemática
do caso em estudo.
3.4 - EQUAÇÕES GOVERNANTES E CONDIÇÕES DE
CONTORNO
3.4.1 – Equações Governantes
O fenômeno a ser estudado é governado pelas expressões matemáticas que representam
os Princípios de Conservação (Princípio de Conservação da Massa e Princípio de
Conservação da Quantidade de Movimento Linear). De acordo com as hipóteses
simplificadoras anteriormente apresentadas, pode-se escrever:
a) Princípio de Conservação da Massa (P.C.M.)
Fazendo-se um balanço de massa em um volume de controle elementar, cartesiano e
fixo, verifica-se que o divergente do campo de velocidades é nulo e, na forma diferencial, a
Equação da Continuidade é expressa por:
0=u (3.1)
onde
u representa o campo de velocidades do escoamento.
b) Princípio de Conservação da Quantidade de Movimento Linear (P.C.Q.M.L.)
27
Aplicando-se a Segunda Lei de Newton sobre um volume de controle elementar,
cartesiano e fixo, e desprezando-se os efeitos de força de campo gravitacional, obtém-se a
seguinte forma para as Equações de Navier-Stokes (N-S):
uuu
u
2
ρ
µ
p
ρ
1
t
+=+
(3.2)
onde, p é o campo de pressões, ρ é a massa específica e
µ
é o coeficiente de viscosidade. O
Princípio de Conservação da Energia não foi considerado devido à utilização da hipótese de
escoamento isotérmico.
3.4.2 – Condições de Contorno
Considerando que a superfície
1
S do corpo seja não porosa, pode-se através da análise
da Figura 3.2, entender as duas condições de contorno que serão apresentadas.
Figura 3.2 – Especificação das condições de contorno sobre as fronteiras que delimitam o
domínio fluido.
onde:
τ
e vetor unitário tangente às superfícies
1
S
e
2
S
em cada ponto;
n
e vetor unitário normal às superfícies
1
S e
2
S em cada ponto.
28
A condição de impenetrabilidade, suficiente no Modelo Potencial, exige que o
componente normal da velocidade da partícula de fluido (
n
u ) seja igual ao componente
normal das velocidades das superfícies
1
S e
2
S (
n
v ), ou seja:
0v-u
nn
= , em
1
S e
2
S (3.3)
No entanto, o Modelo Viscoso (Newtoniano) exige, além da condição de
impenetrabilidade definida pelo Modelo Potencial, também a especificação da condição de
escorregamento nulo, a qual exige que o componente tangencial da velocidade da partícula de
fluido (
τ
u ) seja igual ao componente tangencial da velocidade da superfície
1
S
(
τ
v ), ou seja:
0v-u
ττ
= , apenas em
1
S
(3.4)
uma vez que a superfície plana (
2
S ) se move com a mesma velocidade do escoamento
incidente para que não haja formação de camada limite junto ao solo (Nishino, 2007). A
conseqüência de se impor a condição de escorregamento nulo apenas sobre a superfície do
corpo reside no fato de que vórtices discretos de Lamb serão gerados apenas no corpo. Esta
metodologia permite representar o movimento da superfície plana.
Na verdade, as condições de impenetrabilidade e de escorregamento nulo podem ser
expressas em uma única condição – a condição de aderência – a qual exige que as partículas
fluidas não migrem para o interior da superfície sólida e nem possuam velocidade relativa a
esta.
Por fim, a última condição de contorno a ser imposta refere-se à condição de
escoamento não-perturbado, ou seja, a grandes distâncias do corpo e do solo o escoamento em
estudo deve tender para o valor da velocidade do escoamento incidente, isto é:
3
S em ,Uu
(3.5)
3.5 – ADIMENSIONALIZAÇÃO DO PROBLEMA
A técnica da Análise Dimensional consiste em simplificar (e reduzir) a quantidade de
variáveis necessárias para se descrever um dado fenômeno físico. As variáveis obtidas devem
reter todas as informações importantes do fenômeno analisado. Tal técnica é extremamente
29
usada em vários campos da engenharia, uma vez que projetos de natureza específica
normalmente são caros e de construção complexa. Assim, a adimensionalização do problema
permite que os resultados do programa computacional sejam comparados com resultados de
ensaios experimentais sem a necessidade de se ter conhecimento dos valores das grandezas
envolvidas no experimento, bastando-se conhecer parâmetros adimensionais afetos ao
fenômeno analisado (White, 2002).
Para que se possa realizar a adimensionalização do problema em estudo, devem-se
adotar as seguintes escalas:
d escala de comprimento; neste caso, o diâmetro do cilindro circular;
U
escala de velocidade; neste caso, a velocidade do escoamento incidente;
U
d
escala de tempo.
Definidas as escalas representativas do problema, a adimensionalização das demais
grandezas é imediata, como se segue:
d
x
x
*
= , que representa a medida da abscissa x;
d
y
y
*
= , que representa a medida da ordenada y;
= d
*
, que representa o operador “Nabla”;
=
2*2
d , que representa o operador Laplaciano;
U
u
u
*
= , que representa o componente da velocidade na direção x;
U
v
v
*
= , que representa o componente da velocidade na direção y;
d
tU
t
*
= , que representa um instante de tempo da simulação numérica;
µ
ρUd
Re = , que representa o número de Reynolds;
30
2
*
ρU
p
p = , que representa o campo de pressões;
dU
Γ
Γ
*
=
, que representa a intensidade de vórtice;
dU
σ
σ
*
= , que representa a densidade de fontes;
d
σ
σ
0
*
0
= , que representa o raio do núcleo do vórtice de Lamb;
U
ωd
ω
*
=
, que representa o único componente não-nulo do vetor vorticidade no plano;
U
d
fSt = , que representa o número de Strouhal.
O significado de algumas grandezas que aparecem acima definidas ficará mais bem
compreendido com o desenvolvimento do texto. Assim, as equações adimensionalizadas
tomam a forma:
0 :P.C.M.
**
= u (3.6)
*2******
*
*
Re
1
p
t
:P.C.Q.M.L.
uuu
u
+−∇=+
(3.7)
E as condições de contorno, ficam:
0vu
*
n
*
n
= , em
1
S e
2
S (3.8)
0vu
*
τ
*
τ
=
, apenas em
1
S (3.9)
3
*
S em ,1u (3.10)
O asterisco (*), que designa grandeza adimensionalizada, será agora omitido por
comodidade de digitação.
A Figura 3.3 ilustra o problema adimensionalizado.
31
Figura 3.3 – Representação do problema adimensionalizado.
Embora a superfície plana apareça com uma velocidade na Figura 3.3, vale lembrar
novamente que a simulação da situação de movimento relativo no presente trabalho se
caracteriza pelo simples fato de não se realizar a geração de vórtices discretos junto à
superfície. Esta estratégia numérica inibe o desenvolvimento de camada limite no solo, sendo
de simples implementação numérica.
3.6 – EQUAÇÃO DO TRANSPORTE DA VORTICIDADE
Tomando-se o rotacional nas equações de Navier-Stokes e com o auxílio da Equação da
Continuidade, obtém-se a Equação do Transporte da Vorticidade (Batchelor, 1967):
()( )
ωuωωu
ω
2
Re
1
t
+=+
(3.11)
onde:
t
ω
taxa de variação local da vorticidade;
()
ωu taxa de variação convectiva da vorticidade;
32
()
uω taxa de deformação dos tubos de vorticidade; este termo só é aplicado à
escoamentos tridimensionais;
ω
2
Re
1
taxa de variação difusiva da vorticidade.
Conforme dito anteriormente, o caso em estudo trata-se de um escoamento
bidimensional. Portanto, a Equação do Transporte da Vorticidade para este caso (uma
equação escalar) toma a forma:
()
ω
Re
1
ω
t
ω
2
=+
u (3.12)
Chorin (1973) através do “Viscous Splitting Algorithm” (Algoritmo de Separação da
Parte Viscosa da Equação do Transporte da Vorticidade), propôs um algoritmo que em um
mesmo incremento de tempo os efeitos convectivos são calculados independentemente dos
efeitos difusivos, simplificando assim a implementação numérica do Método de Vórtices, o
qual é caracterizado por acompanhar todas as partículas de fluido ao longo do tempo
(descrição lagrangiana).
Dessa forma, a equação da convecção da vorticidade fica:
()
0ω
t
ω
D
t
Dω
=+
=
u (3.13)
a qual será resolvida utilizando-se a sua versão lagrangiana (
0DtDω = ), uma vez que o
termo convectivo
()
ωu , que é não-linear, é de difícil solução.
Por outro lado, a equação da difusão da vorticidade é dada por:
ω
Re
1
t
ω
2
=
(3.14)
Como mostra a presença do número de Reynolds, o processo de difusão da vorticidade
leva em conta os efeitos da viscosidade.
Quando o intervalo de tempo
t
tender a zero, as Equações 3.13 e 3.14 tenderão para a
Equação 3.12.
33
Capítulo 4
MÉTODO DE SOLUÇÃO: O MÉTODO DE VÓRTICES
4.1 – INTRODUÇÃO
A presença dos efeitos viscosos no problema formulado no Capítulo 3 não permite que
se obtenha solução analítica; assim, este capítulo apresenta uma explanação acerca das
potencialidades do Método de Vórtices Discretos, bem como uma solução numérica para o
problema proposto no capítulo anterior via Método de Vórtices.
4.2 – CARACTERÍSTICAS MARCANTES DO MÉTODO DE
VÓRTICES
O Método de Vórtices consiste em uma ferramenta numérica de grande importância
para simulações de escoamentos de fluidos viscosos ao redor de geometrias complexas com
ou sem movimento relativo entre elas (como exemplos de geometrias complexas aqui se
entendem: o escoamento ao redor de um arranjo de cilindros, o escoamento ao redor de uma
asa se movendo nas proximidades de uma superfície plana fixa, o escoamento ao redor das
pás no interior de uma voluta de máquina de fluxo, etc.), especialmente quando tais
escoamentos são dependentes do tempo (presença de um campo de vorticidades).
34
A simulação numérica do escoamento de um fluido viscoso ao redor de um corpo pode
ser feita basicamente com a utilização de dois enfoques: métodos que utilizam a descrição
euleriana e métodos que utilizam a descrição lagrangiana.
A descrição euleriana caracteriza-se pela necessidade de geração de uma malha no
domínio fluido a ser analisado e associar a cada malha as grandezas de interesse para o
escoamento. Tal característica, no entanto, faz com que as grandezas que governam o
escoamento sejam calculadas em pontos onde estas não são importantes; por outro lado, nas
regiões de interesse deve-se promover um refinamento da malha para que sejam obtidos
resultados mais precisos.
O Método de Vórtices Discretos utiliza a descrição lagrangiana para simular o
escoamento e tem, como principal característica, o fato de acompanhar cada vórtice discreto
gerado durante todo o tempo de simulação. Isto faz com que as condições de contorno no
infinito, ao contrário da descrição euleriana, sejam satisfeitas automaticamente. Tal método é
apropriado para simular o escoamento ao redor de corpos rombudos ou esbeltos providos de
um considerável ângulo de ataque, corpos estes que podem estar em repouso, em movimento,
ou ainda, sofrendo oscilações, uma vez que há descolamento de camada limite e formação de
uma esteira viscosa a jusante destes corpos.
No presente trabalho, estuda-se o escoamento ao redor de um corpo rombudo que se
encontra nas proximidades de uma superfície plana móvel, uma vez que se deseja investigar o
fenômeno do efeito solo sem a presença da vorticidade gerada na superfície do solo. Como se
faz uso do Método de Vórtices, o escoamento será analisado pela dinâmica da vorticidade e,
consequentemente, pela generosa esteira que se forma a jusante do corpo.
No contexto do Método de Vórtices, o cálculo do campo de velocidades é feito apenas
nas posições ocupadas pelos vórtices discretos sendo composto pelas seguintes contribuições:
- Contribuição do escoamento incidente: este cálculo não apresenta dificuldades, porque o
escoamento incidente induz a mesma velocidade em todos os vórtices discretos da nuvem;
- Contribuição das fronteiras sólidas: as fronteiras sólidas presentes em um problema também
induzem velocidade na nuvem de vórtices; no caso deste trabalho estas contribuições são
devidas às presenças do cilindro circular e da superfície plana móvel;
- Contribuição da nuvem: cada vórtice discreto induz, em cada incremento de tempo,
velocidade em todos os outros vórtices presentes na nuvem. Aqui fica evidente o alto custo
35
computacional empregado em tais simulações, uma vez que o tempo de simulação é
proporcional ao quadrado do número de vórtices. A estrutura computacional disponível para
as simulações numéricas deste trabalho é inviável para que se possa pensar em uma nuvem de
vórtices com mais de 300.000 vórtices discretos. O grau de refinamento conseguido neste
trabalho só foi possível devido à não geração de vórtices discretos a partir do solo; assim, toda
a carga computacional ficou concentrada nos vórtices gerados a partir do corpo. Conforme
será discutido no Capítulo 6, o tempo médio de CPU para cada experimento numérico
utilizando-se um processador PENTIUM IV de 1700 MHz fica em torno de 133 horas
ininterruptas.
No decorrer deste capítulo serão apresentadas discussões sobre o fenômeno da geração
de vorticidade, bem como sobre os problemas da convecção e da difusão desta grandeza
cinemática. Além disso, será apresentada a formulação proposta por Shintani & Akamatsu
(1994) para o cálculo das cargas aerodinâmicas.
4.3 – A CONVECÇÃO DA VORTICIDADE
Com a utilização do “Viscous Splitting Algorithm” proposto por Chorin (1973) os
efeitos convectivos são calculados independentemente dos efeitos difusivos, em cada
incremento de tempo.
Sabe-se que a etapa convectiva independe da ação da viscosidade, o que caracteriza o
escoamento nesta fase como sendo potencial (veja a Equação 3.13).
Neste trabalho utilizam-se singularidades do tipo fontes, impondo-se a condição de
impenetrabilidade sobre a superfície do corpo e sobre a superfície plana móvel (Equação 3.3).
A condição de conservação global da circulação é imposta ao longo do domínio fluido
, e é
expressa por:
==
c
0Γ dSu (4.1)
4.3.1 – Contribuição do Escoamento Incidente
36
Conforme dito na Seção 4.2, o escoamento incidente induz velocidade nos vórtices
discretos de Lamb. Este é dado nas direções x e y, respectivamente, por (veja a Figura 4.1):
1cosαUcosαu =
=
=
(4.2)
0senαUsenαv =
=
=
(4.3)
4.3.2 – Contribuição do Corpo: O Método dos Painéis
O Método dos Painéis vem a contribuir de forma fundamental na solução do problema
potencial no contexto do Método de Vórtices.
Figura 4.1 – Condições de Contorno.
O Método dos Painéis (Katz & Plotkin, 1991) possibilita o estudo do escoamento ao
redor de corpos das mais diversas geometrias. Este método tem como característica marcante
a discretização da superfície de um corpo de forma qualquer e conhecida em segmentos
planos (os painéis retos) ou curvos (os painéis curvos) e a distribuição de singularidades sobre
tais segmentos de forma a garantir a validade das condições de contorno em cada ponto de
controle dos painéis. Aliás, essa é uma das desvantagens do Método dos Painéis: garantir as
condições de contorno em apenas um único ponto de cada painel, e não ao longo de todo o
comprimento do painel; deve-se ressaltar ainda, que o polígono formado pelos painéis
37
representa uma aproximação de uma superfície sólida real e que, neste trabalho, os painéis
têm comprimentos iguais.
No que concerne à escolha do tipo de singularidade, pode-se optar por uma distribuição
com densidade concentrada, constante ou linear ao longo de cada painel. No presente trabalho
serão utilizados painéis planos, sobre os quais se distribuíram singularidades do tipo fontes
com densidade constante (ou uniforme), como mostra a Figura 4.2 (Katz & Plotkin, 1991).
Sendo assim, neste estudo, impõe-se a condição de velocidade normal nula, ou seja, a
condição de impenetrabilidade (condição de contorno de Neumann) sobre as superfícies
1
S e
2
S .
Figura 4.2 – Distribuição constante de fontes sobre o eixo x.
Os componentes na direção de x e de y da velocidade induzida no ponto
()
yx,P (Figura
4.3) pela distribuição de fontes com densidade constante,
(
)
xσ , ao longo de uma superfície de
comprimento
()
12
xx
valem, respectivamente:
()
()
()()
+
=
2
1
x
x
0
2
0
2
0
0
dx
yyxx
xx
2π
xσ
u (4.4)
()
()
()()
+
=
2
1
x
x
0
2
0
2
0
0
dx
yyxx
yy
2π
xσ
v (4.5)
38
Figura 4.3 – Velocidade induzida no ponto P(x,y) por uma distribuição de fontes com
densidade constante,
()
xσ
, ao longo de uma superfície de comprimento
()
12
xx
.
Resolvendo-se as duas integrais acima (Katz & Plotkin, 1991), obtém-se:
() ()
2
2
2
1
2
1
r
r
ln
4π
xσ
r
r
ln
2π
xσ
u == (4.6)
()
()
12
θθ
2π
xσ
v = (4.7)
onde:
2 1,k ,
xx
y
tgθ
k
1
k
=
=
(4.8)
()
2 1, k , yxxr
2
2
kk
=+= (4.9)
A indução de velocidades na direção de x para
±
0y de um painel sobre ele mesmo,
considerando a Figura 4.3, é dada pela equação:
()
2
xσ
,0
2
xx
v
12
±=
±
(4.10)
Quando aplicada nos M pontos de controle dos painéis (que incluem o corpo e o solo), a
condição de contorno de Neumann pode ser expressa na seguinte forma matricial:
39
=
M
3
2
1
M
3
2
1
s
M3
s
M2
s
M1
s
3M
s
32
s
31
s
2M
s
23
s
21
s
1M
s
13
s
12
RHSS
.
.
.
RHSS
RHSS
RHSS
σ
.
.
.
σ
σ
σ
0,5...KKK
.......
.......
.......
K...0,5KK
K...K0,5K
K...KK0,5
(4.11)
onde:
s
ij
K
é um elemento da matriz de influência que representa a velocidade normal induzida pelo
painel j no ponto de controle do painel i, quando se tem uma distribuição de fontes com
densidade uniforme e unitária sobre o painel j;
j
σ
é uma distribuição uniforme de fontes sobre o painel j (incógnita do problema);
i
RHSS é um vetor coluna lado direito da equação matricial (“Right Hand Sides”), o qual
possui M elementos, e representa a velocidade normal total induzida no ponto de controle do
painel i devido às contribuições do escoamento incidente e da nuvem de vórtices discretos.
Para um painel genérico i, com ângulo de orientação
i
β , na primeira vez em que este
vetor é calculado tem-se:
iii
cosβvsenβuRHSS
= (4.12)
uma vez que ainda não existem vórtices discretos na região fluida.
Na diagonal principal da matriz de influência
s
ij
K
aparece 0,5K
s
ii
= . Isto significa que a
auto-indução de um painel de fontes está sendo levada em consideração.
4.3.3 – Geração de Vorticidade
Uma das conseqüências mais conhecidas da viscosidade é o desenvolvimento da
camada limite junto das fronteiras sólidas. Outra conseqüência, talvez a mais importante,
porém menos mencionada e analisada, é a geração de vorticidade a partir destas fronteiras. A
fim de melhor compreender o efeito da vorticidade presente no escoamento, analisa-se o
comportamento de três flutuadores localizados em uma região fluida limitada por duas
40
paredes planas e paralelas, conforme ilustra a Figura 4.4. Nesta figura, fica evidente que o
perfil de velocidades do escoamento ainda se encontra em desenvolvimento, uma vez que
existe um núcleo potencial. Este escoamento que se apresenta tem como única finalidade
auxiliar no exemplo didático a seguir.
Figura 4.4 – Geração de vorticidade: um enfoque fenomenológico.
Como bem mostra a Figura 4.4, nas proximidades da parede (interior da camada limite)
o perfil de velocidades do escoamento não é retangular devido à verificação da condição de
aderência e à conseqüente presença da ação da viscosidade. A influência da camada limite
causa um efeito de rotação no sentido anti-horário sobre o flutuador superior e um efeito de
rotação no sentido horário sobre o flutuador inferior; estes efeitos mostram claramente a
presença de vorticidade nas proximidades das paredes. Além disso, nota-se que a partir de
uma certa distância das paredes (fora da camada limite) o perfil de velocidades do escoamento
é retangular e o flutuador central não sofre rotação (apenas translação), mostrando que fora da
camada limite o escoamento é irrotacional.
No presente trabalho, a vorticidade é gerada a cada incremento de tempo
t para anular
o componente tangencial da velocidade; para isso, os vórtices discretos de Lamb são
posicionados de forma a tangenciar o ponto de controle de cada painel plano que representa a
superfície discretizada do corpo (cilindro circular). Para ilustrar o processo de geração de
vorticidade no Método de Vórtices considere a Figura 4.5. Nesta figura está representada de
maneira ilustrativa a superfície discretizada de um cilindro de seção circular em quatro painéis
planos. No algoritmo do Método de Vórtices (veja o Capítulo 5) a vorticidade gerada a partir
da superfície discretizada do corpo sofre um processo de aglutinação instantânea
41
transformando-se num vórtice discreto de Lamb (Alcântara Pereira, 1999). A vorticidade é,
portanto, discretizada e representada por uma nuvem de vórtices discretos de Lamb, os quais
serão transportados pelo escoamento de acordo com os processos de convecção e difusão –
veja a seguir.
Figura 4.5 – Geração de vórtices discretos de Lamb sobre a superfície discretizada de um
corpo.
onde:
4321
co e ,co,co,co são os pontos de controle dos painéis 1, 2, 3 e 4, respectivamente;
eps é a distância de geração dos vórtices, de modo que o núcleo dos vórtices nascentes
tangenciem o ponto de controle de cada painel;
4321
pshed e ,pshed,pshed,pshed são os pontos de desprendimento de vórtices ligados aos
painéis 1, 2, 3 e 4, respectivamente;
0
σ é o raio do núcleo do vórtice de Lamb (Apêndice A).
De maneira similar à montagem da equação matricial de fontes, Equação 4.11, tem-se a
seguinte equação matricial para a geração dos mb vórtices discretos de Lamb ao redor da
superfície do corpo (cilindro circular):
42
=
mb
3
2
1
mb
3
2
1
v
mbmb
v
mb3
v
mb2
v
mb1
v
3mb
v
33
v
32
v
31
v
2mb
v
23
v
22
v
21
v
1mb
v
13
v
12
v
11
RHSV
.
.
.
RHSV
RHSV
RHSV
Γ
.
.
.
Γ
Γ
Γ
K...KKK
.......
.......
.......
K...KKK
K...KKK
K...KKK
(4.13)
onde:
v
ij
K
é um elemento da matriz de influência que representa a velocidade tangencial induzida
pelo vórtice discreto de Lamb posicionado no ponto de desprendimento j sobre o ponto de
controle do painel i, considerando que a intensidade do vórtice localizado em j seja unitária;
j
Γ
é a intensidade do vórtice posicionado no ponto de desprendimento j;
i
RHSV é um vetor coluna lado direito, o qual representa a soma da velocidade do
escoamento incidente e da velocidade induzida por cada um dos vórtices da nuvem sobre o
ponto de controle do painel i, todas na direção tangencial ao painel plano.
A condição de conservação global da circulação é imposta acrescentando-se uma linha e
uma coluna na matriz de influência de vórtices:
() ()
0ΓΓ
N
1k
livres vórtices
k
M
1mgm
nascentes vórtices
m
=+
=+=
(4.14)
Para um painel genérico i, na primeira vez em que o vetor RHSV é calculado tem-se:
iii
senβvcosβuRHSV
= (4.15)
Alternativamente, pode-se representar a equação matricial de fontes, Equação 4.11 por:
[]
{}{}
RHSSSIGMACOUPS = (4.16)
e a equação matricial de vórtices, Equação 4.13, por:
[]
{}{}
RHSVGAMMACOUPV = . (4.17)
De acordo com o algoritmo apresentado no Capítulo 5, sempre que vórtices nascentes forem
posicionados no ponto de desprendimento torna-se necessário o cálculo de uma nova
43
distribuição de fontes sobre os painéis. A geração dos vórtices nascentes mais a geração das
fontes garante a condição de aderência sobre o ponto de controle de cada painel plano.
Assim, deve-se considerar a velocidade que cada vórtice da nuvem induz no ponto de
controle de cada painel. Neste trabalho utiliza-se o modelo do vórtice de Lamb (Panton, 1984)
com esta finalidade (veja o Apêndice A). Sendo m um painel genérico e k um vórtice
arbitrário de intensidade positiva
k
Γ localizado na posição
(
)
kkk
y,xP , os componentes em x
e em y da velocidade induzida neste painel devido à presença do vórtice k são dados
respectivamente por (Lewis, 1991):
[][]
2
km
2
km
kmk
mk
yYMxXM
yYM
2π
Γ
u
+
= (4.18)
[]
[][]
2
km
2
km
kmk
mk
yYMxXM
xXM
2π
Γ
v
+
= (4.19)
onde
()
mm
YM,XM são as coordenadas do ponto de controle do painel genérico m.
As Equações 4.18 e 4.19 são válidas para o modelo do vórtice potencial, ou seja,
quando a distância entre um vórtice arbitrário k e o ponto de controle de um painel for maior
ou igual ao raio do núcleo do vórtice de Lamb (
0km
σr ). Quando a distância entre um
vórtice arbitrário k e o ponto de controle de um painel for menor que o raio do núcleo do
vórtice de Lamb (
0km
σr < ), as Equações 4.18 e 4.19 devem ser substituídas pelas equações
abaixo (Alcântara Pereira, 1999):
[]
[][]
+
=
2
0
2
km
2
km
2
km
kmk
mk
σ
r
5,02572exp1
yYMxXM
yYM
2π
Γ
u (4.20)
[][]
+
=
2
0
2
km
2
km
2
km
kmk
mk
σ
r
5,02572exp1
yYMxXM
xXM
2π
Γ
v
(4.21)
onde
[][]
2
mk
2
mkkm
YMyXMxr += .
Logo, a atualização do vetor
m
RHSS no tempo fica:
[]
+=
=
N
1k
mmkmmkmmm
cosβvsenβuVcosβUsenβRHSS (4.22)
44
Por fim, torna-se possível calcular a velocidade induzida por cada um dos painéis sobre
cada um dos vórtices presentes na nuvem (veja as Equações 4.6 e 4.7).
4.3.4 – Contribuição da Nuvem de Vórtices
Uma das etapas que consome o maior tempo de CPU relaciona-se com a interação
vórtice-vórtice. Quando se utiliza a lei de Biot-Savart, o número de operações realizadas por
um processador é da ordem do quadrado do número total de vórtices discretos presentes no
escoamento.
Quando a distância entre dois vórtices discretos for menor do que o valor do raio do
núcleo do vórtice de Lamb (
0
σ ), a velocidade tangencial induzida neste ponto comporta-se de
acordo com a Figura 4.6b; observe que quando
0r
=
não há auto-indução. Existe na literatura
o modelo do vórtice potencial. Este modelo, embora não possua um núcleo viscoso, é bastante
útil neste trabalho, uma vez que, para distâncias superiores ao valor de
0
σ a velocidade
tangencial que este modelo induz é praticamente igual à velocidade induzida pelo modelo de
Lamb – veja a Figura 4.6a.
Numa nuvem de vórtices discretos, a grande maioria destes vórtices discretos fica a uma
distância superior ao valor do raio do núcleo do vórtice discreto em análise (
0
σ ); deste modo,
é possível utilizar-se do modelo do vórtice potencial, evitando-se o uso do exponencial
presente na expressão do vórtice de Lamb. Este exponencial (veja nas Equações 4.20 e 4.21)
aumenta o tempo de cálculo da velocidade tangencial induzida.
De maneira geral, os componentes nas direções x e y da velocidade total induzida no
vórtice k pelos demais vórtices discretos da nuvem são calculados pelas expressões:
kj,UΓu
N
1j
Vjk
jk,N
=
=
(4.23)
kj,VΓv
N
1j
Vjk
jk,N
=
=
(4.24)
onde
jk,
V
U
e
jk,
V
V
são as velocidades induzidas nas direções x e y, respectivamente, no
vórtice arbitrário k pelo vórtice arbitrário j, se este último possuir intensidade unitária.
45
(a) Vórtice Potencial (b) Vórtice de Lamb
Figura 4.6 – Comportamento da velocidade tangencial induzida.
Para o modelo do vórtice potencial, considerando que a distância do centro do vórtice k
ao centro do vórtice j é maior ou igual ao raio do núcleo do vórtice de Lamb (
0kj
σr ) tem-se:
[
]
[][]
2
jk
2
jk
jk
V
yyxx
yy
2π
1
U
jk,
+
=
(4.25)
[
]
[][]
2
jk
2
jk
jk
V
yyxx
xx
2π
1
V
jk,
+
=
(4.26)
Se
0kj
σr < , tem-se:
[]
[][]
+
=
2
0
2
kj
2
jk
2
jk
jk
V
σ
r
5,02572exp1
yyxx
yy
2π
1
U
jk,
(4.27)
[]
[][]
+
=
2
0
2
kj
2
jk
2
jk
jk
V
σ
r
5,02572exp1
yyxx
xx
2π
1
V
jk,
(4.28)
O raio do núcleo do vórtice de Lamb,
0
σ , é calculado de acordo com os detalhes apresentados
no Apêndice A.
É interessante notar que as expressões 4.23 e 4.24 revelam que um vórtice discreto não
induz velocidade sobre ele mesmo.
46
Um algoritmo eficiente foi desenvolvido para a aceleração da interação vórtice-vórtice.
Detalhes deste algoritmo estão apresentados no Apêndice A4 do trabalho de Alcântara Pereira
(1999), onde é mostrado que o componente x da velocidade induzida em um vórtice discreto k
por um vórtice discreto j (de intensidade unitária), é igual ao componente x da velocidade
induzida no vórtice j pelo vórtice k com sinal contrário, o que revela uma propriedade de anti-
simetria. De forma análoga, o componente y da velocidade induzida em um vórtice k por um
vórtice j (de intensidade unitária), é igual ao componente y da velocidade induzida no vórtice
j pelo vórtice k com sinal contrário. Observando esta particularidade, Mustto
et al. (1998) e
Alcântara Pereira (1999) desenvolveram o algoritmo mencionado anteriormente que calcula
apenas os componentes em x e em y da velocidade induzida no vórtice k pelo vórtice j, ou
seja, o programa faz
jk,kj,
VV
UU = e
jk,kj,
VV
VV
=
.
4.3.5 – Esquemas de Avanço Convectivo
Para se saber a velocidade total induzida no vórtice arbitrário k, basta somar as
contribuições do escoamento incidente, da nuvem de vórtices discretos e das fronteiras
sólidas. Desta forma, obtido o campo de velocidades, a posição de cada vórtice em cada
incremento de tempo na etapa convectiva pode ser calculada por diversos esquemas de
avanço: esquemas de primeira ordem (esquema de Euler), esquemas de segunda ordem
(Adams-Bashforth ou Runge-Kutta), etc. Neste trabalho o avanço convectivo é calculado pelo
esquema de primeira ordem de Euler (Ferziger, 1981). Assim, utilizando este esquema, a nova
posição de um vórtice arbitrário k, após um incremento de tempo
t , nas direções x e y é
dada, respectivamente, por:
()()()
ttutxttx
k
tkk
+=+ (4.29)
()()()
ttvtytty
k
tkk
+=+
(4.30)
onde
k
t
u e
k
t
v são as velocidades totais induzidas no vórtice arbitrário k nas direções x e y,
respectivamente.
47
4.4 – A DIFUSÃO DA VORTICIDADE
A equação que rege o transporte difusivo da vorticidade é dada por:
ω
Re
1
t
ω
2
=
(4.31)
Existem várias técnicas que simulam a difusão da vorticidade, entre elas, o Método de
Vórtices em Células (Christiansen, 1973), Método da Variação do Raio Núcleo do Vórtice
(Leonard, 1980; Kamemoto
et al., 1990), Método da Velocidade de Difusão (Ogami &
Akamatsu, 1991) e Método de Avanço Randômico, inicialmente proposto por Chorin (1973),
baseado na teoria de movimento Browniano desenvolvida por Einstein (1956), e modificado
por Lewis (1991).
Neste trabalho será implementado o Método de Avanço Randômico utilizando a
metodologia de cálculo proposta por Lewis (1991).
Assim, o avanço temporal de um vórtice k da nuvem é dado por:
()()
kk
DIFUSÃOCONVECÇÃO
kk
xxtxttx
+
+=+ (4.32)
()()
kk
DIFUSÃOCONVECÇÃO
kk
yytytty
+
+=+ (4.33)
No Item 4.3.5 foi obtido o avanço convectivo. Agora, será apresentado o movimento
difusivo, caracterizado pelo movimento aleatório das partículas de vorticidade.
O avanço de cada vórtice k na direção radial e no intervalo
()
2π0 é dado
respectivamente por (Alcântara Pereira, 1999):
P
1
ln
Re
t4
r
k
=
(4.34)
2Qπ∆θ
k
= (4.35)
onde P e Q são números randômicos entre 0 e 1. Observe que a presença do número de
Reynolds na Equação 4.31 comprova que os efeitos da viscosidade são levados em
consideração na etapa da difusão da vorticidade.
48
Assim, a difusão de um vórtice discreto k da nuvem após um incremento de tempo
t ,
possui um deslocamento na direção x e um deslocamento na direção y, dados respectivamente
por:
k
DIFUSÃO
rcos∆θx
k
= (4.36)
k
DIFUSÃO
rsen∆θy
k
= (4.37)
4.5 – CARGAS AERODINÂMICAS
Entende-se por cargas aerodinâmicas os esforços mecânicos resultantes da ação que um
fluido em movimento exerce sobre um corpo. Esta ação é exercida pela tensão que atua sobre
a superfície do corpo, a qual pode ser decomposta nas direções normal (pressão) e tangencial.
A integração da tensão atuante sobre a superfície de um corpo, às vezes denominada de
carga aerodinâmica distribuída, dá origem às conhecidas cargas integradas, representadas
pelas forças e momentos aerodinâmicos.
As forças aerodinâmicas podem ser decompostas em dois outros tipos: a força de
sustentação e a força de arrasto. A força de sustentação é resultante da integração dos
componentes normal e tangencial da tensão, e atua em uma direção normal ao escoamento
incidente, dependendo da forma e da orientação do corpo. A força de arrasto, por sua vez, é
resultante da integração dos componentes normal e tangencial da tensão, atuando sempre na
direção do escoamento e opondo-se ao mesmo; esta força depende da forma e orientação do
corpo e das características do escoamento incidente (número de Reynolds).
Em corpos esbeltos que operam com pequeno ângulo de incidência, não ocorre a
separação do escoamento e a esteira que se forma a jusante deste corpo é desprezível, sendo a
força de arrasto originada da integração da tensão tangencial, ou seja, depende principalmente
dos efeitos viscosos. Por outro lado, como foi estudado neste trabalho, em corpos rombudos a
separação do escoamento é evidente e a esteira que se forma é apreciável, sendo dominante o
componente de forma da força de arrasto. Portanto, cabe fixar que neste estudo apenas o
arrasto de forma é calculado.
49
De posse do campo de velocidades do escoamento, procede-se ao cálculo da pressão e
das cargas aerodinâmicas.
Dentre os métodos numéricos para o cálculo da pressão, dois serão destacados neste
texto: o de Lewis (1991), que calcula as cargas aerodinâmicas com base apenas nos vórtices
nascentes em um dado incremento de tempo discreto e o método desenvolvido por Shintani &
Akamatsu (1994), que leva em consideração a influência de todos os vórtices livres presentes
na esteira.
Devido à evidente precisão dos resultados obtidos quando se calcula as cargas por meio
da segunda formulação, tal modelo foi utilizado no cálculo da pressão deste trabalho.
Aplicando-se o operador divergente nas Equações de Navier-Stokes (Equação 3.7), e com
auxílio da Equação da Continuidade (Equação 3.6) chega-se numa equação de Poisson para
pressão. No trabalho de Ricci (2002) pode-se encontrar a dedução detalhada que resultou na
formulação integral abaixo e que permite determinar o valor da pressão em um ponto genérico
do domínio fluido i:
(
)( )
()()
(
)
(
)
()()
()()
()()
+
+
+
=
+
+
+
As
2
i
2
i
ixiy
2
i
2
i
ii
As
2
i
2
i
iyix
i
ωdS
yyxx
yynxxn
π
1
Re
1
ωd
yyxx
yyuxxv
π
1
YdS
yyxx
yynxxn
2π
1
ςY
(4.38)
que discretizada para ser resolvida numericamente, toma a forma:
(
)
(
)
()()
()()
()()
()()
()()
==
=
+
+
+
=
=
+
+
+
mb
ij;1j
jj
2
ij
2
ij
ixjiyj
N
1j
j
2
ij
2
ij
ijjijj
mb
ij;1j
jj
2
ij
2
ij
ijyjijxj
i
γS
yyxx
yynxxn
πRe
1
Γ
yyxx
yyuxxv
π
1
YS
yyxx
yynxxn
2π
1
ςY
(4.39)
onde ς assume o valor 0,5 se o ponto i pertencer a uma fronteira sólida e assume o valor 1,0
se o ponto i não pertencer a uma fronteira sólida, e
j
γ é a densidade de vórtices distribuída
uniformemente sobre o painel j. Neste trabalho a densidade é obtida, num processo inverso,
dividindo-se a intensidade do vórtice de Lamb nascente em cada instante pelo comprimento
do painel:
j
j
j
S
Γ
γ = (4.40)
50
O primeiro somatório do lado direito da Equação 4.39 leva em consideração todos os
vórtices presentes na nuvem e o segundo somatório do lado direito da Equação 4.39 leva em
consideração apenas os vórtices nascentes em cada instante de tempo, que ainda não fazem
parte da nuvem de vórtices.
A Equação 4.39 foi resolvida utilizando-se um Método de Elementos de Contorno,
agrupando o primeiro e o último somatório nas matrizes de influência COUPP (matriz de
pressão) e Ad (matriz lado direito), respectivamente. Assim, a Equação 4.39 pode ser escrita
do seguinte modo:
(
)
(
)
()()
===
+
+
=
mb
1j
jji,
N
1j
j
2
ij
2
ij
ijjijj
mb
1j
jji,
γAdΓ
yyxx
yyuxxv
π
1
YCOUPP
2π
1
(4.41)
O lado direito da equação anterior, por sua vez, pode ser escrito da seguinte forma:
(
)
(
)
()()
==
+
+
=
mb
1j
jji,
N
1j
j
2
ij
2
ij
ijjijj
i
γAdΓ
yyxx
yyuxxv
π
1
RHSP
(4.42)
A aplicação da Equação 4.41 sobre os mb painéis que discretizam a superfície do corpo
conduz ao seguinte sistema linear de equações:
[]
{} { }
RHSPYCOUPP
=
(4.43)
Uma vez conhecidos os valores da incógnita
i
Y para os mb painéis, obtêm-se os valores
do coeficiente de pressão para cada segmento:
1YC
iP
i
+= (4.44)
já que a velocidade sobre a superfície do corpo é nula.
As forças aerodinâmicas são obtidas pela integração da pressão ao longo do corpo. A
força de arrasto atua em cada painel na direção do escoamento incidente, ao passo que a força
de sustentação atua em cada painel na direção normal ao escoamento incidente. Assim, as
referidas forças são dadas por:
()
=
=
mb
1j
j,1jj
nSppD (4.45)
51
()
=
=
mb
1j
j,2jj
nSppL (4.46)
onde,
j
p é a pressão no ponto de controle do painel j;
p é a pressão do escoamento não
perturbado (escoamento incidente);
j,1
n representa o seno do ângulo do painel j;
j,2
n
representa o cosseno do ângulo do painel j. A adimensionalização das equações anteriores
leva à obtenção dos coeficientes de arrasto e sustentação, respectivamente dados por:
()
==
==
mb
1j
j,1jP
mb
1j
j,1jjD
nSCnSpp2C (4.47)
()
==
==
mb
1j
j,2jP
mb
1j
j,2jjL
nSCnSpp2C (4.48)
52
Capítulo 5
IMPLEMENTAÇÃO DO ALGORITMO DO MÉTODO DE
VÓRTICES
5.1 – INTRODUÇÃO
Após a formulação geral do problema proposto partiu-se de uma biblioteca de rotinas
disponível (Moura, 2007) para o desenvolvimento do programa computacional
“MOVINGROUND.FOR” em linguagem FORTRAN. O objetivo deste programa é fornecer
os subsídios para os cálculos de pressão e das forças aerodinâmicas atuantes na superfície do
corpo imerso no domínio fluido
da Figura 3.1.
A seguir, apresenta-se um esquema que mostra a forma como as rotinas são “acionadas”
pelo programa principal (Figura 5.1) e a função de cada uma destas rotinas no código
computacional.
5.2 – DESCRIÇÃO DO FUNCIONAMENTO DAS ROTINAS
Na seqüência apresenta-se uma descrição da função de cada rotina que é acionada pelo
programa principal “MOVINGROUND.FOR”.
53
Figura 5.1 – Estrutura do programa computacional “MOVINGROUND.FOR”.
5.2.1 – Rotina INPUT.FOR
Abre os arquivos de entrada INPUT.DAT e PANELS.DAT e realiza a leitura de todos
os dados da superfície discretizada do corpo e realiza a montagem da superfície plana móvel.
O arquivo PANELS.DAT possui as coordenadas dos pontos extremos de cada painel,
fornecendo a geometria discretizada do corpo. O arquivo INPUT.DAT contém os seguintes
dados de entrada:
stop – número final de incrementos de tempo;
save – freqüência em que se deve “salvar” a posição dos vórtices que formam a nuvem e os
resultados das cargas aerodinâmicas integradas e distribuídas;
start – instante a partir do qual serão calculadas a distribuição média de pressão e a evolução
das cargas aerodinâmicas integradas (força de arrasto e força de sustentação);
lm – comprimento de cada módulo que compõe a superfície plana móvel; cada módulo tem
comprimento igual ao diâmetro do cilindro (veja a Figura 5.2);
54
Figura 5.2 – Representação esquemática do corpo e da superfície plana móvel.
nm –
número total de módulos que compõem a superfície plana móvel;
np – número de painéis utilizados para discretizar a superfície de cada módulo;
mb – número de painéis utilizados para discretizar a superfície do cilindro circular;
eps – distância de geração dos vórtices com relação ao ponto de controle de cada painel;
core – raio do núcleo do vórtice de Lamb, isto é,
0
σ ;
pro – camada protetora: esta camada é utilizada para envolver o corpo de forma a determinar
uma região dentro da qual não é permitida a permanência de vórtices discretos. Este
procedimento deve ser utilizado, uma vez que os vórtices discretos localizados muito
próximos dos painéis de fontes (Método de Painéis) recebem uma indução de velocidade
destes painéis que não refletem a realidade (Ricci, 2002). Esta camada é localizada dentro de
uma região retangular e qualquer vórtice discreto que ultrapassar esta região tem sua posição
investigada com a finalidade de averiguar se o mesmo se encontra no interior da camada
protetora ou não; em caso positivo, o vórtice é deslocado para fora da camada protetora, sendo
devolvido ao domínio fluido.
gap – porcentagem de deslocamento do ponto de controle em relação à superfície
discretizada. Devido ao efeito de curvatura da superfície de um corpo, pode-se deslocar o
ponto de controle de um painel em direção à superfície real do corpo (Ricci, 2002).
55
vel – velocidade do escoamento incidente;
alpha – ângulo do escoamento incidente;
delt – valor do incremento de tempo;
re – número de Reynolds;
h – distância entre a superfície plana móvel e o centro do cilindro circular (veja a Figura 5.2).
5.2.2 – Rotina DATAPRGP.FOR
Recebe o valor dos pontos extremos dos painéis utilizados na discretização da superfície
plana móvel. Em seguida a rotina define uma região retangular e calcula a camada protetora
ao redor da superfície; mais tarde esta camada será utilizada na investigação da posição dos
vórtices discretos e na conseqüente reflexão ou não destes a partir da superfície plana móvel.
Por fim, a rotina calcula o valor do ponto de controle, do ângulo de orientação, do
comprimento e do ponto de desprendimento do vórtice discreto referente a cada painel. Esses
valores são impressos no arquivo SURFACE.DAT. A Figura 5.3 apresenta, como exemplo,
um cilindro de seção circular discretizado em 10 painéis planos, assim como a superfície
plana, a qual está discretizada em 20 painéis planos (5 módulos); na figura, identificam-se os
pontos extremos dos painéis (pontos pretos), os pontos de controle dos painéis (pontos azuis)
e os pontos de desprendimento dos vórtices discretos de Lamb (pontos vermelhos). Veja que
como não há a geração de vorticidade a partir da superfície plana móvel, não há representação
de pontos de desprendimento de vórtices discretos (pontos vermelhos) sobre esta superfície.
5.2.3 – Rotina DATAPRSB.FOR
Recebe o valor dos pontos extremos dos painéis utilizados na discretização da superfície
do cilindro circular. Em seguida, a rotina define uma região retangular e calcula a camada
protetora ao redor do corpo; mais tarde esta camada será utilizada na investigação da posição
dos vórtices discretos e na conseqüente reflexão ou não destes a partir da superfície
discretizada do cilindro circular.
Na seqüência, a rotina calcula o valor do ponto de controle, do ângulo de orientação, do
comprimento e o valor apropriado das coordenadas da posição de geração do vórtice referente
a cada painel; cada vórtice discreto tem seu centro posicionado sobre uma normal que passa
56
pelo ponto de controle do painel e vale
0
σeps
=
. Esses valores são impressos no arquivo
BODY.DAT (Figura 5.3).
De posse das coordenadas dos pontos de controle de cada painel, a rotina realiza um
deslocamento destes pontos no sentido de aproximá-los da superfície real do cilindro, de
maneira a melhorar os resultados da simulação numérica.
Figura 5.3 – Representação dos pontos extremos dos painéis, dos pontos de controle dos
painéis e dos pontos de desprendimento de vórtices discretos.
5.2.4 – Rotina COUPS.FOR
Esta rotina calcula os coeficientes da matriz de influência,
s
ij
K , da Equação 4.11. Para
isso, ela segue os seguintes passos:
- calcula os componentes da velocidade nas direções x e y que dependem da distribuição de
fontes
()
xσ
no sistema de coordenadas local; esta distribuição de fontes, no momento, é igual
à unidade;
- passa esses componentes do sistema de coordenadas local no painel para o sistema de
coordenadas global (fixo ao corpo);
- decompõe esses componentes na direção normal, impondo a condição de contorno de
Neumann (Equação 3.8).
Vale lembrar que os coeficientes da matriz de influência de fontes não sofrem variação
ao longo da simulação, porque dependem apenas da geometria do problema e, no caso em
estudo, tanto o cilindro circular como a superfície plana têm suas posições fixas durante a
57
simulação numérica. O movimento da superfície plana é simulado deixando-se de gerar
vórtices discretos sobre ela.
5.2.5 – Rotina COUPV.FOR
Essa rotina calcula os coeficientes da matriz de influência,
v
ij
K , da Equação 4.13. Para
isso, ela calcula a velocidade tangencial que um vórtice de intensidade unitária localizado no
painel j induz no ponto de controle de um painel i.
Vale lembrar que, tal como ocorre na rotina “acionada” anteriormente, os coeficientes
da matriz de influência de vórtices não sofrem variação ao longo da simulação, uma vez que
dependem apenas da geometria do problema e, no caso em estudo, tanto o cilindro circular
como a superfície plana móvel têm suas posições invariáveis ao longo do tempo; o
movimento da superfície plana é simulado deixando-se de gerar vórtices discretos sobre esta.
5.2.6 – Rotina RHSS.FOR
Calcula o vetor coluna lado direito de fontes da Equação 4.11.
No instante
0t
=
da simulação numérica o escoamento é potencial, uma vez que ainda
não há vórtices discretos no domínio fluido. Nessa situação, o vetor coluna lado direito de
fontes é calculado conforme a Equação 5.1:
mmm
cosβvsenβuRHSS
= (5.1)
Nos instantes seguintes haverá vórtices discretos no escoamento formando a esteira
viscosa. Deste modo, a cada instante deve-se proceder à atualização do vetor RHSS de forma
a garantir que a distribuição de fontes esteja sendo capaz de anular a velocidade normal
(condição de impenetrabilidade). A referida atualização do vetor RHSS é feita segundo a
Equação 5.2:
()
=
+=
N
1k
mmkmmkmmm
cosβvsenβucosβvsenβuRHSS (5.2)
onde N leva em consideração todos os vórtices discretos, incluindo-se os vórtices nascentes.
58
5.2.7 – Rotina RHSV.FOR
Calcula o vetor coluna lado direito de vórtices da Equação 4.13.
No instante
0t
=
da simulação numérica o escoamento é potencial, uma vez que ainda
não há vórtices no domínio fluido. Nessa situação, o vetor coluna lado direito de vórtices
sofre apenas a influência do escoamento incidente:
mmm
senβvcosβuRHSV
= (5.3)
Nos instantes seguintes, esta rotina é acionada após a convecção e difusão dos vórtices,
ou seja, deve-se atualizar o vetor coluna lado direito de vórtices para garantir a condição de
escorregamento nulo, através da seguinte equação:
()
=
+=
N
1k
mmkmmkmmm
senβvcosβusenβvcosβuRHSV (5.4)
Diferentemente da rotina anterior, a rotina RHSV.FOR só leva em conta a influência
dos vórtices sobre os painéis que discretizam o cilindro circular, uma vez que a condição de
escorregamento nulo deve ser satisfeita apenas no corpo e não na superfície plana móvel.
5.2.8 – Rotina GAUSSPIV.FOR
As equações matriciais de fontes, de vórtices e de pressão são resolvidas através do
Método de Eliminação de Gauss com Condensação Pivotal Parcial.
5.2.9 – Rotina MODCOUP.FOR
Acrescenta uma linha e uma coluna na matriz de influência de vórtices (Equação 4.13)
para impor a condição de conservação global da circulação. Esta condição obedece a seguinte
equação:
() ()
0ΓΓ
N
1k
livres vórtices
k
M
1mgm
nascentes vórtices
m
=+
=+=
(5.5)
onde mg (
npn
m
mg
= ) é o número de painéis que discretizam a superfície plana móvel; para
m variando de 1 até mg, tem-se os painéis que discretizam a superfície plana móvel.
59
O acréscimo de apenas uma linha e uma coluna na matriz de influência de vórtices
deve-se ao fato de neste trabalho usar-se o artifício de gerar vórtices discretos apenas na
superfície discretizada do cilindro circular, de forma a representar a situação de movimento
relativo, conforme dito anteriormente. Caso fossem gerados vórtices discretos também na
superfície plana, dever-se-ia acrescentar duas linhas e duas colunas na matriz de influência de
vórtices.
Por fim, essa rotina calcula os coeficientes da matriz de pressão (COUPP), Equação
4.43.
5.2.10 – Rotina GENERAT.FOR
Realiza em cada incremento de tempo, a geração de mb novos vórtices discretos com
intensidades dadas pela solução do sistema de equações 4.13. Ressalta-se que são gerados
apenas mb vórtices discretos por instante de tempo devido ao fato da situação de movimento
relativo ser simulada através da não geração de vórtices discretos sobre a superfície plana.
5.2.11 – Rotina COMPUMVM.FOR
Realiza o cálculo dos componentes nas direções x e y da velocidade induzida nos
pontos de controle de cada painel plano por cada um dos vórtices discretos de Lamb presentes
na nuvem. De acordo com o processo iterativo apresentado na Figura 5.1, esta rotina é
primeiro acionada quando houver geração de vórtices novos. Esta primeira chamada da rotina
é necessária para a atualização do vetor coluna lado direito de fontes. A rotina é novamente
acionada porque os vórtices presentes na nuvem deslocaram-se por convecção e por difusão.
Esta segunda chamada da rotina é necessária para a atualização do vetor coluna lado direito de
vórtices.
5.2.12 – Rotina COMPUCVC.FOR
Rotina que calcula o campo de velocidades levando-se em consideração a contribuição
do escoamento incidente, das fronteiras sólidas (cilindro circular e solo) e da nuvem de
vórtices discretos sobre cada vórtice discreto que compõe a nuvem.
60
5.2.13 – Rotina PRESSURE.FOR
Calcula a distribuição de pressão instantânea sobre a superfície discretizada do corpo
(exatamente sobre o ponto de controle de cada painel plano) e os coeficientes instantâneos de
arrasto e de sustentação. Esta rotina começa a ser acionada quando o número de incrementos
de tempo for igual ao valor de start.
Esses cálculos são feitos de acordo com a formulação proposta por Shintani &
Akamatsu (1994).
5.2.14 – Rotina DIFFUS.FOR
Realiza o processo de difusão da vorticidade pelo Método de Avanço Randômico, onde
são gerados para cada vórtice discreto, em cada incremento de tempo, dois números
randômicos P e Q entre 0 e 1; veja a Seção 4.4.
5.2.15 – Rotina CONVEC.FOR
A convecção da nuvem de vórtices é realizada por esta rotina através do esquema de
avanço de primeira ordem de Euler; veja o Item 4.3.5.
5.2.16 – Rotina REFLECT.FOR
Esta rotina aciona a rotina REFLEC.FOR para que os vórtices que migrarem para o
interior do corpo (e os vórtices que se encontrarem no interior da camada protetora) sejam
devolvidos para o domínio fluido.
No arquivo de saída chamado CONSERV.DAT armazena-se em cada instante de
tempo, o valor da Equação 5.5, ou seja, verifica-se se a conservação da circulação está sendo
satisfeita.
5.2.17 – Rotina REFLECS.FOR
Tem o mesmo papel da rotina anterior, promovendo a reflexão dos vórtices que
migrarem para baixo da superfície plana móvel ou que estiverem no interior da camada
protetora ao redor dela.
61
5.2.18 – Rotina CONSERV.FOR
Para realizar a conservação da circulação do escoamento, a rotina CONSERV.FOR faz
1mbmb += e atribuí o valor igual a zero para RHSV(mb) no vetor coluna lado direito de
vórtices.
5.2.19 – Rotina MODM.FOR
Para a resolução da equação matricial de vórtices, com as novas condições comentadas
anteriormente, a rotina MODM.FOR é acionada para realizar a correção do número mb de
painéis do corpo, que foi alterado pela rotina CONSERV.FOR. Esta rotina também imprime
no arquivo CONTROL.DAT os valores do instante de tempo (step) e os valores extras obtidos
na resolução da equação matricial de vórtices, devido à imposição da conservação da
circulação (veja a Equação 5.5). Como a solução do sistema linear de equações (veja na
Equação 4.13) é obtida pelo Método de Eliminação de Gauss com Condensação Pivotal
Parcial, utilizou-se como estratégia para a verificação da Equação 5.5 a inclusão de uma linha
e uma coluna adicionais na Equação 4.13. Esta estratégia resulta na obtenção de um vetor
incógnita de vórtices com um valor a mais; este valor (
ε ) deve tender a zero. No arquivo
CONTROL.DAT, portanto, encontra-se o valor de
ε para cada instante de tempo da
simulação numérica.
5.2.20 – Rotina PRINT.FOR
Imprime a posição, a intensidade, a velocidade induzida e o valor do raio do núcleo para
cada vórtice discreto que compõe a nuvem no arquivo de saída do tipo WAKE XXXX.DAT.
Os valores instantâneos das forças de arrasto e de sustentação são impressos no arquivo
FORBOD.DAT e os valores instantâneos da distribuição de pressão ao longo da superfície
discretizada do corpo, no arquivo CPBOD XXXX.DAT. O símbolo XXXX se refere a um
incremento de tempo arbitrário da simulação; por exemplo, em
10t
=
e com 0,05delt
=
,
200XXXX = , logo, têm-se WAKE0200.DAT e CPBOD0200.DAT. O arquivo
FORBOD.DAT é sempre atualizado.
5.2.21 – Rotina AVERAGE.FOR
62
Calcula os valores médios para o coeficiente de pressão
m
P
C ao longo da superfície
discretizada do corpo. Os valores médios são calculados entre
deltstartt
=
e deltstopt
= e
são impressos no arquivo de saída CPBODY.DAT.
Os valores médios dos coeficientes de arrasto e de sustentação são impressos no arquivo
AVERAGE.DAT.
Para as conclusões mais importantes deste trabalho, que serão apresentadas no Capítulo
7, necessitou-se de análises minuciosas de pós-processamento. Para cumprir esta etapa, são
apresentados e analisados gráficos da evolução ao longo do tempo das cargas aerodinâmicas
integradas, gráficos das distribuições de pressão média e instantânea e visualizações do
comportamento do campo de velocidades nas proximidades do corpo. Estes tratamentos
gráficos utilizaram os “softwares” GRAPHER 4 e TECPLOT 360.
63
Capítulo 6
ANÁLISE DE RESULTADOS
6.1 - INTRODUÇÃO
Este capítulo apresenta os principais resultados obtidos para a simulação numérica do
escoamento bidimensional, incompressível e em regime não-permanente de um fluido
newtoniano com propriedades constantes que incide sobre um cilindro circular. As condições
geométricas impostas ao cilindro circular o definem como isolado e estacionado nas
proximidades de uma superfície plana móvel. O algoritmo do Método de Vórtices Discretos
descrito no Capítulo 5 e que foi implementado em linguagem FORTRAN resultou no
programa “MOVINGROUND.FOR”. Na Figura 6.1 identificam-se as superfícies do cilindro
circular e da parede plana móvel. Estas superfícies foram discretizadas e representadas por
painéis planos sobre os quais se distribuíram singularidades do tipo fontes com densidade
constante, impondo-se, assim, a condição de impenetrabilidade sobre o ponto de controle de
cada painel. Vórtices discretos de Lamb foram gerados nas vizinhanças da superfície
discretizada do cilindro circular para satisfazer a condição de escorregamento nulo. Como o
propósito deste trabalho é o de eliminar a influência da camada limite criada junto à superfície
plana móvel, a condição de escorregamento nulo somente foi imposta sobre o ponto de
controle de cada painel montado para representar a superfície do corpo. Ainda na Figura 6.1
pode-se perceber a estrutura da esteira de vórtices discretos de Lamb que se forma a jusante
do corpo.
64
O primeiro passo constituiu-se na aferição do novo código computacional
implementado e, para isso, foram realizadas várias simulações numéricas para o caso do
cilindro circular isolado, ou seja, corpo posicionado a uma distância da superfície plana móvel
onde a influência do efeito solo não mais se faz presente; neste caso:
1000dh = .
Após a aferição inicial do código, procedeu-se à análise do escoamento ao redor do
cilindro circular na presença da parede plana se movendo em relação ao cilindro com a
mesma velocidade do escoamento incidente. Os resultados a serem apresentados se referem
ao comportamento das cargas aerodinâmicas integradas (coeficientes de arrasto e de
sustentação) e distribuídas (coeficiente de pressão) e do número de Strouhal. Também se
utiliza recursos gráficos de visualização da formação da esteira a jusante do corpo.
6.2 – PARÂMETROS UTILIZADOS NA SIMULAÇÃO
NUMÉRICA
Nas simulações numéricas foram consideradas duas classes de parâmetros: aqueles
relacionados ao método numérico e aqueles afetos ao fenômeno físico. A seguir serão
esclarecidas as influências que tais parâmetros têm sobre os resultados da simulação
numérica.
6.2.1 – Parâmetros Relacionados com o Método Numérico
a) Número de painéis planos (M)
1
S
1
S
2
S
3
+θ
Figura 6.1 – Cilindro circular localizado próximo a uma superfície plana móvel.
Superfície Plana Móvel
65
O número de painéis planos utilizados para discretizar e representar as duas fronteiras
sólidas do problema tem influência direta na precisão dos resultados da simulação numérica;
sendo assim, aumentando-se o número de painéis utilizados para discretizar as fronteiras
sólidas, aumenta-se a precisão dos resultados obtidos. Uma das grandes dificuldades impostas
pelo Método de Vórtices consiste no tempo de CPU; sabe-se que quanto maior o número de
painéis utilizados, maior será o número total de vórtices discretos presentes na nuvem. O
tempo de CPU, quando se utiliza a lei de Biot-Savart, torna-se proporcional ao quadrado do
número de vórtices discretos presentes na nuvem para o cálculo da interação vórtice-vórtice.
Neste trabalho, devido às limitações computacionais, a superfície do cilindro circular foi
discretizada em 300 painéis planos de comprimentos iguais (
300mb
=
), ao passo que a
superfície plana móvel foi representada por 10 módulos (
10nm
=
) de comprimento igual ao
diâmetro do cilindro circular, sendo cada módulo discretizado em 30 painéis (
30np
=
),
totalizando 300 painéis; veja a Figura 5.2. A estratégia numérica de se gerar vórtices apenas a
partir da superfície do cilindro permitiu que se concentrasse todo o esforço computacional
sobre os vórtices nascentes no corpo.
b) Incremento de tempo (t )
O incremento de tempo foi escolhido a partir de vários testes, tendo como primeira
aproximação a expressão utilizada no trabalho de Mustto et al. (1998), ou seja:
mb
2kπ
t = ,
1k0
<
(6.1)
onde mb é o número de painéis planos utilizados para discretizar a superfície do corpo.
Na escolha do incremento de tempo foram avaliados esquemas de avanço convectivo de
primeira ordem (esquema de Euler) e de segunda ordem (esquema de Adams-Bashforth),
verificando que o valor de 0,05 para o incremento de tempo era satisfatório quando se usava o
esquema de primeira ordem de Euler. A utilização do esquema de segunda ordem de Adams-
Bashforth (Ferziger, 1981) exigiria o conhecimento, também, do campo de velocidades no
instante imediatamente anterior.
c) Posição de desprendimento dos vórtices discretos e raio do núcleo do
vórtice de Lamb (eps e
0
σ ).
66
Em todas as simulações numéricas realizadas neste trabalho os vórtices nascentes, em
cada incremento de tempo, foram desprendidos a partir de uma distância dos pontos de
controle dos painéis igual ao raio do núcleo do vórtice de Lamb ( 0,001epsσ
0
== ). Esta
técnica de geração de vórtices discretos para simular a vorticidade impõe que o núcleo de
cada vórtice nascente sempre tangencie o ponto de controle de cada painel onde ele é gerado
(Alcântara Pereira, 1999). O valor nominal calculado para o raio do núcleo do vórtice de
Lamb é determinado através da Equação A.17 (ver Apêndice A). Yang & Huang (1999)
apresentaram uma técnica para um tratamento mais preciso da geração dos vórtices nascentes,
onde se impõe uma variação do raio do núcleo do vórtice dependente do comportamento da
camada limite local que se desenvolve. Mustto (2004) em sua Tese de Doutorado apresentou
uma argumentação semelhante à anteriormente comentada, onde a escolha do parâmetro eps
dependia, também, da espessura da camada limite.
Deve-se ressaltar que ao se utilizar o Método de Painéis existe um erro quando se
calcula a velocidade induzida por um painel (com singularidade distribuída) sobre um vórtice
discreto localizado muito próximo deste painel. Ricci (2002) mostrou que ao se utilizarem
painéis de fontes com densidade constante, inscritos na superfície real do corpo, os erros se
conservam da ordem de 1% quando o ponto de controle é deslocado para 23% da distância
entre o ponto de controle do painel e a superfície real do corpo ( 0,23gap
=
). Como no chão a
superfície discretizada coincide com a superfície real, não se deve lançar mão de tal artifício.
Além do deslocamento acima citado, o presente trabalho também reposiciona alguns vórtices
discretos caso estes se encontrem muito próximos aos painéis. Esta proximidade é
estabelecida através de uma camada protetora que envolve todos os painéis. Esta camada
protetora corresponde a uma distância do vórtice até o painel igual ou inferior a 40% do
comprimento do painel (
0,40pro =
). Os vórtices que se encontrarem no interior desta camada
devem ser reposicionados para fora desta.
6.2.2 – Parâmetros Relacionados com o Fenômeno Físico
a) Número de Reynolds (Re)
O valor adotado para o número de Reynolds foi o de
5
10 . Este valor foi escolhido
porque permite que os resultados das simulações numéricas possam ser comparados com os
resultados experimentais disponíveis na literatura. Para a situação do cilindro isolado
comparam-se os resultados numéricos da presente simulação com aqueles de Blevins (1984);
67
para a situação do cilindro na presença de uma parede plana que se move com a mesma
velocidade do fluido, comparam-se os resultados numéricos da presente simulação com
aqueles de Nishino (2007).
b) Distância vertical entre o centro do cilindro circular e a superfície plana
móvel (h/d)
Este parâmetro é fundamental para indicar os valores de dh que submetem o corpo à
presença do efeito de bloqueio da parede plana móvel.
c) Velocidade da superfície plana móvel
Os resultados experimentais encontrados na literatura se referem ao escoamento sobre
um cilindro circular na presença de uma parede plana que se move com a mesma velocidade
do escoamento incidente. Nenhum outro resultado experimental diferente desta situação foi
encontrado, mas sabe-se que se houver movimento relativo entre a superfície plana e o
escoamento incidente, haverá formação de camada limite a partir do solo.
d) Parâmetros relacionados com a modelagem de turbulência
Os escoamentos turbulentos caracterizam-se pelo alto grau de complexidade à medida
que o número de Reynolds aumenta. Desta maneira, as dificuldades para a solução das
Equações 3.6 e 3.7 tornam-se muito grandes. Surge, portanto, a necessidade de se obter
soluções aproximadas através de diversas técnicas estabelecidas de modelagem da
turbulência. Entre estas técnicas citam-se: (a) SND – Simulação Numérica Direta (“Direct
Numerical Simulation” - DNS); (b) SMR – Simulação via Equações Médias de Reynolds
(“Reynolds-Averaged Navier-Stokes Equations” - RANS); (c) SGE – Simulação de Grandes
Escalas (“Large Eddy Simulation” – LES).
Assim, o número de Reynolds investigado neste trabalho requer que se leve em conta a
transferência de energia que ocorre entre as grandes escalas do escoamento e as escalas sub-
malha. Se fosse utilizado neste trabalho, por exemplo, LES, haveria a referida transferência de
energia; lembrar que para a utilização de um modelo de turbulência devem-se filtrar as
equações governantes e modelar as escalas que não são resolvidas. Matematicamente, os
termos difusivos da equação da quantidade de movimento são afetados pela presença de uma
viscosidade turbulenta, que tem como finalidade realizar a referida transferência de energia.
Um modelo de turbulência se presta para o cálculo desta viscosidade turbulenta. O modelo da
Função Estrutura de Velocidade de Segunda Ordem (Alcântara Pereira
et al., 2002) foi
68
adaptado para o Método de Vórtices tendo como grande vantagem aquela de se trabalhar com
o conceito de diferenças de velocidades ao invés do conceito de derivadas. Deve-se, portanto,
conhecer os parâmetros afetos a esta modelagem para corretamente simular novas situações
de escoamentos. Este modelo, embora não utilizado neste trabalho, será incorporado nos
próximos desenvolvimentos a serem realizados.
6.3 – ESCOAMENTO AO REDOR DE UM CILINDRO
CIRCULAR ISOLADO
No caso deste trabalho, um cilindro circular sujeito a um escoamento incidente é dito
isolado quando se encontra a uma distância suficientemente grande da superfície plana móvel,
de modo que os mecanismos do efeito solo não sejam mais sentidos por este corpo. Sendo
assim, quando a superfície plana móvel não exercer mais influência sobre as cargas
aerodinâmicas que atuam sobre o corpo considera-se este último como isolado. No código
computacional adota-se
dh =1000 para caracterizar o cilindro circular como isolado.
Desta forma, o programa “MOVINGROUND.FOR” foi utilizado para a análise do
comportamento aerodinâmico de um cilindro circular isolado (
dh =1000) e, principalmente,
quando estacionado nas proximidades de uma superfície plana móvel (
dh variável e
representando efeito solo). Os resultados numéricos obtidos para
dh =1000 foram
comparados com os resultados experimentais de Blevins (1984) e os resultados numéricos
obtidos para h/d variável e representando efeito solo foram comparados com os resultados
experimentais de Nishino (2007).
A Tabela 6.1 apresenta uma comparação entre os resultados numéricos e experimentais
disponíveis na literatura e os resultados obtidos pela presente simulação. As aferições
realizadas permitiram simular o escoamento ao redor do cilindro circular isolado com 1000
avanços no tempo (
1000stop =
), sendo que os valores médios dos coeficientes de arrasto e
sustentação, bem como do número de Strouhal, foram calculados entre 28,3t
inicial
=
( 566step = da simulação e 169.800 vórtices discretos presentes na nuvem) e 48,0t
final
=
(
960step =
da simulação e 288.000 vórtices presentes na nuvem).
69
Tabela 6.1 – Valores médios dos coeficientes de arrasto e de sustentação e número de
Strouhal para um cilindro circular isolado.
5
10Re =
D
C
L
C
St
Blevins (1984) 1,20 - 0,19
Mustto et al. (1998) 1,22 - 0,22
Alcântara Pereira et al. (2002) 1,21 0,04 0,22
Presente Simulação 1,25 0,02 0,21
Analisando a Tabela 6.1, deve-se comentar que Mustto
et al. (1998) representaram a
superfície real do cilindro circular utilizando o Teorema do Círculo (Milne-Thompson, 1955),
o qual tem a grande vantagem de garantir a condição de impenetrabilidade em todos os pontos
da superfície real do corpo. Neste caso, a condição de escorregamento nulo deve ser satisfeita
apenas sobre um determinado número de pontos de geração distribuídos ao redor da superfície
do cilindro. Alcântara Pereira
et al. (2002) utilizaram painéis planos sobre os quais se
distribuíram vórtices com densidade constante para representar a superfície discretizada do
cilindro circular (Katz & Plotkin, 1991); além disso, conforme já mencionado, incorporaram
um modelo de turbulência na formulação matemática do Método de Vórtices.
Comparando-se os resultados numéricos da presente simulação com os valores
experimentais de Blevins (1984), os quais possuem uma incerteza de
10%±
, observa-se uma
boa concordância entre eles. Verifica-se que o resultado obtido para o coeficiente de arrasto é
um pouco superior ao experimental, o que é uma característica das simulações numéricas
bidimensionais. Portanto, o coeficiente de arrasto se mostrou bem mais sensível aos efeitos
tridimensionais do que o número de Strouhal. Deve ser mencionado que na formulação
matemática do problema (Capítulo 3) foi desprezado o termo responsável pela deformação
dos tubos de vorticidade,
()
uω , o qual deve ser implementado nas simulações numéricas
tridimensionais.
A Figura 6.2 apresenta o comportamento da distribuição média de pressão sobre a
superfície discretizada do cilindro circular. É feita uma comparação entre os resultados da
teoria potencial, os resultados experimentais de Blevins (1984) e os resultados numéricos de
Mustto (1998).
Verifica-se que a separação do escoamento para a presente simulação ocorre quando
° 68θ
, o que é um resultado aceitável. A literatura apresenta valores experimentais distintos
70
para este ângulo:
°
82θ (Milne-Thompson, 1955) e °
78θ (Son & Hanratty, 1969), ambos
para
5
101,0Re = . O ângulo θ foi definido na Figura 6.1.
Os resultados de Mustto
et al. (1998) para o coeficiente de pressão diferem bastante do
resultado experimental de Blevins (1984) devido à metodologia empregada no cálculo da
pressão. Mustto
et al. (1998) utilizaram a formulação proposta por Lewis (1991), a qual leva
em consideração apenas os vórtices nascentes em um dado incremento de tempo para o
cálculo da pressão. O presente trabalho apresenta um resultado mais preciso, uma vez que faz
uso da formulação proposta por Shintani & Akamatsu (1994), que leva em consideração não
só os vórtices nascentes em um dado incremento de tempo, mas também aqueles presentes na
esteira viscosa.
Figura 6.2 – Distribuição média de pressão ao longo da superfície discretizada do cilindro
circular isolado (Euler,
300mb
=
, 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ).
Nota-se, ainda, que há uma diferença entre o resultado da presente simulação e o
resultado experimental na região de descolamento da camada limite. Esta discrepância pode
ser minimizada mediante a utilização de painéis menores nesta região. Uma outra fonte de
erro está associada à posição de geração dos vórtices. Conforme dito anteriormente, os raios
θ
71
dos núcleos dos vórtices possuem o mesmo valor; além disso, todos os vórtices são
igualmente posicionados em relação ao ponto de controle. Sendo assim, o presente trabalho
viola a física real do problema, na medida em que impõe o mesmo efeito viscoso ao longo de
toda a superfície do corpo, quando na verdade sabe-se que existe uma camada limite laminar
que evolui para uma camada limite turbulenta (Yang & Huang, 1999). Além disso, a
convecção da vorticidade pode ser mais bem simulada utilizando-se incrementos de tempo
menores associados ao uso do Método de Expansão em Multipólos (Greengard & Rohklin,
1987). Com o objetivo de melhorar ainda mais os resultados, a etapa difusiva pode ser
realizada através do Método do Crescimento do Raio do Núcleo do Vórtice Modificado
(Rossi, 1996, 1997, 2005, 2006). Este método se encontra em fase inicial de implementação
numérica no Grupo de Método de Vórtices da UNIFEI e irá garantir a repetibilidade dos
resultados obtidos, característica não apresentada pelo Método de Avanço Randômico.
A Figura 6.3 mostra a evolução das cargas aerodinâmicas integradas ao longo do tempo.
Observa-se que após a passagem do transiente numérico ( 20,0t
) o coeficiente de
sustentação oscila em torno do valor nulo (uma vez que o cilindro circular é um corpo
rombudo e simétrico) e o coeficiente de arrasto oscila em torno de 1,25. O número de
Strouhal, que mede a freqüência com que os pares contra-rotativos de vórtices são
desprendidos a partir dos pontos de separação, apresenta o valor médio de 0,21. Observa-se
ainda, que o período de oscilação do coeficiente de arrasto é duas vezes maior do que o
período de oscilação do coeficiente de sustentação, o que é uma característica intrínseca do
cilindro circular isolado. Este comportamento ocorre porque o coeficiente de arrasto oscila
uma vez para cada vórtice que se desprende ou no lado superior ou no lado inferior da
superfície do cilindro.
Ainda na Figura 6.3, são identificados cinco pontos importantes que serão utilizados
para entender os mecanismos de desprendimento das estruturas vorticosas ao longo de um
período da curva que representa o coeficiente de sustentação.
Uma outra observação importante está relacionada à amplitude de oscilação do
coeficiente de sustentação na Figura 6.3. Verifica-se que a amplitude de oscilação média
positiva é de aproximadamente 1,01, ao passo que a amplitude de oscilação média negativa é
de aproximadamente 1,11, ou seja, o valor líquido é de aproximadamente zero, comprovando
que não há sustentação no cilindro circular isolado. O valor do coeficiente de sustentação
obtido não é nulo devido a erros numéricos.
72
Figura 6.3 – Evolução das cargas aerodinâmicas integradas ao longo do tempo para o cilindro
circular isolado (Euler,
300mb
=
, 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ).
A Figura 6.4 apresenta a distribuição de pressão instantânea que atua sobre a superfície
discretizada do cilindro circular em cada um dos cinco instantes marcados na Figura 6.3.
Figura 6.4 – Distribuição instantânea de pressão ao longo da superfície discretizada do
cilindro circular isolado (Euler,
300mb
=
,
0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ).
1,01
1,11
cálculo da média
θ
73
Nota-se que no instante representado pelo ponto A existe uma zona de baixa pressão
compreendida entre 67º e 175°, medida a partir do bordo de ataque do cilindro circular no
sentido horário. Esta zona de baixa pressão corresponde a uma estrutura vorticosa horária
sendo desprendida na parte superior do cilindro (veja a Figura 6.5), e implica em uma força de
sustentação positiva (veja o ponto A da Figura 6.3). Observa-se ainda, na Figura 6.3, que o
coeficiente de arrasto está aumentando quando a estrutura vorticosa está sendo desprendida.
Isto acontece porque o vórtice é desprendido e o coeficiente de arrasto aumenta até atingir um
valor máximo, quando a estrutura vorticosa começará a ser incorporada pela esteira viscosa
(ponto B).
Na Figura 6.5(a) se encontram distribuídos no domínio fluido 172.500 vórtices discretos
de Lamb.
No instante representado pelo ponto B da Figura 6.4 há uma região de baixa pressão
aproximadamente constante entre 67° e 288°, a qual corresponde à passagem do coeficiente
de sustentação de positivo para negativo. Neste instante, o vórtice desprendido anteriormente
começa a ser incorporado pela esteira formada a jusante do corpo (Figura 6.6) e o coeficiente
de arrasto diminui.
(a) Ponto A: posição dos vórtices na esteira.
(b) Ponto A: campo de velocidades nas vizinhanças do corpo.
Figura 6.5 – Detalhes do desprendimento de estruturas vorticosas em um cilindro circular
isolado no instante
39,4t =
(Euler, 300mb
=
,
0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ).
74
Na Figura 6.6(a) se encontram distribuídos no domínio fluido 179.400 vórtices discretos
de Lamb. Identifica-se, também, a formação da esteira de von Kármán de uma maneira muito
bem definida. Os três primeiros pares contra-rotativos de estruturas vorticosas unidos por
folhas de vorticidade e visualizados a partir do cilindro na Figura 6.6 (a) já definem a correta
formação da esteira pulsante.
Observa-se que no instante representado pelo ponto C da Figura 6.4 existe uma zona de
baixa pressão compreendida entre 188° e 279°. Esta zona de baixa pressão corresponde a uma
estrutura vorticosa anti-horária desprendida na parte inferior do cilindro (Figura 6.7), e
implica em uma sustentação negativa (veja o ponto C da Figura 6.3).
(a) Ponto B: Posição dos vórtices na esteira.
(b) Ponto B: Campo de velocidades nas vizinhanças do corpo.
Figura 6.6 - Detalhes do desprendimento de estruturas vorticosas em um cilindro circular
isolado no instante
6,40t = (Euler, 300mb
=
, 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ).
Na Figura 6.7(a) se encontram distribuídos no domínio fluido 185.700 vórtices discretos
de Lamb.
Verifica-se ainda que os vórtices desprendidos nos instantes representados pelos pontos
A e C possuem um movimento de rotação com sentidos opostos.
75
(a) Ponto C: Posição dos vórtices na esteira.
(b) Ponto C: Campo de velocidades nas vizinhanças do corpo.
Figura 6.7 - Detalhes do desprendimento de estruturas vorticosas em um cilindro circular
isolado no instante 6,41t = (Euler,
300mb
=
, 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re =
).
Tal como ocorre no instante representado pelo ponto B, o ponto D da Figura 6.4
apresenta uma região de baixa pressão aproximadamente constante entre 67° e 292°, a qual
corresponde à passagem do coeficiente de sustentação de negativo para positivo. Neste
instante, o vórtice desprendido anteriormente começa a ser incorporado pela esteira formada a
jusante do corpo; veja a Figura 6.8.
(a) Ponto D: Posição dos vórtices na esteira.
(b) Ponto D: Campo de velocidades nas vizinhanças do corpo.
Figura 6.8 - Detalhes do desprendimento de estruturas vorticosas em um cilindro circular
isolado no instante
9,42t = (Euler, 300mb
=
, 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ).
76
O período se completa no instante representado pelo ponto E, o qual apresenta uma
região de baixa pressão compreendida entre 50° e 158°. Esta região de baixa pressão
corresponde a uma estrutura desprendida na parte superior do cilindro (Figura 6.9), e implica
em uma força de sustentação positiva (veja o ponto E da Figura 6.3).
Os mecanismos acima descritos se repetem, ou seja, pares contra-rotativos de vórtices
são gerados alternadamente a partir dos pontos de separação do cilindro de seção circular, o
que faz com que a esteira viscosa formada a jusante do corpo tenha um caráter oscilatório,
caracterizando a formação da esteira de von Kármán apresentada na Figura 6.10.
Neste momento é importante destacar o elevado tempo dispendido nas simulações
numéricas. Para se obter a esteira com 300.000 vórtices mostrada na Figura 6.10, levou-se
132 horas e 17 minutos de tempo de CPU em um processador PENTIUM IV de 1700 MHz.
(a) Ponto E: Posição dos vórtices na esteira.
(b) Ponto E: Campo de velocidades nas vizinhanças do corpo.
Figura 6.9 - Detalhes do desprendimento de estruturas vorticosas em um cilindro circular
isolado no instante
0,44t = (Euler, 300mb
=
, 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ).
Considera-se, portanto, que o código computacional está aferido e apto a simular o
escoamento ao redor do cilindro circular estacionado nas proximidades da superfície plana
móvel, com a utilização dos mesmos parâmetros usados no caso do cilindro isolado.
77
Figura 6.10 – Posição dos vórtices na esteira ao final da simulação numérica do cilindro
circular isolado (Euler,
300mb
=
, 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ).
6.4 – ESCOAMENTO AO REDOR DE UM CILINDRO
CIRCULAR SUBMETIDO AO EFEITO SOLO
O esclarecimento da física presente no escoamento ao redor de um corpo estacionado
nas proximidades de uma superfície plana parada é bastante difícil de ser obtido, uma vez que
há um forte entrelaçamento entre a esteira do corpo e a camada limite formada junto à
superfície plana fixa.
A maioria dos casos práticos de engenharia ocorre quando o corpo se movimenta em
relação à superfície. Nishino (2007) observou que quando o fundo móvel de um túnel de
vento se desloca com a mesma velocidade do escoamento incidente, praticamente não há
formação de camada limite sobre a esteira rolante, o que ocorre na prática com os veículos.
Verificando que existiam poucos estudos envolvendo corpos rombudos nesta situação,
Nishino (2007) realizou testes em túneis de vento para elucidar o escoamento ao redor de um
cilindro circular estacionado próximo a uma esteira rolante que se movia com a mesma
velocidade do escoamento incidente, eliminando a dificuldade causada pela camada limite
que se formaria junto à esteira rolante, caso ela estivesse parada.
Assim, o presente trabalho se compromete em simular numericamente principalmente
esta situação, de forma a demonstrar as potencialidades do Método de Vórtices Discretos e
sua capacidade de reproduzir resultados experimentais. Todo o esforço computacional ficará
disponível para se gerar mais vórtices sobre o cilindro circular, já que na superfície plana não
serão gerados vórtices para que não haja formação de camada limite sobre ela, de forma a
simular a situação em que um corpo se move relativamente a uma superfície com a mesma
velocidade do escoamento incidente.
78
A análise do problema se inicia com a observação da Figura 6.11. Nela estão
representados os resultados experimentais obtidos por Roshko
et al. (1975) e por Nishino
(2007), bem como os resultados numéricos obtidos neste trabalho.
Os resultados de Roshko
et al. (1975) correspondem ao caso em que o cilindro circular
está estacionado nas proximidades de uma superfície plana fixa. Observa-se que o coeficiente
de arrasto sofre uma queda brusca à medida que o corpo se aproxima da superfície a uma
distância menor do que
1,00dh = . Associa-se esta queda do arrasto à injeção da vorticidade
da pista na esteira do corpo; veja mais detalhes nas análises de Ricci (2002) e de Moura
(2007).
Figura 6.11 – Comparação entre os resultados numéricos e experimentais obtidos para o
coeficiente de arrasto com relação aos efeitos tridimensionais.
Entre os testes realizados em túnel de vento por Nishino (2007), a situação de maior
interesse para efeito de comparação com os resultados deste trabalho é aquela referente ao
comportamento das cargas aerodinâmicas quando os efeitos tridimensionais são parcialmente
retirados do problema. Para isso, Nishino (2007) utilizou placas (“end-plates”) nas
extremidades do cilindro para inibir a manifestação de efeitos de ponta, como mostra a Figura
6.12, retirada de sua Tese de Doutorado.
79
No entanto, na situação de movimento relativo, Nishino (2007) verificou que, para o
caso em que as placas não eram utilizadas (caso essencialmente tridimensional), à medida que
o corpo se aproximava do solo o coeficiente de arrasto aumentava suavemente, contrariando o
que se pensava até então.
Figura 6.12 – Aparato experimental usado nos testes em túnel de vento por Nishino (2007).
Como o presente trabalho estuda o fenômeno do efeito solo através de simulações
numéricas bidimensionais, os resultados obtidos neste trabalho serão comparados apenas com
testes feitos por Nishino (2007) para os casos em que foram usadas as placas. Nesta condição,
nota-se que o código bidimensional representa satisfatoriamente a física do problema, uma
vez que os resultados numéricos se aproximam bastante dos resultados experimentais
correspondentes ao caso em que os efeitos tridimensionais são minimizados.
Da mesma maneira como ocorreu com o coeficiente de arrasto, a Figura 6.13 mostra
que o código também fornece bons resultados acerca do coeficiente de sustentação. Verifica-
se que à medida que o corpo se aproxima do solo o cilindro passa a ter uma sustentação
positiva, devido à ocorrência de uma região de baixa pressão na parte superior do corpo.
Observa-se, ainda, que diferentemente do coeficiente de arrasto, o coeficiente de
sustentação não apresenta grandes diferenças entre os casos em que as placas são usadas e
aqueles em que elas são retiradas; nota-se ainda que o comportamento do coeficiente de
sustentação para um cilindro circular estacionado nas vizinhanças de uma superfície plana
móvel é bastante semelhante ao comportamento deste coeficiente quando o cilindro se
encontra próximo a uma superfície plana fixa.
80
Figura 6.13 – Comparação entre resultados numéricos e experimentais obtidos para o
coeficiente de sustentação envolvendo o cilindro circular submetido ao efeito solo.
A Tabela 6.2 apresenta o valor médio do coeficiente de arrasto e o número de Strouhal
para o cilindro circular situado a uma distância 0,95dh
=
da superfície plana móvel,
situação esta em que o efeito solo se faz presente. Resultados para outras relações
dh são
fornecidos no Apêndice B.
Tabela 6.2 – Valor médio do coeficiente de arrasto e número de Strouhal para o cilindro
circular submetido ao efeito solo para
0,95dh
=
.
5
10Re =
D
C
St
Nishino (2007) – com placas 1,31 -
Presente Simulação 1,40 0,21
Nota-se que os resultados da simulação numérica estão bem próximos dos resultados
experimentais. É preciso dizer também, que em um estudo experimental a força de arrasto
medida pelos instrumentos equivale ao arrasto total, ao passo que o código utilizado neste
trabalho calcula apenas o componente de forma da força de arrasto (aquele resultante da
81
integração da pressão e da decomposição desta na direção x). Portanto, este fato já introduz
um pequeno erro nos resultados do presente estudo.
A Figura 6.14 (a) mostra a evolução no tempo das cargas aerodinâmicas integradas para
o cilindro circular localizado a uma distância
0,95dh
=
da superfície plana móvel. Observa-
se que a amplitude de oscilação do coeficiente de sustentação é maior para o caso do cilindro
sujeito ao efeito solo do que para o cilindro isolado (Figura 6.14 (b)). Além disso, nota-se que
na situação de efeito solo o período de oscilação do coeficiente de arrasto não possui
amplitude aproximadamente constante tal como ocorre para o cilindro isolado; verifica-se que
o período do
D
C possui amplitudes intercaladas: ora maiores, ora menores. Ainda nesta
figura, são identificados cinco pontos importantes que serão utilizados para comparar o
comportamento aerodinâmico do cilindro circular sujeito ao efeito solo, com o
comportamento aerodinâmico deste mesmo corpo, porém, isolado. O Apêndice B apresenta a
evolução das cargas aerodinâmicas para outras relações
dh .
Ao contrário do que ocorreu com o cilindro circular isolado, a situação de efeito solo
provocou uma sustentação positiva no cilindro circular. Vê-se que a amplitude de oscilação
média positiva é de aproximadamente 1,56, ao passo que a amplitude de oscilação média
negativa é de aproximadamente 1,33, ou seja, o valor líquido é de 0,23, provendo o cilindro
de uma força de sustentação dirigida para cima.
(a)
95,0dh = (b) 1000dh =
Figura 6.14 – Comparação da evolução das cargas aerodinâmicas para os casos do cilindro
circular (a) submetido ao efeito solo e (b) isolado (Euler,
300mb
=
, 0,05t = ,
0,001epsσ
0
=
=
,
5
10Re = ).
1,56
1,33
1,01
1,11
cálculo da média
c
cálculo da média
82
A Figura 6.15 apresenta a distribuição de pressão que atua sobre a superfície
discretizada do cilindro circular em cada um dos cinco instantes marcados na Figura 6.14(a).
Figura 6.15 – Distribuição instantânea de pressão ao longo da superfície discretizada do
cilindro circular submetido ao efeito solo (Euler,
300mb
=
, 0,05t
=
, 0,001epsσ
0
== ,
5
10Re = , 0,95dh
=
).
Nota-se que no instante representado pelo ponto A existe uma zona de baixa pressão
compreendida entre 67° e 171°, medida a partir do bordo de ataque do cilindro circular no
sentido horário. Esta zona de baixa pressão corresponde a uma estrutura vorticosa sendo
desprendida no sentido horário na parte superior do cilindro (Figura 6.16) e implica em uma
força de sustentação positiva (veja o ponto A da Figura 6.14(a)). Este ponto de sustentação
máxima está associado ao valor máximo do coeficiente de arrasto correspondente à amplitude
menor.
Na Figura 6.16(a) se encontram distribuídos no domínio fluido 214.500 vórtices
discretos de Lamb.
Compare as Figuras 6.10 e 6.16(a) e observe a formação das esteiras de von Kármán.
θ
83
(a) Ponto A: Posição dos vórtices na esteira.
(b) Ponto A: Campo de velocidades nas vizinhanças do corpo.
Figura 6.16 - Detalhes do desprendimento de estruturas vorticosas em um cilindro circular
submetido ao efeito solo no instante 4,46t
=
(Euler, 300mb
=
, 0,05t = ,
0,001epsσ
0
=
= ,
5
10Re =
, 0,95dh
=
).
A formação da esteira de von Kármán na Figura 6.16(a) se apresenta com uma
influência bastante significativa do efeito da parede plana móvel (efeito de bloqueio).
No instante representado pelo ponto B da Figura 6.15 há uma região de baixa pressão
aproximadamente constante entre 58° e 292°, a qual corresponde à passagem do coeficiente
de sustentação de positivo para negativo. Neste instante, a estrutura vorticosa desprendida no
ponto A começa a ser incorporada pela esteira formada a jusante do corpo (Figura 6.17). É
interessante notar ainda, que este momento está associado a um valor mínimo do coeficiente
de arrasto correspondente à amplitude menor.
Aliás, observa-se que o coeficiente de arrasto aumenta até atingir um valor máximo nos
instantes em que estruturas vorticosas são desprendidas (pontos A, C e E) e que depois sofre
uma queda à medida que as estruturas que foram desprendidas vão sendo incorporadas pela
esteira viscosa que se forma a jusante do corpo.
84
(a) Ponto B: Posição dos vórtices na esteira.
(b) Ponto B: Campo de velocidades nas vizinhanças do corpo.
Figura 6.17 - Detalhes do desprendimento de estruturas vorticosas em um cilindro circular
submetido ao efeito solo no instante 6,47t
=
(Euler, 300mb
=
, 0,05t = ,
0,001epsσ
0
=
= ,
5
10Re =
, 0,95dh
=
).
Analisando-se o instante representado pelo ponto C, identifica-se uma zona de baixa
pressão compreendida entre 175° e 275º. Esta zona de baixa pressão corresponde a uma
estrutura vorticosa desprendida no sentido anti-horário na parte inferior do cilindro (Figura
6.18) e implica em uma força de sustentação negativa (veja o ponto C da Figura 6.14(a)). Este
ponto de sustentação mínima está associado ao valor máximo do coeficiente de arrasto
correspondente à amplitude maior. Vê-se que as pressões nos pontos de descolamento para o
cilindro sob efeito solo são menores do que as pressões nos pontos de descolamento do
cilindro isolado (veja as Figuras 6.4 e 6.15). Entretanto, comparando apenas os pontos de
descolamento superior e inferior do cilindro sujeito ao efeito solo, verifica-se que o ponto de
descolamento inferior (ponto C) possui uma pressão bastante baixa, cuja conseqüência é a
criação de uma zona onde não se verifica a presença de vorticidade entre os vórtices superior
e inferior (Figura 6.18).
No instante representado pelo ponto D da Figura 6.15 há uma região de baixa pressão
aproximadamente constante entre 67° e 283°, a qual corresponde à passagem do coeficiente
de sustentação de negativo para positivo. Neste instante, a estrutura desprendida no ponto C
começa a ser incorporada pela esteira que se forma a jusante do corpo (Figura 6.19). É
interessante notar ainda, que este momento está associado a um valor mínimo do coeficiente
de arrasto correspondente à amplitude maior.
85
(a) Ponto C: Posição dos vórtices na esteira.
(b) Ponto C: Campo de velocidades nas vizinhanças do corpo.
Figura 6.18 - Detalhes do desprendimento de estruturas vorticosas em um cilindro circular
submetido ao efeito solo no instante 7,48t
=
(Euler,
300mb
=
, 0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re = ,
0,95dh
=
).
Em função da baixa pressão ocorrida no momento em que a estrutura vorticosa do ponto
C foi desprendida, verifica-se na Figura 6.19(b) a presença bem mais nítida da zona sem
vorticidade na parte traseira do cilindro circular.
(a) Ponto D: Posição dos vórtices na esteira.
(b) Ponto D: Campo de velocidades nas vizinhanças do corpo.
Figura 6.19 - Detalhes do desprendimento de estruturas vorticosas em um cilindro circular
submetido ao efeito solo no instante 1,50t
=
(Euler, 300mb
=
, 0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re = , 0,95dh
=
).
zona sem vorticidade
86
O período do coeficiente de sustentação se completa no instante representado pelo ponto
E, o qual apresenta uma região de baixa pressão compreendida entre 67° e 175°. Esta região
de baixa pressão, tal como ocorreu no instante representado pelo ponto A, corresponde a um
vórtice desprendido no sentido horário na parte superior do cilindro (Figura 6.20), e implica
em uma força de sustentação positiva (veja o ponto E da Figura 6.14(a)). Este ponto de
sustentação máxima também está associado ao valor máximo do coeficiente de arrasto
correspondente à amplitude menor.
A Figura 6.21 mostra a esteira formada ao final da simulação do escoamento ao redor
do cilindro circular localizado a uma distância
0,95dh
=
da superfície plana móvel. Como
não foram gerados vórtices discretos sobre a superfície plana, o tempo de simulação numérica
foi o mesmo do caso do cilindro isolado, ou seja, foram necessárias 132 horas e 17 minutos de
tempo de CPU utilizando um processador PENTIUM IV de 1700 MHz para se chegar a uma
esteira de 300.000 vórtices. Esteiras formadas ao final de simulações numéricas de outras
relações
dh são fornecidas no Apêndice B.
(a) Ponto E: Posição dos vórtices na esteira.
(b) Ponto E: Campo de velocidades nas vizinhanças do corpo.
Figura 6.20 - Detalhes do desprendimento de estruturas vorticosas em um cilindro circular
submetido ao efeito solo no instante 51,2t
=
(Euler,
300mb
=
, 0,05t
=
, 0,001epsσ
0
=
= ,
5
10Re = ,
0,95dh
=
).
87
É importante mencionar o fato de que o comprimento da superfície plana móvel nas
simulações numéricas é finito. Verifique na Figura 6.21 que a superfície plana está
representada por dez módulos com comprimentos unitários e a esteira viscosa se estende até
aproximadamente sessenta módulos a partir do cilindro circular. A superfície plana não foi
simulada com um comprimento maior devido a limitações computacionais, já que tal tarefa
exigiria o uso de um número muito grande de painéis planos na discretização da referida
superfície. A condição de escorregamento nulo não foi imposta sobre a parede plana móvel,
porém a condição de impenetrabilidade, que foi imposta, se verifica até o décimo módulo. No
algoritmo desenvolvido, todos os vórtices ao longo do eixo dos x que migrarem para o interior
da parede plana móvel são reposicionados acima da camada protetora.
Figura 6.21 – Posição dos vórtices na esteira ao final da simulação numérica do cilindro
circular submetido ao efeito solo (Euler,
300mb
=
, 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re = ,
0,95dh
=
).
88
Capítulo 7
CONCLUSÕES E SUGESTÕES
7.1 – CONCLUSÕES
A influência da camada limite gerada a partir da superfície plana foi desconsiderada
neste trabalho, uma vez que o solo se move com a mesma velocidade do escoamento
incidente (Nishino, 2007). Esta supressão da camada limite foi de suma importância para que
se concentrassem todos os esforços computacionais sobre os vórtices desprendidos a partir da
superfície discretizada do cilindro circular. Esta estratégia permitiu que se obtivesse ao final
da simulação numérica 300.000 vórtices discretos de Lamb. Este número foi compatível com
a estrutura computacional disponível para a realização dos experimentos numéricos. Cada
simulação numérica levou 132 horas e 17 minutos num processador PENTIUM IV de 1700
MHz. Houve, contudo, uma melhora significativa dos resultados; situações anteriores para a
análise do efeito solo dentro do mesmo grupo de pesquisa onde este trabalho foi desenvolvido
obtiveram esteiras finais com 100.000 vórtices discretos (Ricci (2002), Silva de Oliveira
et al.
(2005) e Moura (2007)).
No entanto, a aferição de outros parâmetros como, por exemplo, incremento de tempo,
posição de geração dos vórtices discretos e raio do núcleo do vórtice de Lamb necessitam de
maiores investigações.
A acurácia e os tempos de CPU atuais devem ser melhorados. Deseja-se num futuro
muito próximo repetir estas análises utilizando-se o Método de Expansão em Multipólos
89
(Greengard & Rohklin, 1987) associado ao Método do Crescimento do Raio do Núcleo do
Vórtice Modificado (Rossi, 2006) chegando a um número final de vórtices da ordem de
6
10 .
Estas duas mudanças significativas exigem a utilização de computação de alto desempenho. A
estas duas mudanças deseja-se incluir, também, modelagem sub-malha de turbulência
(Alcântara Pereira
et al., 2002).
Se os propósitos anteriores forem atingidos, a extensão do código para o espaço
tridimensional será a próxima meta a ser cumprida.
Os resultados numéricos obtidos neste trabalho se apresentam muito promissores. As
análises das Figuras 6.11 e 6.13 mostram que para uma situação bidimensional o código
computacional responde bem, conforme era de se esperar. Para o caso do cilindro circular
isolado, observou-se que o coeficiente de arrasto foi pouco sensível aos efeitos
tridimensionais. Entretanto, quando o cilindro circular está sob efeito solo, os efeitos
tridimensionais influenciaram muito mais no valor deste coeficiente.
Na Figura 6.11, pode-se notar que o resultado da simulação numérica bidimensional
difere em muito no que diz respeito ao comportamento do coeficiente de arrasto obtido por
Nishino (2007), sem a presença das placas que inibem os efeitos tridimensionais. Os
resultados experimentais de Roshko
et al. (1975), situação para a superfície plana parada,
foram colocados na mesma figura para efeito de comparação.
Observando a Figura 6.14, verifica-se que tanto para o caso do cilindro submetido ao
efeito solo como para o caso do cilindro isolado, quando o coeficiente de arrasto completa um
período, uma estrutura vorticosa foi desprendida e, ou ela é incorporada à esteira logo em
seguida, ou a estrutura formada imediatamente antes dela é que é incorporada. Esta
verificação explica o fato dos casos do cilindro sob efeito solo e isolado apresentarem o
mesmo valor para o número de Strouhal, ou seja, as duas situações têm pares contra-rotativos
de vórtices sendo desprendidos com a mesma freqüência. É interessante notar que é possível
observar este fato mediante a análise do comportamento do coeficiente de arrasto quando a
sustentação assume os valores máximo e mínimo: para o caso do cilindro circular isolado, um
valor máximo ou mínimo do coeficiente de sustentação está sempre associado a uma ascensão
do coeficiente de arrasto; de forma parecida, para o caso do cilindro circular sujeito ao efeito
solo, um valor máximo ou mínimo do coeficiente de sustentação está sempre associado a um
ponto de máximo do coeficiente de arrasto.
90
Nesta linha de raciocínio, é interessante comentar que no caso do cilindro sob efeito
solo, toda vez que uma estrutura vorticosa está sendo incorporada à esteira, verifica-se que
neste momento o coeficiente de arrasto passa por um ponto de mínimo.
As características do escoamento ao redor de um cilindro circular estacionado nas
proximidades de uma superfície plana móvel variam em função da distância do cilindro até o
solo. Entretanto, Nishino (2007) observou que existem três faixas da relação
dh onde
fenômenos importantes acontecem. Para distâncias maiores da superfície em relação ao
cilindro (
1,0dh > ), grandes pares contra-rotativos de vórtices são gerados na parte traseira
do cilindro circular, o que produz um alto valor para o coeficiente de arrasto (
1,30C
D
); isto
pôde ser comprovado através da simulação numérica deste trabalho. Numa faixa intermediária
(
1,0dh0.85 << ), o desprendimento dos pares contra-rotativos de vórtices passa a ser
intermitente; como o caso mais bidimensional testado por Nishino (2007) em túnel de vento
não analisou esta faixa, não há como comparar os resultados experimentais com os obtidos
através da simulação numérica. Por fim, a pequenas distâncias do solo (
0,85dh < ) Nishino
(2007) observou que o desprendimento de vórtices é interrompido, criando-se uma zona sem a
presença da vorticidade.
Ainda, a não inclusão dos efeitos tridimensionais no problema analisado mostrou que as
simulações numéricas produzem valores para o número de Strouhal que são insensíveis à
presença de tais efeitos. Porém, o coeficiente de arrasto apresenta-se com valores maiores que
os experimentais apresentados por Nishino (2007) quando este autor faz uso das placas
inibidoras dos efeitos tridimensionais. Nota-se, também, que as maiores discordâncias nos
resultados numéricos obtidos para o coeficiente de pressão estão concentradas na região de
descolamento da camada limite. Uma possível causa para este comportamento inadequado
pode estar relacionada com a discretização da superfície do cilindro circular, pois os painéis
planos têm o mesmo comprimento nesta região. Portanto, um refinamento na distribuição de
painéis nesta região associado a uma adequada modelagem de turbulência podem minimizar
esta discrepância. Uma outra alternativa para melhorar os resultados na região de
descolamento da camada limite consiste em se realizar um tratamento mais adequado para os
vórtices nascentes (Yang & Huang, 1999). No presente trabalho os vórtices nascentes são
posicionados com os seus núcleos distantes de um mesmo valor do ponto de controle do
painel de onde se desprendem.
Finalmente, considera-se o código computacional capaz de reproduzir a física do efeito
solo (parede plana móvel) sobre um corpo para as condições de escoamento bidimensional,
91
incompressível e em regime não-permanente de um fluido Newtoniano com propriedades
constantes.
7.2 – SUGESTÕES
Como sugestões para a continuação do desenvolvimento do presente algoritmo,
primeiramente, seria interessante testar o comportamento do programa para perfis
aerodinâmicos (corpos esbeltos) sem a presença da camada limite gerada a partir do solo
(Ahmed, 2005) e (Ahmed & Sharma, 2005).
O número de Reynolds investigado neste trabalho para o caso do cilindro circular está
na faixa sub-crítica. Outras faixas de número de Reynolds podem ser investigadas tanto para
corpos rombudos quanto para corpos esbeltos.
Uma extensão necessária do presente código consiste na inclusão dos efeitos de
interação térmica a partir de um Método de Partículas de Calor (Ogami, 2001) e (Alcântara
Pereira & Hirata, 2003).
Alcântara Pereira & Hirata (2003) incluíram na biblioteca de rotinas utilizada neste
trabalho efeitos de interação térmica. A metodologia considerou apenas efeitos de convecção
forçada. Foi mostrado que existe uma similaridade entre a Equação do Transporte da
Vorticidade:
()
ω
Re
1
ω
t
ω
2
=+
u (7.1)
e a Equação da Energia:
()
T
ReP
r
1
T
t
T
2
=+
u (7.2)
Em analogia ao Algoritmo de Separação da Parte Viscosa da Equação do Transporte da
Vorticidade (Chorin, 1973) pode-se separar a Equação 7.2 nos efeitos de convecção e de
difusão de calor. A nuvem de vórtices discretos deve ser utilizada para realizar a convecção
do calor.
92
Ogami (2001) propôs dois métodos para simular a convecção natural. Nesta abordagem
o campo de vorticidades é modificado pela presença do campo de temperaturas. A criação de
vorticidade a partir do calor pode ser analisada aumentando-se a intensidade dos vórtices
discretos ou criando-se um par de vórtices a partir de uma partícula de calor.
A inclusão de efeitos harmônicos de oscilação no corpo pode ser feita segundo a
abordagem de Moura (2007) ou aquela de Recicar (2007).
Para a inclusão da abordagem de Moura (2007), efeitos de oscilação do corpo de
pequena amplitude, o presente código pode ser facilmente adaptado.
No entanto, para a inclusão da abordagem de Recicar (2007) no presente trabalho,
efeitos de oscilação do corpo de amplitude qualquer, o presente código necessita de mudanças
radicais. Estas mudanças se devem ao movimento relativo de oscilação que o corpo
apresentará em relação à superfície plana. Torna-se necessário o cálculo de todas as matrizes
de influência do presente código como função do tempo.
93
REFERÊNCIAS BIBLIOGRÁFICAS
AHMED, M. R. (2005), “Aerodynamics of a Cambered Airfoil in Ground Effect”,
International Journal of Mechanics Research, v 32, n 2, pp 157-183.
AHMED, M. R., SHARMA, S. D. (2005), “An Investigation on the Aerodynamics of a
Symmetrical Airfoil in Ground Effect”, Experimental Thermal and Fluid Science, v 29,
pp 633-647.
ALCÂNTARA PEREIRA, L. A. (1999), Simulação Numérica do Escoamento em torno de
um Corpo de Forma Arbitrária Utilizando o Método de Vórtices Discretos, Dissertação
de Mestrado, EFEI/IEM/DME.
ALCÂNTARA PEREIRA, L. A., RICCI, J. E. R., HIRATA, M. H., SILVEIRA NETO,
A. (2002),
“Simulation of the Vortex-Shedding Flow about a Circular Cylinder with
Turbulence Modeling”,
CFD Journal, v 11, n 3, October, pp 315-322.
ALCÂNTARA PEREIRA, L. A., HIRATA, M. H. (2003), “Heat Transfer in the Wake
Behind a Body Using a Particle Method”,
17
th
International Congress of Mechanical
Engineering, Proceedings of COBEM 2003, November 10-14, São Paulo, SP, Brazil.
ANDERSON, C. R., GREENGARD, C. (1991), “Vortex Dynamics and Vortex Methods”,
AMS Lectures in Appl. Math., v 28.
BARBA, L. A., LEONARD, A., ALLEN, C. B. (2004), “Vortex Method with Fully Mesh-
Less Implementation for High-Reynolds Number Flow Computations”, European
Congress on Computational Methods in Applied Sciences and Engineering, ECCOMAS
2004, Jyväskylä.
94
BATCHELOR, G. K. (1967), An Introduction to Fluid Dynamics, Cambridge University
Press.
BEALE, J. T., MAJDA, A. (1981), “Rates of Convergence for Viscous Splitting of the
Navier-Stokes Equations”, Math. Comp., 37:243-259.
BEALE, J. T., MAJDA, A. (1982a), “Vortex Methods I: Convergence in Three
Dimensions”, Math. Comp., 39:1-27.
BEALE, J. T., MAJDA, A. (1982b), “Vortex Methods II: High-Order Accuracy in Two and
Three Dimensions”,
Math. Comp., 39:29-52.
BEALE, J. T., MAJDA, A. (1985), “High-Order Accurate Vortex Methods with Explicit
Velocity Kernels”, J. Comp. Phys., 58:188-208.
BEARMAN, P. W., ZDRAVKOVICH, M. M. (1978), “Flow around a Circular Cylinder
Near a Plane Boundary, Journal of Fluid Mechanics, v 89, pp 33-47.
BIRKHOFF, G., FISHER, J. (1959), “Do Vortex Sheets Roll UP?”, Rend. Circ. Mat.
Palermo Ser. 2, 8:77-90.
BIRKHOFF, G. (1962), “Helmholtz and Taylor Instability”, Proc. Symp. Appl. Math., Amer.
Math. Soc.. v XIII, pp 55-76.
BLEVINS, R. D. (1984), Applied Fluid Dynamics Handbook, Van Nostrand Reinhold, Co.
BRECHT, S. H., FERRANTE, J. R. (1990), “Vortex-In-Cell Calculations in Three
Dimensions”,
Comp. Phys. Comm., 58:25-54.
BURESTI, G., LANCIOTTI, A. (1979), “Vortex Shedding from Smooth and Roughened
Cylinders in Cross-Flow Near a Plane Surface”, Aeronautical Quarterly, v 30, pp 305-
321.
CHACALTANA, J. T. A., BODSTEIN, G. C. R., HIRATA, M. H. (1994), “2-D
Interaction of a Point-Vortex with a Thin Airfoil Near a Ground Plane”, Anais do V
ENCIT, São Paulo, Brazil.
95
CHACALTANA, J. T. A., BODSTEIN, G. C. R., HIRATA, M. H. (1995), “Analytical
Study of the Time-Dependent 2D Interaction of a Thin Airfoil and a Vortex in the
Presence of a Ground Plane”,
XIII Brazilian Congress of Mechanical Engineering
XIII COBEM, Belo Horizonte, MG, Brazil.
CHORIN, A. J. (1973), “Numerical Study of Slightly Viscous Flow”, Journal of Fluid
Mechanics, v 57, pp 785-796.
CHRISTIANSEN, J. P. (1973), “Numerical Simulation of Hydrodynamics by the Method of
Point Vortices”,
J. Comp. Phys., v 13, pp 363-379.
COTTET, G. H. (1987), “Convergence of a Vortex-In-Cell Method for the Two-
Dimensional Euler Equations”, Math. Comp., 49(180):407-425.
COTTET, G. H., KOUMOUTSAKOS, P. (2000), Vortex Methods. Theory and Practice,
Cambridge University Press.
COTTET, G. H., PONCET, P. (2004), “Advances in Direct Numerical Simulations of 3D
Wall-Bounded Flows by Vortex-In-Cell Methods”, J. Comp. Phys., 193(1):136-158.
COUET, B., BUNEMAN, O., LEONARD, A. (1981), “Simulation of Three-Dimensional
Incompressible flows with a Vortex-In-Cell Method”, J. Comp. Phys., 39(2):305-328.
DEGOND, P., MAS-GALLIC, S. (1989), “The Weighted Particle Method for Convection
Diffusion Equations Part 1: The Case of an Isotropic Viscosity”, Math Comput, 53:485-
507.
EBIANA, A. B., BARTHOLOMEW, R.W. (1996), “Design Considerations for Numerical
Filters Used in Vortex-In-Cell Algorithms”, Comp. & Fluids, 25(1):61-75.
EINSTEIN, A. (1956), Investigation on the Theory of Brownian Motion, Dover, NY.
FERZIGER, J. H. (1981), Numerical Methods for Engineering Application, John Wiley &
Sons Inc..
96
FISHELOV, D. (1990), “A New Vortex Scheme for Viscous Flows”, J. Comput. Phys,
86:211-224.
GREENGARD, C. (1985), “The Core Spreading Vortex Method Aproximates the Wrong
Equation”, Journal of Computational Physics, v 61, pp 345-348.
GREENGARD, L., ROKHLIN, V. (1987), “A Fast Algorithm for particle simulations”, J.
Comp. Phys., 73:325-348.
HALD, O., MAUCERI DEL PRETE, V. (1978), “Convergence of Vortex Methods for
Euler’s Equations”,
Math. Comp., 32:791-801.
HALD, O. (1979), “Convergence of Vortex Methods for Euler’s Equations II”, SIAM J. Num.
Anal., 16:726-755.
HE, F., SU, T. C. (1999), “An Improved Offset Model for Vorticity Shedding from a Solid
Boundary in Discrete Vortex Element Method”, Department of Mechanical
Engineering, FAU: Florida Atlantic University, Boca Raton, FL 33431, pp 1-24.
HESS, J. I., SMITH, A. M. O. (1967), “Calculation of Potential Flow about Arbitrary
Bodies”, Progress in Aeronautical Sciences, v 8, pp 1-138.
HIWADA, M., MABUCHI, I., KUMADA, M., IWAKOSHI, H. (1986), “Effect of the
Turbulent Boundary Layer Thickness on the Flow Characteristics around a Circular
Cylinder Near a Plane Surface”, (in Japanese),
Transactions of the Japan Society of
Mechanical Engineers B
, v 52, pp 2566-2574.
KAMEMOTO, K., NAKAHARA, N., KAWATA, Y., IMAMURA, K., KANEKO, T.
(1990),
“Numerical Simulation of Vortex Flows Interacting with Vibrations under Flow
Gates”, IAHR Symposium, Belgrade, Yugoslave, N4, pp 1-8.
KAMEMOTO, K. (1993), “Procedure to Estimate Unstead Pressure Distribution for Vortex
Method” (In Japanese), Trans. Jpn. Soc. Mech. Eng., v 59, n 568B, pp 3708-3713.
97
KAMEMOTO, K., ZHU, B., OJIMA, A. (2000), “Attractive Features of an Advanced
Vortex Method and its Subjects as a Tool of Lagrangian LES”, 14
th
Japan Society of
CFD Symposium, Tokyo, December, pp 1-10.
KAMEMOTO, K. (2004), “On Contribution of Advanced Vortex Element Methods Toward
Virtual Reality of Unsteady Vortical Flows in the New Generation of CFD”,
Proceedings of the 10
th
Brazilian Congress of Thermal Sciences and Engineering,
ENCIT 2004, Rio de Janeiro, Brazil, Nov. 29 – Dec. 03, Invited Lecture – CIT04 –
IL04.
KATZ, J., PLOTKIN, A. (1991), Low Speed Aerodynamics: from Wing Theory to Panel
Methods, McGraw Hill Inc..
KEMPKA, S. N., STRICKLAND, J. H. (1993), A Method to Simulate Viscous Diffusion of
Vorticity by Convective Transport of Vortices at a Non-Solenoidal Velocity, SAND93-
1763 Report, Sandia Laboratories.
KOUMOUTSAKOS, P. D. (1993), Direct Numerical Simulations of Unsteady Separated
Flows Using Vortex Methods, Ph.D. Thesis, California Institute of Technology.
KRASNY, R. (1983), A Numerical Study of Kelvin-Helmholtz Instability by the Point Vortex
Method, Ph.D. Thesis, University of California Berkley.
KRASNY, R. (1986a), “A Study of Singularity Formation in a Vortex Sheet by the Point-
Vortex Approximation”,
J. Fluid Mech., 167:65-93.
KRASNY, R. (1986b), “Desingularization of Periodic Vortex Sheet Roll-Up”, J. Comp.
Phys., 65:292-313.
KUNDU, P. K. (1990), Fluid Mechanics, Academic Press.
LEI, C., CHENG, L., KAVANAGH, K. (1999), “Re-Examination of the Effect of a Plane
Boundary on Force and Vortex Shedding of a Circular Cylinder”, Journal of Wind
Engineering and Industrial Aerodynamics, v 80, pp 263-286.
98
LEONARD, A. (1980), “Vortex Methods for Flow Simulation”, Journal of Computational
Physics, v 37, pp 289-335.
LEONARD, A. (1985), “Computing Three-Dimensional Incompressible Flows with Vortex
Elements”, Ann. Rev. Fluid Mech., 17:523-559.
LEWIS, R. I. (1986), “Recent Developments and Engineering Applications of the Vortex
Cloud Method”. Presented to the First World Congress on Computational Mechanical,
Austin, Texas and also Computer Methods in Applied Mechanics and Enginnering, v
64, pp 153-176 (1987).
LEWIS, R. I. (1991), Vortex Element Method for Fluid Dynamic Analysis of Engineering
Systems, Cambridge Univ. Press, Cambridge, England, U.K.
LIN, C., LIN, W-J., LIN, S-S. (2005), “Flow Characteristics around a Circular Cylinder
Near a Plane Boundary”, Proceedings of the Sixteenth International Symposium on
Transport Phenomena, ISTP-16, Prague, Czech Republic, 29 August – 1 September,
(CD-Rom).
LIU, C. H., DOORLY, D. J. (2000), “Vortex Particle-In-Cell Method for Three-
Dimensional Viscous Unbounded Flow Computations”, Int. J. Num. Meth. In Fluids,
32:29-50.
LIU, C. H. (2001), “A Three-Dimensional Vortex Particle-In-Cell Method for Vortex
Motions in the Vicinity of a Wall”,
Int. J. Num. Meth. In Fluids, 37:501-523.
MILNE-TOMPSON, L. M. (1955), Theoretical Hydrodynamics, Macmillan & Co, London.
MOURA, W. H. (2007), Análise do Escoamento ao Redor de um Corpo Oscilante na
Presença de uma Superfície Plana, Dissertação de Mestrado, IEM/UNIFEI.
MUSTTO, A. A. (1998), Simulação Numérica do Escoamento em torno de um Cilindro
Circular com e sem Rotação Utilizando o Método de Vórtices
, Dissertação de Mestrado,
PEM – COPPE/UFRJ.
99
MUSTTO, A. A., HIRATA, M. H., BODSTEIN, G. C. R. (1998), “Discrete Vortex Method
Simulation of the Flow Around a Circular Cylinder with and without Rotation”,
Proceedings of the 16
th
A.I.A.A. Applied Aerodynamics Conference, Albuquerque, NM,
USA, June, A.I.A.A. Paper 98-2409.
MUSTTO, A. A. (2004), Simulação Numérica do Escoamento Turbulento em torno de um
Cilindro Circular Via o Método de Vórtices, Tese de Doutorado, PEM – COPPE/UFRJ.
NISHINO, T. (2007), Dynamics and Stability of Flow Past a Circular Cylinder in Ground
Effect
, Ph.D. Thesis, Faculty of Engineering, Science and Mathematics, University of
Southampton, 145p.
OGAMI, Y., AKAMATSU, T. (1991), “Viscous Flow Simulation Using the Discrete Vortex
Model – The Diffusion Velocity Method”, Computers & Fluids, v 19, pp 433-441.
OGAMI, Y. (2001), “Simulation of Heat-Fluid Motion by Vortex Method”, International
Journal, Series B, v 44, n 4, pp 513-519.
PANTON, R. L. (1984), Incompressible Flow, John Wiley & Sons.
PRÄGER, W. (1928), “Die Druckverteilung na Körpern in ebener PotentialStrömung”,
Physik. Zeitschr., XXIX: 865.
PRICE, S. J., SUMNER, D., SMITH, J. G., LEONG, K., PAIDOUSSIS, M. P. (2002),
“Flow Visualization around a Circular Cylinder Near to a Plane Wall, Journal of Fluids
and Structures
, v 16, pp 175-191.
PUCKETT, E. G. (1993), An Introduction and Survey of Selected Research Topics. In M. D.
Gunzburger and R. A. Nicolaides, editors, Incompressible Computational Fluid
Dynamics: Trends and Advances, Cambridge University Press, pages 335-408.
RECICAR, J. N. (2007), Oscilações de Grandes Amplitudes num Corpo que se Move com
Velocidade Constante
, Dissertação de Mestrado, IEM/UNIFEI.
RICCI, J. E. R. (2002), Simulação Numérica do Escoamento ao redor de um Corpo de
Forma Arbitrária, Estacionado nas Imediações de uma Superfície Plana, com o
100
Emprego do Método de Vórtices, Tese de Doutorado, Departamento de Mecânica,
Escola Federal de Engenharia de Itajubá.
ROSENHEAD, L. (1931), “The Formation of Vortices from a Surface of Discontinuity”,
Proc. R. Soc. Lond. A, A134:170-192.
ROSHKO, A., STEINOLFSON, A., CHATTOORGOON, V. (1975), “Flow Forces on a
Cylinder Near a Wall or Near Another Cylinder”, Proceedings of the 2
nd
US National
Conference on Wind Engineering Research, Fort Collins, Paper IV-15.
ROSSI, L. F. (1996), “Resurrecting Core Spreading Vortex Methods: A New Scheme that is
both Deterministic and Convergent”,
SIAM J. Sci. Comput., 17:370-397.
ROSSI, L.F. (1997), “Merging Computational Elements in Vortex Simulations”, SIAM J. Sci.
Comput., 18:1014-1027.
ROSSI, L. F. (2005), “Achieving High Order Convergence Rates with Deforming Basis
Functions”, SIAM J. Sci. Comput., v 26, n 3, pp 885-906.
ROSSI, L. F. (2006), “A Comparative Study of Lagrangian Methods Using Axisymmetric
and Deforming Blobs”, SIAM J. Sci. Comput., v 27, n 4, pp 1168-1180.
ROTTMAN, J. W., SIMPSON, J. E., STANSBY, P. K. (1987), “The motion of a cylinder
of Fluid Released from Rest in a Cross Flow”,
J. Fluid Mech., 177:307-337.
SALMON, J. K., WARREN, M. S., WINCKELMANS, G. S. (1994), “Fast Parallel Tree
Codes for Gravitational and Fluid Dynamical N-Body Problems”,
Intl. J. Supercomput.
Appl. High Perf. Comp., 8(2):129-142.
SARPKAYA, T. (1989), “Computational Methods with Vortices”, J. Fluids Eng., 11:5-52.
SCHLICHTING, H. (1979), Boundary-Layer Theory, 7
th
edition, McGraw-Hill.
SETHIAN, J. A. (1991), “A Brief Overview of Vortex Method, Vortex Methods and Vortex
Motion”, SIAM Philadelphia, pp 1-32.
101
SHANKAR, S., VAN DOMMELEN, L. (1996), “A New Diffusion Procedure for Vortex
Methods”, Journal of Computational Physics, v 127, pp 88-109.
SHERMAN, F. S. (1990), Viscous Flow, McGraw-Hill, International Editions Mechanical
Engineering Series.
SHINTANI, M., AKAMATSU, T. (1994), “Investigation of Two Dimensional Discrete
Vortex Method with Viscous Diffusion Model”, Computational Fluid Dynamics
Journal, v 3, n 2, pp 237-254.
SILVA DE OLIVEIRA, L., ALCÂNTARA PEREIRA, L. A., HIRATA, M. H. (2005),
“Numerical Simulation of Airfoil-Vortex Cloud Interaction in Ground Effect”,
Proceedings of the 18
th
International Congress of Mechanical Engineering, COBEM
2005, Ouro Preto, Brazil, Nov. 6-11, Paper COB-0106.
SMITH, P. A., STANSBY, P. K. (1988), “Impulsively Started Flow around a Circular
Cylinder by the Vortex Method”, J. Fluid Mech., 194:45-77.
SON, J. S., HANRATTY, T. J. (1969), “Velocity Gradients at the Wall for Flow Around a
Cylinder at Reynolds Number from to”, J. Fluid Mech, v 35 (2), pp 353-368.
SPALART, P. R. (1988), Vortex Methods for Separated Flows, Technical Memorandum
100068, NASA.
STOCK, M. J. (2007), Summary of Vortex Methods Literature (A lifting document rife with
opnion)
, April, 18.
TAKEDA, K., TUTTY, O. R., NICOLE, D. A. (1999), “Parallel Discrete Vortex Method on
Commodity Supercomputers; an Investigation into Bluff Body Far Wake Behaviour”,
ESAIM:
Proceedings, Third International Workshop on Vortex Flows and Related
Numerical Methods, v 7, pp 418-428.
TANEDA, S. (1965), “Experimental Investigation of Vortex Streets”, Journal of the Physical
Society of Japan, v 20, pp 1714-1721.
102
TANIGUCHI, S., MIYAKOSHI, K. (1990), “Fluctuating Fluid Forces Acting on a Circular
Cylinder and Interference with a Plane Wall”, Experiments in Fluids, v 9, pp 197-204.
TRITTON, D. J. (1988), Physical Fluid Mechanics, 2
nd
edition, Oxford Univ. Press.
UHLMAN, J. S. (1992), “An Integral Equation Formulation of the Equation of Motion of an
Incompressible Fluid”, Naval Undersea Warfare Center T. R., pp 10-86.
VAN DYKE, M. (1982), An Album of Fluid Motion, Parabolic Press, Stanford.
WHITE, F. M. (2002), Mecânica dos Fluidos, McGraw-Hill, 4ª edição, 570p.
YANG, R. J., HUANG, Y. G. (1999), “Simulation of High Reynolds Number Flow Over a
Flat Plate by a Deterministic Vortex Method”,
Comput. Fluid Dynamics Journal, v 8, n
2.
ZDRAVKOVICH, M. M. (1985a), “Observation of Vortex Shedding behind a Towed
Circular Cylinder Near a Wall, Flow Visualization III: Proceedings of the Third
International Symposium on Flow Visualization, ed. W. J. Yang, Hemisphere,
Washington DC, pp 423-427.
ZDRAVKOVICH, M. M. (1985b), “Forces on a Circular Cylinder Near a Plane Wall”,
Applied Ocean Research, v 7, pp 197-201.
ZDRAVKOVICH, M. M. (1987), “The Effects of Interference Between Circular Cylinders
in Cross Flow”, Journal of Fluids and Structures, v 1, pp 239-261.
ZDRAVKOVICH, M. M. (2003), Flow around Circular Cylinders: Vol 2: Applications,
Oxford University Press, Oxford, UK.
103
Apêndice A
DISTRIBUIÇÃO DA VORTICIDADE E DA
VELOCIDADE INDUZIDA
A.1 – O VÓRTICE POTENCIAL
O vórtice potencial tem o seguinte potencial complexo (Batchelor, 1967):
()
θ
2π
Γ
lnr
2π
iΓ
lnre
2π
iΓ
lnZ
2π
iΓ
Zf
θi
=== (A.1)
A partir da Equação A.1 define-se a função potencial de velocidades para um vórtice
girando no sentido horário, como sendo do tipo:
θ
2π
Γ
Φ = (A.2)
Na Figura A.1 tem-se um vórtice j, de intensidade positiva
j
Γ , localizado no ponto
(
)
jjj
y,xP = e um vórtice k, de intensidade negativa
k
Γ , localizado no ponto
()
kkk
y,xP = .
A velocidade que este vórtice j induz no ponto
(
)
kkk
y,xP
=
possui apenas um
componente tangencial, definido por:
kj
j
kj
kj
θ
r
1
2π
Γ
θ
Φ
r
1
u =
= (A.3)
104
Figura A.1 – Velocidade induzida (retirada de Alcântara Pereira, 1999).
A função potencial complexo para o conjunto escoamento uniforme (
°= 0α ), vórtice j e
vórtice k, vale:
() ()
[]
()
[]
kZZln
2π
Γ
ijZZln
2π
Γ
iUZeZf
k
j
iα
++=
(A.4)
Em termos dos componentes na direção x e na direção y, tem-se:
(
)
() ()
()
()()
()()
()
()()
()
()()
()
()()
()
()()
+
+
+
+
+
+
+
+
+=
=
+
+
+
+==
+
+==
2
k
2
k
kk
2
j
2
j
jj
2
k
2
k
kk
2
j
2
j
jj
kjkj
kk
k
jj
j
k
j
kjkj
yyxx
xx
2π
Γ
yyxx
xx
2π
Γ
i
yyxx
yy
2π
Γ
yyxx
yy
2π
Γ
Uivu
yyixx
1
2π
Γ
i
yyixx
1
2π
Γ
iU
dZ
Zdf
kZZ
1
2π
Γ
i
jZZ
1
2π
Γ
iU
dZ
Zdf
ivu
(A.5)
Para calcular a velocidade induzida no vórtice localizado no ponto k, faz-se
k
xx
=
e
k
yy = .
Deste modo:
(
)
()()
2
jk
2
jk
jkj
kj
yyxx
yy
2π
Γ
Uu
+
+=
(A.6)
(
)
()()
2
jk
2
jk
jkj
kj
yyxx
xx
2π
Γ
0v
+
=
(A.7)
U
=
1
y
Γ
k
<
0
u
kj
k
α
= 0
D
r
kj
v
kj
θ
Γ
j
> 0
θ
u
kj
θ
j
x
105
Como se pode notar, o vórtice potencial apresenta um comportamento singular na sua
origem, ou seja, onde ele está localizado. Com isto, na utilização de dois vórtices potenciais
muito próximos podem surgir problemas de instabilidade numérica, já que o campo de
velocidades é singular para
0r . Na Figura A.2 pode-se notar graficamente que para
vórtices potenciais muito próximos, o valor de
kj
r
torna-se pequeno, fazendo com que a
velocidade tenda para um valor infinito. Estas dificuldades impedem a efetiva utilização do
modelo do vórtice potencial.
Figura A.2 – Velocidade induzida pelo vórtice potencial (retirada de Alcântara Pereira, 1999).
Neste trabalho, para eliminar tal singularidade, o modelo do vórtice de Lamb, que
possui um núcleo viscoso com uma distribuição de vorticidade gaussiana no seu interior e
velocidade finita para todos os valores de r, mostra-se apropriado; veja a seguir.
A.2 – O VÓRTICE DE LAMB
O vórtice de Lamb possui distribuições de vorticidade (ω ) e velocidade tangencial
induzida (
θ
u ), mais suas derivadas, contínuas em todo o domínio, porque é solução da
equação da difusão da vorticidade:
=
r
ω
r
rr
ν
t
ω
(A.8)
A solução da Equação A.8, em uma região infinita (Kundu, 1990) é dada por:
0.00 1.00 2.00 3.00 4.00 5.00
0.00
0.40
0.80
1.20
1.60
u
kj
θ
r
kj
106
()
=
2
2
kj
2
j
σ
r
exp
πσ
Γ
tr,ω (A.9)
Na Equação A.9
σ é expresso como:
4tνσ = (A.10)
O componente tangencial da velocidade induzida pelo vórtice de Lamb, com a
distribuição de vorticidade acima, é definido pela seguinte equação:
=
2
2
kj
kj
j
kj
θ
σ
r
exp1
r
1
2π
Γ
u
(A.11)
O ponto
máx
r , onde a velocidade
(
)
ru
kj
θ
é máxima, é encontrado derivando-se a Equação
A.11 em relação à
kj
r
e igualando-a a zero (Mustto, 1998):
1,12091
σ
r
kj
= 1,12091σr
máx
= (A.12)
Deste modo:
=
2
máx
2
kj
kj
j
kj
θ
r
r
1,25643exp1
r
1
2π
Γ
u
(A.13)
Para
máx
rr = , tem-se:
π2r
Γ
0,71533u
máx
j
kj
θ
máx
= (A.14)
O raio do núcleo do vórtice de Lamb, definido de modo que a diferença entre as
velocidades induzidas calculadas pelo vórtice de Lamb e pelo vórtice potencial seja pequena,
vale:
máx0
2rσ = (A.15)
Nesta situação a diferença é de 0,6%.
107
Resolvendo-se a Equação A.12 com auxílio da Equação A.10, adimensionalizada, e
levando-se em consideração um incremento de tempo
t , encontra-se:
Re
t
2,24182r
máx
(A.16)
A equação final para o cálculo do raio do núcleo do vórtice de Lamb é obtida
relacionando-se as Equações A.15 e A.16:
Re
t
4,48364σ
0
= (A.17)
Para calcular a velocidade induzida pelo vórtice de Lamb, com
()
0θθ
σuu = , substitui-
se a Equação A.15 na Equação A.13:
=
2
0
2
kj
kj
j
kj
θ
σ
r
5,02572exp1
r
1
2π
Γ
u
(A.18)
O modelo do vórtice de Lamb, Figura A.3, não apresenta os problemas de
singularidade, mas na sua equação nota-se a presença do exponencial, que tem um cálculo
computacional demorado. Deste modo, todos os vórtices existentes na nuvem são vórtices de
Lamb inicialmente; a Equação A.18 pode ser usada quando dois vórtices estiverem muito
próximos um do outro, de modo que
0kj
σr
<
. Quando
0kj
σr , a Equação A.3 pode ser
aplicada sem os problemas de singularidade.
Figura A.3 – Velocidade induzida pelo vórtice de Lamb (retirada de Alcântara Pereira, 1999).
1. 00 2. 00 3. 00 4. 00 5. 0
0
0.04
0.08
r
kj
u
kj
θ
108
Apêndice B
RESULTADOS ADICIONAIS DA SIMULAÇÃO
NUMÉRICA
Tabela B.1 – Valores médios dos coeficientes de arrasto e de sustentação e número de
Strouhal para um cilindro circular submetido ao efeito solo para vários valores da relação
dh .
dh
D
C
L
C
St
0,55 1,15 0,41 0,11
0,60 0,83 0,38 0,11
0,65 1,29 0,06 0,16
0,70 1,38 -0,08 0,19
0,75 1,41 -0,07 0,20
0,80 1,39 -0,05 0,21
0,85 1,41 -0,13 0,21
0,90 1,42 -0,06 0,21
0,95 1,40 -0,01 0,21
1,00 1,39 -0,03 0,21
1,10 1,38 0,04 0,21
1,30 1,36 0,04 0,21
1,50 1,35 0,02 0,21
2,00 1,28 -0,01 0,21
2,50 1,27 0,10 0,21
Presente Simulação: (Euler, 300mb
=
, 0,05t
=
, 0,001epsσ
0
=
=
,
5
10Re =
).
109
(a)
0,55dh = (b) 0,60dh =
(c)
0,65dh = (d) 0,70dh =
(e)
0,75dh =
(f)
0,80dh =
Figura B.1 – Evolução das cargas aerodinâmicas ao longo do tempo para um cilindro circular
estacionado a pequenas distâncias da superfície plana móvel (Euler,
300mb = ,
0,05t
=
,
0,001epsσ
0
=
=
,
5
10Re = ).
110
(a)
0,85dh = (b) 0,90dh =
Figura B.2 – Evolução das cargas aerodinâmicas ao longo do tempo para um cilindro circular
estacionado a distâncias intermediárias da superfície plana móvel (Euler,
300mb = ,
0,05t = , 0,001epsσ
0
=
=
,
5
10Re =
).
(a)
0,55dh
=
(b)
0,75dh
=
(c)
30,1dh
=
Figura B.3 - Posição dos vórtices na esteira ao final da simulação numérica do cilindro
circular submetido ao efeito solo para alguns valores de
dh (Euler, 300mb = , 0,05t
=
,
0,001epsσ
0
=
=
,
5
10Re = ).
Livros Grátis
( http://www.livrosgratis.com.br )
Milhares de Livros para Download:
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas
Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo