ads:
MINISTÉRIO DA EDUCAÇÃO
UNIVERSIDADE FEDERAL DO RIO GRANDE DO SUL
PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA MECÂNICA
DESENVOLVIMENTO DE UM SISTEMA DE DINÂMICA DOS FLUIDOS
COMPUTACIONAL EMPREGANDO O MÉTODO DE ELEMENTOS FINITOS E
TÉCNICAS DE ALTO DESEMPENHO
por
João Américo Aguirre Oliveira Jr.
Dissertação para obtenção do Título de
Mestre em Engenharia
Porto Alegre, outubro de 2006
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
DESENVOLVIMENTO DE UM SISTEMA DE DINÂMICA DOS FLUIDOS
COMPUTACIONAL EMPREGANDO O MÉTODO DE ELEMENTOS FINITOS E
TÉCNICAS DE ALTO DESEMPENHO
por
João Américo Aguirre Oliveira Jr.
Engenheiro Mecânico
Dissertação submetida ao Corpo Docente do Programa de Pós-Graduação em
Engenharia Mecânica, PROMEC, da Escola de Engenharia da Universidade Federal do Rio
Grande do Sul, como parte dos requisitos necessários para a obtenção do Título de
Mestre em Engenharia
Área de Concentração: Fenômenos de transporte
Orientador: Profª. Drª. Adriane Prisco Petry
Comissão de Avaliação:
Prof. Dr. Argimiro Resende Secchi
Prof. Dr. Armando Miguel Awruch
Prof. Dr. Horácio Antônio Vielmo
Prof. Dr. Sérgio Luiz Frey
Prof. Dr. Flávio José Lorini
Coordenador do PROMEC
Porto Alegre, 06 de outubro de 2006.
ads:
RESUMO
O presente trabalho apresenta o desenvolvimento de um código numérico computacional,
baseado no método de elementos finitos, para simulação de grandes escalas de escoamentos bi- e
tridimensionais, transientes, incompressíveis, isotérmicos e turbulentos.
O código aproxima as equações médias espaciais de Navier-Stokes e da continuidade
escritas utilizando a hipótese de quase-incompressibilidade. O uso dessa hipótese permite a
aplicação de aproximações de mesma ordem para os campos de pressão e velocidade. Junto com
a hipótese de quase-incompressibilidade, é implementada uma forma proposta pelo grupo de
diagonalização das matrizes formadas no problema, de forma a simplificar e estabilizar a solução
do sistema linear.
Para a modelagem sub-malha é utilizado o modelo clássico de Smagorinsky para
contabilizar os efeitos das menores estruturas turbulentas (estruturas sub-malha) sobre o
escoamento resolvido (com grandes escalas turbulentas).
Toda a abordagem numérica utilizada é detalhada e apresentada na forma em que foi
implementada no código. Em resumo é utilizado o método de elementos finitos, aplicando os
métodos de Newton ou de Picard para aproximar o sistema não linear de equações e um esquema
de avanço no tempo variável, permitindo o uso das formas explícita e semi-implícita de avanço.
Uma importante parte do trabalho é dedicada à estruturação e aplicação de técnicas de alto
desempenho na elaboração do código. O objetivo é obter um código que permita melhorias e
implementação de novos métodos e modelos de forma simples e que, ao mesmo tempo, tenha um
bom desempenho computacional. São aplicados métodos de armazenagem compacta de matrizes
para a redução da quantidade de memória necessária e para viabilizar o uso de esquema semi-
implícitos de avanço no tempo. Uma técnica de paralelização do processamento para permitir o
uso de clusters foi implementada, sendo a mesma adequada para uso em máquinas de memória
compartilhada.
Dois casos de teste são apresentados. O escoamento em uma cavidade, bi- e tridimensional,
em regimes laminar e turbulento e o escoamento bidimensional, também laminar e turbulento,
sobre um degrau (expansão em um canal).
Por fim são apresentadas as conclusões obtidas e sugestões para melhorias, novas
implementações e novas linhas oriundas desse trabalho são feitas.
iii
ABSTRACT
“Development of a Computational Fluid Dynamics System Applying the Finite
Element Method and High Performance Techniques”
The present work presents the development of a computational code, based on the finite
elements method, for large eddy simulation of two and three dimensional, transient,
incompressible, isothermal and turbulent flows.
The code approximates the spatial average of the Navier-Stokes and continuity equations.
The hypothesis of slight-compressibility is used allowing the application of equal order
approximations for pressure and velocity fields. A form of matrix diagonalization is
implemented, in order to simplify and stabilize the solution of the linear system.
For the sub-grid modeling, the classical Smagorinsky model is employed, in order to
account for the effects of small turbulent structures (sub-grid structures).
The numerical approach employed in this work is detailed and presented in the form that
was implemented in the code. In summary, the finite elements method is employed, applying
Newton or Picard techniques to solve the non-linear system of equations and a variable time
marching scheme, allowing the use of explicit and semi-implicit forms schemes.
An important part of this work is dedicated to the structuring and application of high
performance techniques in the code elaboration. The objective is to obtain a code which allows
improvements and implementation of new methods and models in a simple way and conserving,
at the same time, a good computational performance. Methods of compact matrix storage are
applied to reduce the required memory when semi-implicit schemes are used. A parallelization
technique is also employed in order to enable the use of clusters, the implemented technique is
suitable to be used in shared memory machines.
Two test cases are presented. The two and three dimensional laminar and turbulent flow in
a lid-driven cavity and the two-dimensional laminar and turbulent flow over a backward-facing
step (a channel with an expansion).
Finally, the conclusions are presented, and suggestions for improvements, new
implementations and new subjects of future researches are given.
iv
ÍNDICE
1. Introdução....................................................................................................................................1
2. Equações básicas e Escoamentos turbulentos.............................................................................6
2.1. Equações básicas da dinâmica dos fluidos...........................................................................6
2.2. Formulação de quase-incompressibilidade ..........................................................................8
2.2.1. Equações básicas quase-incompressíveis na forma adimensional..............................10
2.3. Escoamentos turbulentos....................................................................................................12
2.3.1. Grandes escalas do escoamento turbulento e estruturas coerentes .............................15
3. Simulação numérica de escoamentos turbulentos.....................................................................17
3.1. Simulação numérica direta de escoamentos turbulentos....................................................18
3.2. Modelagem clássica de escoamentos turbulentos..............................................................19
3.3. Simulação de grandes escalas ............................................................................................21
3.3.1. Filtragem espacial para simulação de grandes escalas................................................23
3.3.2. Modelos sub-malha.....................................................................................................31
4. Método numérico ......................................................................................................................38
4.1. Forma variacional das equações de conservação de massa e quantidade de movimento ..40
4.2. Aproximação de elementos finitos e formação do problema matricial..............................42
4.3. Discretização do avanço no tempo.....................................................................................47
4.4. Método de Newton para sistemas de equações algébricas não-lineares ............................50
4.5. Oscilações no campo de pressão e diagonalização seletiva da matriz de massa................52
4.6. Método de Picard para sistemas de equações algébricas não-lineares...............................55
4.7. Método explícito simplificado e explícito diagonalizado ..................................................57
4.8. Solução do sistema linear...................................................................................................58
5. ASPECTOS COMPUTACIONAIS..........................................................................................59
5.1. Estruturação do novo código..............................................................................................59
5.2. Armazenagem compacta de matrizes.................................................................................60
5.3. Implementação de processamento paralelo utilizando a técnica OpenMP ........................63
6. RESULTADOS.........................................................................................................................65
6.1. Escoamento em uma cavidade bidimensional....................................................................65
6.2. Escoamento em uma cavidade tridimensional ...................................................................74
6.2.1. Escoamento em uma cavidade tridimensional a Re
cavidade
= 3200 ..............................74
6.2.2. Escoamento em uma cavidade tridimensional a Re
cavidade
= 10000 ............................82
v
6.3. Escoamento em um degrau bidimensional.........................................................................88
7. Conclusão..................................................................................................................................93
8. Referências bibliográficas.........................................................................................................95
vi
LISTA DE SÍMBOLOS
A Matriz de termos advectivos [-]
′
A
Matriz da derivada da matriz de termos advectivos [-]
B Dimensão da cavidade tridimensional no eixo de coordenadas x [m]
c Velocidade do som no meio [m/s]
C Termos cruzados da SGE [m
2
/s
2
]
C Parâmetro de viscosidade do modelo dinâmico [-]
C
Matriz de termos viscosos [-]
C
s
Constante de Smagorinsky [-]
D Dimensão da cavidade tridimensional no eixo de coordenadas y [m]
D
Matriz de termos do divergente de velocidade [-]
f Componente do vetor de forças externas por unidade de massa [N/kg]
F
Vetor de forças externas por unidade de volume [N/m
3
]
F
Vetor de carregamento (forças externas) [N]
F Componente do vetor de forças externas por unidade de volume [N/m
3
]
g Variável genérica [-]
G Função filtro (para simulação de grandes escalas) [-]
G
Matriz de termos do gradiente de pressão [-]
h Tamanho característico da malha [m]
h Altura do degrau [m]
H Altura do canal de entrada [m]
H Dimensão da cavidade tridimensional no eixo de coordenadas z [m]
ija
Vetor de índices [-]
J
Matriz jacobiana, ou matriz tangente [-]
k Número de onda do espaço de Fourier [m
-1
]
K
Matriz de rigidez não linear [-]
l Comprimento característico das menores escalas turbulentas [m]
l Comprimento da seção de entrada do degrau [m]
L Comprimento do canal após o degrau [m]
L Comprimento característico da geometria [m]
L Termos de Leonard da SGE [m
2
/s
2
]
Mc Número de Mach [-]
vii
M
Matriz de massa [-]
n número de variáveis [-]
N Número de nós que compõem o elemento [-]
p Pressão estática [Pa]
q Função peso associada à pressão [-]
Re Número de Reynolds [-]
R
Vetor de valores residuais da equação matricial aproximada do problema [-]
sa
Vetor de valores [-]
S Razão entre a altura do canal de entrada (H) e a altura do degrau (h)
S
Tensor taxa de deformação de grandes escalas [-]
S
Tensor taxa de cisalhamento [-]
S
Vetor de forças de tração sobre a fronteira do elemento [N]
t Tempo [s]
T Temperatura [K] ou tempo característico [s]
T Tensor de tensões sub-malha [Pa]
u
Vetor velocidade [m/s]
u Componente genérica ou componente na direção x do vetor velocidade [m/s]
U Velocidade de referência do problema [m/s]
U Velocidade média do perfil parabólico na entrada do canal [m/s]
U
Vetor de incógnitas [-]
v Componente genérica ou componente na direção y do vetor de velocidade [m/s]
V
Matriz de termos viscosos volumétricos [-]
w Componente genérica ou componente na direção z do vetor de velocidade [m/s]
x Direção genérica, direção x do sistema cartesiano de coordenadas ou componente
genérica do vetor de posição [m]
x
R
Distância de recolamento da recirculação principal [m]
x
Vetor posição [m]
y Direção genérica, direção y do sistema cartesiano de coordenadas ou componente
genérica do vetor de posição [m]
z Direção genérica, direção z do sistema cartesiano de coordenadas ou componente
genérica do vetor de posição [m]
∆t Passo de tempo [s]
w Função peso associada à velocidade [-]
viii
γ
Constante da aproximação dos termos de Leonard e Cruzados [-]
δ
Função delta de Kroneker [-]
Γ Contorno do domínio do problema [-]
λ Viscosidade volumétrica do fluido [Pa·s]
µ Viscosidade absoluta do fluido [Pa·s]
ν
Viscosidade cinemática do fluido [m
2
/s]
θ
Parâmetro de avanço no tempo [-]
ρ
Massa específica [kg/m
3
]
∆
Dimensão característica de um filtro [m]
∆
Tamanho característico de um filtro de teste [m]
τ Tensões sub-malha por unidade de massa específica da SGE [m
2
/s
2
]
ω
ωω
ω Vorticidade [s
-1
]
ψ
Função de forma [-]
Ω Domínio do problema [-]
Subscritos
0
Valor inicial ou valor definido de uma variável
ef
Valor efetivo de uma variável
i, j, k
Componentes das grandezas vetoriais e matriciais
c
Variável característica de uma estrutura coerente
D
Variável com condição de primeira espécie (condição de Dirichlet)
L
Variável baseada no comprimento característico da geometria L
n, n+1
Indicador do instante de tempo de cada equação
N
Variável com condição de segunda espécie (condição de Neumann)
RMS
Valor flutuante médio da variável
( )
S
Processo à entropia constante
λ
Variável baseada na viscosidade volumétrica do fluido
T
Variável associada localmente ao escoamento
T
Variável associada à turbulência
ix
Sobrescritos
T
Transposto (matriz ou vetor)
e
Representa o elemento
m
Variável associada ao nó m do elemento
p, p+1
Indicador do número da iteração no método de Newton
•
Derivada no tempo da incógnita
Variável média temporal ou filtrada (média espacial)
ɶ
Variável associada ao filtro teste
′
Valor flutuante da variável
*
Variável adimensional
ˆ
Variável definida no espaço de Fourier
Variável referente à matriz totalmente diagonalizada
ɶ
Variável referente à matriz diagonalizada seletivamente
⌢
Variável referente à matriz diagonalizada na pressão
x
ÍNDICE DE FIGURAS
Figura 3.1. a) Forma da função filtro tipo “caixa”; b) Forma da função filtro gaussiana.............25
Figura 4.1. Esquema do domínio Ω e das fronteiras Γ
D
e Γ
N
para o método de elementos finitos.
.......................................................................................................................................................39
Figura 5.1. Esquema da estrutura do programa dividido por módulos.........................................60
Figura 5.2. Trecho do código mostrando as diretivas de compilação in-line para o uso de
paralelização..................................................................................................................................64
Figura 6.1. Esquema com descrição do problema da cavidade bidimensional, valores de
condições de contorno, iniciais e propriedades do fluido. ............................................................66
Figura 6.2. Malha mais grosseira (25x25 elementos) e malha mais refinada (200x200 elementos)
utilizadas no caso da cavidade para comparação. .........................................................................67
Figura 6.3. Perfil de velocidade u* (adimensional) ao longo da linha vertical central da cavidade
(em x = 0.5m), comparação com outros resultados publicados. Re
cavidade
= 100. .........................68
Figura 6.4. Perfil de velocidade v* (adimensional) ao longo da linha vertical central da cavidade
(em y = 0.5m), comparação com outros resultados publicados. Re
cavidade
= 100. .........................69
Figura 6.5. Perfil de velocidade u* (adimensional) ao longo da linha vertical central da cavidade
(em x = 0.5m), comparação com outros resultados publicados. Re
cavidade
= 400. .........................70
Figura 6.6. Perfil de velocidade v* (adimensional) ao longo da linha vertical central da cavidade
(em y = 0.5m), comparação com outros resultados publicados. Re
cavidade
= 400. .........................70
Figura 6.7. Perfil de velocidade u* (adimensional) ao longo da linha vertical central da cavidade
(em x = 0.5m), comparação com outros resultados publicados. Re
cavidade
= 1000. .......................71
Figura 6.8. Perfil de velocidade v* (adimensional) ao longo da linha vertical central da cavidade
(em y = 0.5m), comparação com outros resultados publicados. Re
cavidade
= 1000. .......................71
Figura 6.9. Campo de velocidade (módulo) para Re
cavidade
= 1000...............................................72
Figura 6.10. Isolinhas de pressão para Re
cavidade
= 1000...............................................................72
Figura 6.11. Linhas de corrente obtidas sobre o campo resolvido para Re
cavidade
= 1000.............73
Figura 6.12. Linhas de corrente publicadas por Elias et al., 2004, para Re
cavidade
= 1000. ...........73
Figura 6.13. Descrição do problema da cavidade 3D com domínio, condições de contorno e
propriedades do fluido...................................................................................................................75
Figura 6.14. Vista frontal e tridimensional da malha para o caso da cavidade 3D, Re
cavidade
=
3200...............................................................................................................................................76
xi
Figura 6.15. Perfil de velocidade u* (adimensional) ao longo da linha vertical central da
cavidade (em x = 0,5) no plano central, comparação com outros resultados publicados. Re
cavidade
= 3200............................................................................................................................................77
Figura 6.16. Perfil de velocidade v* (adimensional) ao longo da linha horizontal central da
cavidade (em y = 0,5) no plano central, comparação com outros resultados publicados. Re
cavidade
= 3200............................................................................................................................................77
Figura 6.17. Perfil de velocidade u’
RMS
* (adimensional) ao longo da linha vertical central da
cavidade (em x = 0,5) no plano central, comparação com outros resultados publicados. Re
cavidade
= 3200............................................................................................................................................78
Figura 6.18. Perfil de velocidade v’
RMS
* (adimensional) ao longo da linha vertical central da
cavidade (em y = 0,5) no plano central, comparação com outros resultados publicados. Re
cavidade
= 3200............................................................................................................................................79
Figura 6.19. Campo de velocidade média (magnitude do vetor) no plano central da cavidade (z =
0,5m). Re
cavidade
= 3200. ................................................................................................................80
Figura 6.20. Isolinhas de pressão estática média no plano central da cavidade, com escala de
valores reduzida entre -1,0⋅10
-4
Pa e 1,0⋅10
-4
Pa. Re
cavidade
= 3200.................................................80
Figura 6.21. Linha de corrente partindo do ponto central da cavidade colorida pelo módulo da
velocidade média (local). Re
cavidade
= 3200. ..................................................................................81
Figura 6.22. Isosuperfícies de pressão dentro da cavidade. Pressão estática com escala reduzida
entre -1,0⋅10
-4
Pa e 1,0⋅10
-4
Pa. Re
cavidade
= 3200............................................................................82
Figura 6.23. Isosuperfícies de pressão dentro da cavidade. Pressão estática com escala reduzida
entre -1,0⋅10
-4
Pa e 1,0⋅10
-4
Pa. Re
cavidade
= 10000..........................................................................83
Figura 6.24. Perfil de velocidade u* (adimensional) ao longo da linha vertical central da
cavidade (em x = 0,5) no plano central, comparação com outros resultados publicados. Re
cavidade
= 10000..........................................................................................................................................84
Figura 6.25. Perfil de velocidade v* (adimensional) ao longo da linha horizontal central da
cavidade (em y = 0,5) no plano central, comparação com outros resultados publicados. Re
cavidade
= 10000..........................................................................................................................................84
Figura 6.26. Perfil de velocidade u’
RMS
* (adimensional) ao longo da linha vertical central da
cavidade (em x = 0,5) no plano central, comparação com outros resultados publicados. Re
cavidade
= 10000..........................................................................................................................................85
xii
Figura 6.27. Perfil de velocidade v’
RMS
* (adimensional) ao longo da linha vertical central da
cavidade (em y = 0,5) no plano central, comparação com outros resultados publicados. Re
cavidade
= 10000..........................................................................................................................................85
Figura 6.28. Campo de velocidade média (magnitude do vetor) no plano central da cavidade (z =
0,5m). Re
cavidade
= 10000. ..............................................................................................................86
Figura 6.29. Isolinhas de pressão estática média no plano central da cavidade, com escala de
valores reduzida entre -2,0⋅10
-4
Pa e 2,0⋅10
-4
Pa. Re
cavidade
= 10000...............................................87
Figura 6.30. Linha de corrente partindo do ponto central da cavidade colorida pelo módulo da
velocidade média (local). Re
cavidade
= 10000. ................................................................................87
Figura 6.31. Isosuperfícies de pressão dentro da cavidade. Pressão estática com escala reduzida
entre -2,0⋅10
-4
Pa e 2,0⋅10
-4
Pa. Re
cavidade
= 10000..........................................................................88
Figura 6.32. Diagrama esquemático do problema do degrau, com dimensões, condições de
contorno e propriedades do fluido.................................................................................................89
Figura 6.33. Malha utilizada no caso do degrau bidimensional, região próxima ao degrau.........90
Figura 6.34. Campo de velocidades obtido para Re
degrau
= 800....................................................91
Figura 6.35. Campo de pressão estática obtido para Re
degrau
= 800..............................................91
Figura 6.36. Comparação dos pontos de recolamento obtidos com valores publicados por Armaly
et al., 1893.....................................................................................................................................91
xiii
1
1. INTRODUÇÃO
1.1. Apresentação e motivação
Cada vez mais o uso de ferramentas de simulação de processos fluidodinâmicos vem se
tornando importante no desenvolvimento de novos projetos industriais. Aplicações desse tipo de
ferramenta são conhecidas em diversos campos como, por exemplo, na indústria aeronáutica,
petrolífera, mecânica, civil, etc. Vários códigos comerciais para a simulação numérica de
processos estão disponíveis, cabendo ao engenheiro, através de seu conhecimento dos fenômenos
e diferentes modelos disponíveis, decidir quais ferramentas são mais ajustadas ao caso em
questão.
Assim o estudo de métodos numéricos e de dinâmica dos fluidos computacional, junto com
um conhecimento do problema físico, é de extrema importância para a qualidade da pesquisa
realizada.
No presente trabalho é apresentado um estudo sobre a modelagem numérica de
escoamentos turbulentos empregando o método de elementos finitos e a simulação de grandes
escalas. A partir deste estudo é elaborado um código numérico para a aproximação das equações
de conservação.
1.2. Objetivos
O principal objetivo deste trabalho é o desenvolvimento de um sistema de dinâmica dos
fluidos computacional que será utilizado como base para futuros trabalhos do grupo. Nesta
primeira etapa esse código deverá ser capaz de simular escoamentos bi- e tridimensionais,
transientes, isotérmicos, incompressíveis e turbulentos (via simulação de grandes escalas) de
fluidos newtonianos.
É necessário que este código permita alterações e implementação de novos métodos e
modelos de forma simples, porém ao mesmo tempo tenha um desempenho computacional
adequado ao estudo de problemas de escoamentos transientes e turbulentos sobre geometrias
complexas.
2
1.3. Estrutura do trabalho
O presente trabalho é apresentado utilizando a seguinte estrutura:
No capítulo 2 é detalhada uma revisão bibliográfica sobre a análise de escoamentos
turbulentos. São apresentadas as equações básicas da mecânica dos fluidos, a equação de
consevação de massa (equação da continuidade) e a equação de conservação da quantidade de
movimento para fluidos newtonianos (Equação de Navier-Stokes). Essas equações fundamentais
[Schlichting, 1979; White, 1999; Fox e McDonald, 2001] são baseadas na hipótese do contínuo e
descrevem o movimento de um fluido sob uma perspectiva euleriana.
A equação da conservação da massa é re-escrita utilizando a hipótese de quase-
incompressibilidade [Kawahara e Hirano, 1983; Sonntag et al., 1998; Zienkiewicz e Taylor,
2000] obtendo-se uma forma que é fisicamente mais coerente (pricipalmente para escoamentos
turbulentos) e com vantagens na sua implementação em códigos numéricos.
Uma discussão sobre escoamentos turbulentos, características gerais, definição,
apresentação dos conceitos físicos de escalas turbulentas, grandes escalas e estruturas coerentes
da turbulência também é feita nesse capítulo. Uma lista de referências bibliográficas sobre o
assunto é apresentada [Silva Freire et al., 2002; Schlichting, 1979; Hinze, 1975; Tennekes e
Lumley, 1994; Lesieur et al., 1995; Ferziger, 1993] mas sobre o assunto uma vastíssima gama de
trabalhos está disponível.
No capítulo 3 são exibidas as diferentes metodologias utilizadas para a simulação numérica
de escoamentos turbulentos. As três metodologias são:
Simulação numérica direta de escoamentos turbulentos, onde todas as escalas (de
comprimento e tempo) do escoamento turbulento são resolvidas na aproximação numérica das
equações de Navier-Stokes e da continuidade. Os resultados são mais completos e comparáveis a
qualquer dado experimental de qualidade. A grande dificuldade aqui é o imenso esforço
computacioal requerido, o qual inviabiliza a solução de escamentos muito complexos e a altos
valores do número de Renolds.
Modelagem clássica da turbulência, na qual o escoamento tubulento é dividido em uma
componente média no tempo e uma componente flutuante (decomposição de Reynolds). A parte
média no tempo é então resolvida e a parte flutuante é modelada, aqui entrando diversos modelos
de turbulência para resolver o problema de fechamento (aparecimento de incógnitas adicionais
nas equações médias nas quais está a contribuição da componente flutuante na componente
média). A grande maioria dos modelos de turbulência na abordagem clássica se baseia na
hipótese de Boussinesq (“viscosidade turbulenta”) e alguns trabalham com equações adicionais
3
de transporte das componentes do tensor de Reynolds (tensor de tensões turbulentas) [Wilcox,
1994].
Simulação de grandes escalas de escoamentos turbulentos, aqui apenas as maiores escalas
dos escoamentos turbulentos (mais fisicamente significativas na maioria dos fenômenos e
processos presentes nos escoamentos) são resolvidas e as pequenas escalas (que são mais
uniformes e isotrópicas, mais simples de modelar) são aproximadas. A separação das escalas é
equivalente à decomposição de Reynolds, porém ao invés do uso de médias no tempo se utilizam
filtros espaciais. A simulação de grandes escalas é mais cara computacionalmente em relação à
modelagem clássica e mais barata que a simulação numérica direta. Seus resultados são muito
bons e também bastante detalhados (resultados sempre ao longo do tempo, simulações não-
estacionárias) quando bem conduzidas. A maior quantidade de informação obtida (quando
comparada a simulação usando a metodologia clássica) permite a análise de fenômenos mais
complexos, como, por exemplo, interação fluido-estrutura (vibrações induzidas pelo
escoamento).
Em detalhes são apresetadas as diferentes abordagens para a simulação de grandes escalas
(métodos de filtragem espacial para divisão entre grandes e pequenas escalas e modelos de
turbulência sub-malha) utilizados. Boas revisões sobre simulação de grandes escalas de
escoamentos turbulentos são encontradas nos trabalhos de Rogallo e Moin, 1984, Lesieur et al.,
1995 e Piomelli, 1999, entre outros.
No capítulo 4 toda a modelagem numérica, baseada no método de elementos finitos
[Hughes, 1987; Zienkiewicz e Taylor, 2000; Reddy e Gartling, 1994] é mostrada. Todo o
processo de obtenção da forma de resíduos ponderados das equações de conservação (a partir de
suas formas adimensionais e quase-incompressíveis) é exibido. A partir dessa forma de resíduos
ponderados é feita discretização do domínio em um domíno de elementos finitos e a
aproximação de elementos finitos (na qual são definidas funções teste e funções de forma para
representar os diferentes campos dentro da discretização de elementos finitos) até chegar-se ao
problema matricial (ainda diferencial no tempo e não-linear) é definida.
O esquema de avanço no tempo utilizado para tornar o problema algébrico (esquema de
família θ) permite o uso de diferentes aproximações (totalmente implícito, semi-implícito e
totalemente explícito). Para solução do sistema de equações não-lineares são implementados os
métodos de Picard e Newton. É utilizada também a diagonalização seletiva das matrizes para
estabilização e aceleração da solução do problema, sendo uma das formas implementadas
sugeridas na literatura e outra proposta no trabalho. Por fim são citadas as diferentes soluções
implementadas para atacar o problema matricial algébrico e linear obtido.
4
No capítulo 5 são apresentados aspectos computacionais do trabalho elaborado. Questões
sobre o desenvolvimento do código computacional utilizado e sua estruturação, para futuros
incrementos e melhoias, são apresentadas e discutidas. Dois pontos dentro do código também são
focados, são eles a estrutura de armazenagem compacta de matrizes; utilizada para economia de
memória, redução do número de operações realizadas e conseqüente aumento na velocidade do
código; e a paralelização do processamento; utilziando a técnica OpenMP que possibilitou a
solução de problemas maiores e mais complexos com o mesmo código apenas pela distribuição
de carga computacional aproveitando todos os recursos disponíveis.
O capítulo 6 contém uma série de resultados de casos de referência obtidos utilizando o
código desenvolvido para validação do programa e início do estudo da turbulência utilizando a
simulação computacional.
Por fim no capítulo 7 são mostradas as conclusões e propostas para novos trabalhos e no
capitulo 8 está disponível a lista completa de referências bibliográficas citadas ao longo do texto.
1.4. Desenvolvimento do código
Como brevemente apresentado na seção anterior, o código desenvolvido utiliza o método
de elementos finitos para aproximação das equações de conservação. O método é aplicado sobre
as equações de balanço escritas utilizando a hipótese de quase-incompressibilidade, o que
permite o uso de funções de interpolação de mesma ordem para os campos de velocidade e
pressão dos escoamentos.
Para permitir as simulações bi- e tridimensionais, são criados elementos quadriláteros e
hexaedricos na biblioteca base do programa. Lembrando que outros elementos podem ser
facilmente implementados.
Para a simulação de escoamentos turbulentos a técnica utilizada de aproximação é a
simulação de grandes escalas, onde as maiores escalas turbulentas do escoamento (que carregam
a maior parte da energia e que promovem a parte significativa da mistura de propriedades do
escoamento) são resolvidas e apenas os efeitos das menores escalas (menores que a resolução da
malha utilizada) são modelados. Para modelar esses efeitos sub-malha o modelo de Smagorinsky
é utilizado.
O código é escrito em FORTRAN 95, utilizando técnicas de alto desempenho. Para
garantir sua portabilidade e viabilizar alterações ele é escrito de forma modular (cada módulo é
um conjunto de sub-rotinas armazenadas um arquivo diferente), isso permite também uma
compilação seletiva, acrescentando ou excluindo rotinas de acordo com o caso. Para permitir a
5
solução de grandes problemas são utilizadas técnicas de armazenagem compacta de matrizes e de
processamento em paralelo.
6
2. ANÁLISE DE ESCOAMENTOS TURBULENTOS
2.1. Equações básicas da dinâmica dos fluidos
Para estudar-se o problema de simulação numérica de escoamentos turbulentos, parte-se
sempre das equações básicas da mecânica dos fluidos [Schlichting, 1979; White, 1999; Fox e
McDonald, 2001], que são a equação da conservação de massa e a equação da conservação da
quantidade de movimento (com alguma relação constitutiva adicional).
As equações de conservação são deduzidas baseadas na hipótese do contínuo, na qual um
fluido é considerado uma distribuição contínua da matéria e todas as propriedades do fluido
podem ser expressas por funções contínuas das coordenadas espaciais e do tempo. As
propriedades do fluido são então representadas por campos (visão euleriana).
A equação da conservação de massa é obtida através de um balanço entre o fluxo de massa
que entra e que sai de um volume de controle diferencial e a variação da massa específica no
interior desse volume. Essa equação também é chamada de equação da continuidade, pois não
requer nenhuma hipótese para sua dedução exceto a do meio contínuo. Na forma apresentada na
equação (2.1) ela não leva em conta nenhuma fonte ou sumidouro de massa no interior do
volume de controle.
0 sobre , 0
D
t
Dt t
ρ ρ
ρ ρ
∂
+ ∇ = +∇ = Ω ≥
∂
u ui i (2.1)
onde t é o tempo,
ρ
é a massa específica e u o vetor velocidade.
A equação da conservação da quantidade de movimento é obtida por um balanço entre o
fluxo de quantidade de movimento (também chamado de momentum), que entra e que sai de um
volume de controle diferencial, e as forças externas (de corpo e de superfície) agindo sobre o
volume de controle. Essa relação é uma forma da segunda lei de Newton, onde as forças são
representadas pelas tensões sobre o volume de controle. Os valores das tensões são obtidos a
partir de relações constitutivas (de um fluido newtoniano no caso) envolvendo as deformações
sofridas pelo volume de controle e pelas propriedades do fluido.
Para um fluido newtoniano a relação entre as tensões e as deformações sofridas pelo
volume de controle é linear. Expressando então as relações através das propriedades do fluido e
dos gradientes de velocidade obtém-se a equação (2.2), também chamada de equação de Navier-
Stokes, expressa em termos de variáveis primárias e na forma mais geral.
7
( ) ( ) ( )
sobre , 0
T
D
p t
Dt
ρ
µ λ
= −∇ +∇ ∇ + ∇ +∇ ∇ + Ω ≥
u
u u u Fi i (2.2)
onde
µ
é a viscosidade absoluta e
λ
a viscosidade volumétrica do fluido, p a pressão estática e F
o vetor de forças externas por unidade de volume.
Nesta forma da equação de Navier-Stokes está mantido o termo do divergente de
velocidade que multiplica a viscosidade volumétrica. Para essa viscosidade volumétrica
normalmente se aplica a hipótese de Stokes, que estabelece uma relação desta com a viscosidade
absoluta (cisalhante) do fluido, esta relação é dada pela equação (2.3).
2
3
λ µ
= − (2.3)
Essa hipótese garante uma relação de igualdade entre a tensão normal sobre o fluido e a
pressão termodinâmica. Sendo então bastante conveniente, pouca discussão é feita a esse
respeito. A interpretação física da viscosidade volumétrica é que ela representa uma dissipação
de energia quando um volume de fluido é comprimido ou expandido devido unicamente à
tensões normais. Usualmente não se tem preocupação a respeito de seu significado e seus efeitos
quando nas investigações (numéricas e experimentais) se trabalha com a hipótese (ou condição)
de incompressibilidade, na qual o divergente de velocidade que aparece multiplicando essa
propriedade se anula.
Essas equações, exibidas na forma vetorial, são mais comumente utilizadas na área de
análise numérica na forma indicial [Petry, 2002], como mostradas nas equações (2.4) e (2.5).
0 ( 1,2,3), sobre , 0
j j
j j
u u
D
j t
Dt x t x
ρ
ρ ρ
ρ
∂ ∂
∂
+ = + = = Ω ≥
∂ ∂ ∂
(2.4)
( , , 1,2,3), sobre , 0
j i j
i i k
i
j i j i j i k
u u u
u u u
p
F i j k t
t x x x x x x x
ρ
ρ
µ λ
∂ ∂
∂ ∂ ∂
∂ ∂ ∂
+ = − + + + + = Ω ≥
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
(2.5)
onde u
i
, u
j
e u
k
são as componentes do vetor velocidade e F
i
as componentes do vetor de forças
externas por unidade de volume nas direções x
i
, x
j
ou x
k
, respectivamente.
8
2.2. Formulação de quase-incompressibilidade
A formulação de quase-incompressibilidade, também chamada de leve compressibilidade
(do inglês slight compressibility), se baseia na relação entre a pressão e a massa específica do
fluido. Essa formulação parte da suposição de que, nos escoamentos em questão, as variações de
massa específica são desprezíveis. Usualmente essa suposição (de escoamentos incompressíveis)
implica em que a velocidade de propagação de uma perturbação de pressão no meio (velocidade
do som) seja infinita, o que caracteriza uma condição idealizada. Na formulação de quase-
incompressibilidade trabalha-se com a hipótese de que apesar da massa específica do fluido
permanecer praticamente constante, uma perturbação de pressão se propaga com velocidade
finita no meio (velocidade do som finita); por isso essa formulação também é chamada de
incompressibilidade real.
O uso da formulação de quase-incompressibilidade é bastante conveniente para aplicação
no método de elementos finitos, pois, como discutido mais tarde, evita problemas de
aparecimento de zeros na diagonal principal do sistema linear formado. Isso é devido ao
aparecimento de um termo de derivada no tempo da pressão na equação da continuidade,
derivado como se segue.
Conforme explicitado na literatura [Kawahara e Hirano, 1983; Sonntag et al., 1998;
Zienkiewicz e Taylor, 2000] as equações de balanço são completadas através de relações de
estado para a obtenção das propriedades do fluido. Para a massa específica, a equação de estado
fornece:
(
)
,
p T
ρ ρ
=
(2.6)
Onde T é a temperatura local do fluido. Considerando um escoamento isotérmico e
invertendo-se a relação, chega-se à:
(
)
p p
ρ
=
(2.7)
Derivando-se essa relação e utilizando-se a regra da cadeia obtém-se:
p p
t t
ρ
ρ
∂ ∂ ∂
=
∂ ∂ ∂
(2.8)
9
Utilizando a definição de velocidade do som dada pela equação (2.9) e considerando a
variação de pressão causada por uma onda sonora isentrópica, obtém-se a relação da equação
(2.10).
2
constante
S
p
c
ρ
=
∂
=
∂
(2.9)
2
p
c
t t
ρ
∂ ∂
=
∂ ∂
(2.10)
onde c é a velocidade do som no meio. A hipótese de considerar-se o processo isentrópico é
justificada pelo fato de que uma onda sonora é uma variação de pressão muito pequena como
argumentado em Sonntag et al., 1998. Pela equação (2.9) conclui-se também que para um
escoamento incompressível a velocidade do som assume um valor constante.
Inserindo a relação (2.10) na equação da conservação de massa na forma da equação (2.4)
chega-se a:
2
1
0 ( 1,2,3)
j
j
u
p
j
c t x
ρ
∂
∂
+ = =
∂ ∂
(2.11)
Percebe-se que se considerarmos a velocidade do som como infinita a equação (2.11)
reduz-se à comumente utilizada equação da continuidade para fluidos incompressíveis, dada por:
0 ( 1,2,3)
j
j
u
j
x
∂
= =
∂
(2.12)
Para o escoamento transiente, isotérmico e quase-incompressível de um fluido newtoniano,
considera-se que as propriedades do fluido (massa específica, viscosidades e velocidade do som)
sejam constantes. Com isso chega-se nas formas finais (utilizadas no restante do trabalho) das
equações de conservação de massa, equação (2.13), e de quantidade de movimento, equação
(2.14).
2
0 ( 1,2,3)
j
j
u
p
c j
t x
ρ
∂
∂
+ = =
∂ ∂
(2.13)
10
1
j i j
i i k
i
j i j i j i k
u u u
u u u
p
f
t x x x x x x x
λ
ν
ρ ρ
∂ ∂
∂ ∂ ∂
∂ ∂ ∂
+ = − + + + +
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
(2.14)
Onde f
i
é a componente do vetor de forças externas por unidade de massa e ν a viscosidade
cinemática do fluido. Observa-se que o termo que multiplica a viscosidade volumétrica não
desaparece, já que a equação da continuidade na forma quase-incompressível não anula o
divergente de velocidade. Entretanto, como serão discutidos apenas escoamentos
incompressíveis, o valor da viscosidade volumétrica será tomado como nulo.
A forma mais comumente encontrada das equações de conservação é aquela para
escoamentos incompressíveis (considerando a massa específica do fluido, ρ, constante) mostrada
nas equações (2.15) e (2.16).
0
j
j
u
x
∂
=
∂
(2.15)
1
i i i
j i
j i j j
u u u
p
u f
t x x x x
ν
ρ
∂ ∂ ∂
∂ ∂
+ = − + +
∂ ∂ ∂ ∂ ∂
(2.16)
Onde a equação (2.16) foi obtida através do uso da equação (2.15) na conservação de
quantidade de movimento (2.5).
2.2.1.
Equações básicas quase-incompressíveis na forma adimensional
Para facilitar a análise de diferentes escoamentos através da variação de poucos parâmetros
se utiliza a análise dimensional. Através da adimensionalização das equações básicas se obtém
grupos adimensionais (variáveis) que caracterizam um grupo de escoamentos que apresentam a
mesma solução, guardadas as relações de semelhança. O uso das equações adimensionalizadas
também tem vantagens associadas aos métodos numéricos para solução, como será discutido
adiante.
A adimensionalização das equações (2.13) e (2.14) são baseadas nos parâmetros
característicos do problema (o comprimento característico da geometria L, a velocidade de
referência do problema U e a velocidade do som no meio c) e seguem as seguintes definições:
11
* * * * *
2
L
U
, , , ,
L L U U U
i i i
i i i
x u f
t p
x t u p f
c
ρ
= = = = = (2.17)
Onde t
*
e p
*
são o tempo e a pressão adimensionais e x
i
*
, u
i
*
e f
i
*
são as componentes
adimensionais do vetor posição, da velocidade e da força por unidade de massa respectivamente.
Aplicando essas definições nas equações (2.13) e (2.14) chega-se a:
*
*
* *
1
0
j
j
u
p
t Mc x
∂
∂
+ =
∂ ∂
(2.18)
* * *
* * *
*
*
* * * * * *
1 1 1
j i j
i i k
i
j i L i j k
u u u
u u u
p
f
t x Mc x Re x x Re x
λ
∂ ∂
∂ ∂ ∂
∂
+ = − + + + +
∂ ∂ ∂ ∂ ∂ ∂
(2.19)
Na equação (2.19) aparecem, além das variáveis adimensionais, três números
adimensionais que são o número de Mach (Mc) e o número de Reynolds baseado na viscosidade
cinemática e na dimensão característica (Re
L
) e baseado na viscosidade volumétrica e na
dimensão característica (Re
λ
), obtido pelas expressões (2.20), (2.21) e (2.22).
U
Mc
c
=
(2.20)
UL
L
Re
ν
= (2.21)
UL
Re
λ
ρ
λ
= (2.22)
O número de Mach expressa a relação entre a velocidade característica do escoamento e a
velocidade do som no fluido, utilizada como um parâmetro de compressibilidade. Para números
de Mach até 0,3 (para alguns autores 0,4) o escoamento é considerado incompressível, ou seja, as
variações de massa específica são desprezíveis. Para valores maiores efeitos de
compressibilidade não podem mais ser desprezados.
O número de Reynolds expressa a relação entre a inércia e as forças viscosas do
escoamento; esse valor é utilizado como referência para definir o regime dos escoamentos. Cada
12
tipo de escoamento (cada combinação geometria – condição de operação) possui um valor de
Reynolds crítico a partir do qual o escoamento transita para o regime turbulento. Este valor é
como uma média do valor de transição, já que essa transição não é instantânea. Seu início e fim
dependem das perturbações presentes no escoamento. O número de Reynolds baseado no
comprimento característico da geometria é utilizado para referência em trabalhos e para
comparação com valores experimentais. Por exemplo, para o escoamento em um duto circular o
número de Reynolds crítico é 2300, acima do qual o escoamento é considerado turbulento.
O número de Reynolds baseado na viscosidade volumétrica não tem grande significado no
problema em questão, já que o valor de
λ
será mantido nulo (por considerar o escoamento
incompressível). Existem outras representações do número de Reynolds, que serão discutidas no
momento apropriado.
2.3.
Escoamentos turbulentos
A grande maioria dos escoamentos de interesse em engenharia está no regime turbulento.
O escoamento do ar sobre a asa de um avião, da água em contato com o casco de um navio ou de
derivados de petróleo através de uma linha de transmissão são exemplos de situações onde o
regime turbulento é predominante. A turbulência se origina de pequenas perturbações presentes
nos escoamentos que são amplificadas devido a instabilidades que aparecem quando o
escoamento atinge altos números de Reynolds (entrando nas faixas consideradas de escoamentos
turbulentos). Tais perturbações se propagam e se degeneram fazendo com que o escoamento
entre no regime turbulento [Silva Freire et al., 2002].
Vários autores procuraram definir mais precisamente o conceito de turbulência
[Schlichting, 1979; Hinze, 1975; Tennekes e Lumley, 1994; Lesieur et al., 1995]. A idéia básica
é que o escoamento turbulento consiste no movimento de um fluido no qual uma flutuação
irregular é sobreposta a uma corrente principal. Sobre essa idéia se acrescentam os seguintes
pontos:
•
Irregularidade: a aleatoriedade das flutuações presentes no escoamento turbulento o torna
imprevisível, obrigando o uso de ferramentas estatísticas para a análise;
•
Difusividade aumentada: devido aos movimentos turbulentos as quantidades do escoamento
(como massa, quantidade de movimento, energia) são misturadas de forma mais eficiente;
esse efeito faz com que aparentemente a viscosidade (e outros coeficientes de difusão do
escoamento) seja aumentada;
13
•
Altos números de Reynolds: como já citado as perturbações que geram o escoamento
turbulento só se amplificam a altos números de Reynolds;
•
Flutuações tridimensionais de vorticidade: intensas flutuações de vorticidade estão presentes
nos escoamentos turbulentos, flutuações que não conseguem se manter em escoamentos
bidimensionais, o que força o escoamento turbulento a ser tridimensional por definição;
•
Alta dissipação de energia: o aumento nas taxas de cisalhamento interno do fluido, a constante
produção de energia cinética turbulenta e a dissipação de energia nas menores escalas torna o
escoamento turbulento altamente dissipativo; por esse mesmo motivo a turbulência necessita
de uma fonte constante de energia para ser mantida;
•
Fenômeno contínuo: embora as menores escalas do escoamento turbulento sejam, algumas
vezes, muitas ordens de grandeza inferiores às maiores escalas, elas continuam bem maiores
que as escalas moleculares, assim o fenômeno pode ser tratado como um problema contínuo
governado pelas equações da mecânica dos fluidos;
•
Escoamentos turbulentos são escoamentos e não fluidos: a turbulência é uma propriedade dos
escoamentos, não dos fluidos;
•
Larga faixa de escalas de comprimento e tempo: um escoamento turbulento deve envolver
vórtices de uma larga faixa de escalas de comprimento e tempo (ou números de onda num
domínio de Fourier).
As equações que descrevem escoamentos nesse regime são as mesmas para o escoamento
laminar, especificamente as equações da continuidade (2.4) e de Navier-Stokes (2.5). Devido a
sua complexidade essas equações possuem solução analítica apenas em casos com diversas
simplificações, sendo necessário então o uso de alguma ferramenta numérica para obter uma
solução em casos de interesse prático (mais complexos geométrica ou fisicamente). A
dificuldade aumenta ainda mais quando outros fenômenos que afetam a turbulência estão
presentes como em escoamentos com estratificação, com empuxo, rotacionais, com reações
químicas ou compressíveis [Ferziger, 1993].
A dificuldade nesse tipo de solução é devida justamente ao fato de que os escoamentos
turbulentos são compostos por movimentos em várias escalas de comprimento e tempo. As
maiores escalas de comprimento (L) da ordem da geometria do problema e as menores escalas (l)
estão na ordem da dissipação viscosa. Uma estimativa para a ordem de grandeza da escala
dissipativa e da relação entre as duas escalas são mostradas nas equações (2.23) e (2.24)
respectivamente [Silva Freire et al., 2002].
14
1
3
4
3
L
l
U
ν
=
(2.23)
3
4
L
l
L
Re
= (2.24)
A equação (2.24) mostra que a relação de escalas aumenta com o aumento do número de
Reynolds. Assim para uma mesma geometria as menores escalas ficam cada vez menores até o
ponto em que se torna inviável (para os recursos computacionais disponíveis atualmente)
simular-se numericamente todas as escalas presentes no escoamento devido à imensa capacidade
computacional requerida.
As diversas escalas de movimento presentes num escoamento interagem entre si; algumas
interações de vórtices são já bem conhecidas e também são de interesse particular como o
enrolamento de uma esteira de vórtices, o pareamento, formação de dipolo e vórtices tipo
grampo de cabelo. Seus papéis na transição e na influência nas propriedades do escoamento são
freqüentemente objeto de estudo. Uma das principais características dessas escalas é a “cascata
de energia” formada. Os maiores vórtices recebem energia do escoamento principal e passam
para os menores, que por sua vez a repassam aos menores ainda, até as escalas dissipativas, onde
essa energia é convertida em calor pelo atrito viscoso.
Apesar de que a transferência de energia ser normalmente das maiores para as menores
escalas, o fluxo inverso de energia (transferência de energia das menores escalas para as
maiores) ocorre em diversas situações, especialmente nas primeiras etapas da transição para o
regime turbulento em escoamentos parietais. Um estudo sobre esse efeito e uma melhor
descrição é apresentada no trabalho de Piomelli et al., 1991. Comumente o fenômeno da
transferência inversa de energia é citado pela expressão consagrada do inglês backscatter.
Por sua característica estocástica (randômica), muitas vezes sobreposta a estruturas
organizadas, e fortemente não-linear, muitos trabalhos tentaram gerar teorias estatísticas para
descrição dos escoamentos turbulentos. Todas as tentativas se defrontaram com o problema do
fechamento, ou seja, em algum ponto aparecem variáveis (quantidades estatísticas de alguma
ordem) que não possuem equações governantes. Para resolver esse problema diversas
aproximações baseadas em análise dimensional ou em medições experimentais foram
desenvolvidas. Porém muitas dessas tentativas, como fechamentos de um ou dois pontos, estão
limitados a escoamentos homogêneos (normalmente isotrópicos) onde simetrias permitem a
redução de variáveis.
15
2.3.1.
Grandes escalas do escoamento turbulento e estruturas coerentes
A definição de grandes e pequenas escalas num escoamento turbulento é de grande
utilidade tanto no meio experimental quanto na pesquisa em dinâmica dos fluidos
computacional.
As grandes escalas de um escoamento turbulento são aquelas que dominam a dinâmica do
mesmo. Elas são responsáveis pela maior parte do transporte de quantidade de movimento e pela
produção da turbulência [Ferziger, 1993]. A valores elevados do número de Reynolds essas
grandes escalas compreendem as estruturas da escala da geometria até aquelas na escala inercial
de Kolmogorov. Essas grandes estruturas são diretamente influenciadas pela geometria e
condições de contorno do escoamento. Dessa forma a obtenção de modelos genéricos para
descrição desses movimentos é complexa. Os modelos clássicos disponíveis são dependentes de
parâmetros empíricos ajustados para um ou alguns escoamentos específicos e sob condições bem
controladas. Esses valores podem ser muito diferentes do ideal para escoamentos em condições
diversas, o que prejudica a qualidade dos resultados obtidos.
As menores escalas tendem a ser mais homogêneas e isotrópicas e menos afetadas pelas
condições de contorno. Assim modelos mais universais e independentes em relação àqueles
utilizados na metodologia média clássica podem ser aplicados para sua descrição [Silva Freire et
al., 2002].
A descrição “grandes escalas” não está ligada unicamente ao comprimento característico
dos movimentos, em alguns casos importantes estruturas do escoamento são extremamente
pequenas, como ocorre em regiões próximas de contornos sólidos. Nesse tipo de condição, além
do tamanho reduzido dos vórtices, esses ainda apresentam elevada anisotropia o que dificulta
uma modelagem correta.
Como já citado o movimento turbulento apresenta estruturas organizadas, chamadas
estruturas coerentes. Em seu trabalho, Hussain, 1983, define uma estrutura coerente como sendo
uma massa de fluido turbulenta (de grandes escalas) conectada exibindo vorticidade
correlacionada em fase sobre sua extensão. Lesieur et al., 1995, complementa essa definição
afirmando que além de serem concentrações locais de vorticidade (ω
ωω
ω) essas estruturas também
devem ter forma reconhecível e devem apresentar um “tempo de vida” (T
c
) muito maior que o
tempo de giro local do vórtice. A definição de vorticidade e do tempo de vida estão nas equações
(2.25) e (2.26), respectivamente.
= ∇×
ω u
(2.25)
16
1
c
T
−
=
ω
(2.26)
Uma estrutura coerente também é caracterizada por apresentar (na região ocupada) uma
variação de pressão, usualmente uma queda.
As estruturas coerentes são responsáveis por parte significativa do transporte (de grandes
escalas) de massa, calor e momentum, sem necessariamente serem elas mesmas energéticas. O
que quer dizer que uma estrutura coerente mostra altos níveis de vorticidade, produção e
transporte de calor e massa coerentes, mas não necessariamente altos níveis de energia cinética
(a maior parte da energia cinética turbulenta está então associada à turbulência incoerente). Em
alguns escoamentos, entretanto, o transporte dessas grandezas devido à turbulência incoerente
pode se tornar comparativamente significante, como, por exemplo, em escoamentos plenamente
desenvolvidos [Hussain, 1983].
17
3.
SIMULAÇÃO NUMÉRICA DE ESCOAMENTOS TURBULENTOS
Uma alternativa aos métodos experimentais para análise de problemas envolvendo
escoamentos de fluidos é a dinâmica dos fluidos computacional (do inglês computational fluid
dynamics – CFD). A simulação numérica de escoamentos turbulentos está se difundindo cada
vez mais e se tornando uma ferramenta prática de engenharia. Códigos comerciais com
algoritmos cada vez mais eficientes e incorporando cada vez mais modelos para diferentes
fenômenos já se encontram a disposição.
Apesar da facilidade de se obter bons resultados para problemas bastante complexos, assim
como qualquer outra ferramenta em engenharia, deve-se tomar muito cuidado na hora de se usar
uma ferramenta tão poderosa.
Diversos métodos, tais como os métodos de diferenças finitas, volumes finitos, elementos
finitos e métodos espectrais, podem ser empregados. Cada um apresentando vantagens e
restrições próprias. Uma vez escolhido o método numérico resta ainda definir os modelos a
serem utilizados de acordo com as necessidades do problema em questão. Simulações numéricas
podem prover informações detalhadas sobre os escoamentos e são vantajosas quando muitas
variáveis dos escoamentos são necessárias ao mesmo tempo, ou onde trabalhos experimentais
são de difícil controle ou perigosas. Alguns casos, como por exemplo combustão, requerem o
conhecimento de detalhes finos do campo de escoamento, enquanto outros se satisfazem com
dados relativamente grosseiros. O grau de detalhamento desejado para uma dada análise está
diretamente ligado ao custo computacional necessário.
A simulação numérica da turbulência requer julgamentos adequados com respeito às
equações governantes, condições iniciais e de contorno e esquemas de solução utilizados nos
diferentes métodos. Um das maiores questões é a especificação de condições de contorno,
especialmente na região de entrada dos escoamentos e em contornos sólidos.
Nas regiões de entrada onde o escoamento incidente é turbulento as variáveis do
escoamento dependem do escoamento desconhecido fora do domínio. No caso de
indisponibilidade de condições bem conhecidas nesses contornos deve-se utilizar uma que
melhor represente o escoamento de entrada e ao mesmo tempo minimize a propagação de erros
de contorno. Valores equivocados nesse tipo de condição podem causar transição antecipada na
camada limite e também produzir aumento na turbulência, isso pode mover o ponto de separação
da camada limite resultante, alterando bastante os resultados da simulação. Em outro caso a
turbulência no escoamento de entrada pode causar um aumento nas forças flutuantes agindo
sobre um corpo imerso no escoamento [Ferziger, 1993].
18
Nas regiões próximas a contornos sólidos o principal problema é a resolução necessária
para captar corretamente os efeitos de camada limite. Devido aos altos gradientes das variáveis
do escoamento nessa região, e a grande influência de seus efeitos sobre os resultados uma má
representação da camada limite pode comprometer seriamente a validade dos campos obtidos.
Além das condições de contorno usuais nesse tipo de contorno, alguns autores utilizam
alternativas como condições tipo “lei da parede” no intuito de diminuir a requisição de malha
nessas regiões.
Existem três linhas principais para a modelagem de escoamentos turbulentos: a simulação
numérica direta (sem modelagem), a modelagem clássica (baseada em médias de Reynolds) e a
simulação de grandes escalas. Mais detalhes dessas linhas serão discutidos nas seções seguintes.
Apesar de ainda haver muito a ser feito na análise de escoamentos turbulentos muitos
resultados úteis para elaboração e comprovação de novas bases teóricas para o estudo da
turbulência já foram obtidos. Lesieur et al., 1995, em seu trabalho argumenta de forma otimista
sobre os avanços no estudo da turbulência dizendo: “Graças à dinâmica dos fluidos
computacional e ao grande progresso em técnicas experimentais é claro agora que a estrutura da
turbulência é muito mais simples do que se pensava”.
3.1.
Simulação numérica direta de escoamentos turbulentos
A simulação numérica direta (do inglês direct numerical simulation – DNS) consiste na
resolução direta de todas as escalas do escoamento turbulento. A principal vantagem desse tipo
de simulação é a eliminação da necessidade de modelos ad hoc e a obtenção de resultados
detalhados de todo o movimento turbulento [Rogallo e Moin, 1984].
Para que seja feita uma simulação numérica direta de um escoamento as equações de
conservação são resolvidas sem o uso de qualquer modelagem. Porém, para que isto seja
possível, é necessário que a resolução de malha seja capaz de captar até os menores movimentos
turbulentos (o que coloca limitações também na resolução no tempo da simulação), e que a
precisão dos métodos numéricos utilizados nessas simulações garantam que o erro introduzido
não seja significativo.
Quando devidamente conduzidas, os resultados obtidos através da simulação numérica
direta são comparáveis de todas as formas a dados experimentais de qualidade. Ainda com as
vantagens de que são obtidos dados detalhados de todas as variáveis de interesse no campo de
escoamento e as condições de contorno e iniciais são definidas de forma clara [Ferziger, 1993].
19
Apesar de todas essas vantagens, a aplicação da simulação numérica direta ainda é muito
limitada. O principal fator limitante é a resolução numérica necessária para a solução completa
de todas as escalas do escoamento turbulento. Tal resolução cresce ainda mais com o aumento do
número de Reynolds. Como já discutido, a relação entre as maiores escalas (delimitadas pela
geometria) e as menores aumenta com o aumento do número de Reynolds. O número de
variáveis necessárias para uma simulação numérica direta é da ordem de Re
9/4
em três
dimensões, o que limita seu uso a casos com baixos números de Reynolds. Outro problema é a
precisão do método numérico, o que força o uso de esquemas numericamente eficientes os quais
normalmente conduzem a limitações geométricas.
As principais aplicações da simulação numérica direta atualmente são o estudo do
fenômeno da transição e da formação de estruturas turbulentas e a geração de uma massa de
dados sobre casos de referência para uso na validação de novos métodos, numéricos e
experimentais.
3.2.
Modelagem clássica de escoamentos turbulentos
Devido aos problemas já discutidos relativos à simulação de escoamentos turbulentos, uma
metodologia foi desenvolvida baseada na idéia de O. Reynolds de que as variáveis de campo em
um escoamento turbulento (como velocidade e pressão) podem ser decompostas em uma parte
média e uma parte flutuante.
i
i i
u u u
p p p
′
= +
′
= +
(3.1)
Onde
i
u
e
p
são as parcelas médias da componente
i
da velocidade e da pressão e
i
u
′
, e
p
′
as parcelas flutuantes das mesmas variáveis. Para a obtenção do valor médio é utilizado o
seguinte processo de integração no tempo. Para uma variável genérica
g
o valor médio é obtido
através da integral da equação (3.2) e o valor flutuante pela expressão (3.3), respectivamente.
0
0
1
( )
t t
t t
g g t dt
t
+∆
−∆
=
∆
∫
(3.2)
g g g
′
= −
(3.3)
20
onde
∆
t
é um tempo suficientemente longo para que todas as escalas significativas sejam
consideradas e
t
0
um determinado instante de tempo na integração. O operador de média possui
as seguintes propriedades importantes.
0
i i
g g
g
g g
x x
=
′
=
∂ ∂
=
∂ ∂
(3.4)
Com essa integração no tempo, apenas fenômenos com escalas de tempo maiores do que a
dos maiores vórtices turbulentos poderão ser consideradas (como variações no tempo das
condições de contorno). Aplicando a definição das variáveis dadas pela relação (3.1) nas
equações (2.15) e (2.16) obtém-se as seguintes expressões.
0
j
j
u
x
∂
=
∂
(3.5)
1
i i i
j
i j
i
j i j j
u u p u
u u u f
t x x x x
ν
ρ
∂ ∂ ∂ ∂ ∂
′ ′
+ = − + − +
∂ ∂ ∂ ∂ ∂
(3.6)
Observa-se que a equação da continuidade permanece inalterada, indicando que o campo
médio deve respeitar a conservação de massa, já na equação de conservação de quantidade de
movimento aparece um termo adicional, que representa o efeito das flutuações sobre o campo
médio de escoamento.
Essas equações são conhecidas como as equações médias de Reynolds (do inglês Reynolds
averaged Navier-Stokes equations – RANS) que dão nome a essa metodologia de aproximação,
hoje conhecida também como modelagem clássica da turbulência.
Todos os termos das equações (3.5) e (3.6) representam o movimento do escoamento
médio, exceto o termo que envolve a média do produto das flutuações conhecido como tensor de
Reynolds, ou tensor de tensões turbulentas. O aparecimento desse termo leva ao problema de
fechamento, que é a ausência de equações de conservação bem conhecidas para representar essas
flutuações. Dessa forma aparecem mais incógnitas que equações no sistema a ser resolvido,
modelos devem então ser utilizados para representar essas variáveis adicionais em função das
variáveis médias.
21
A fim de modelar-se os efeitos do tensor de Reynolds sobre o escoamento médio a ser
resolvido, duas principais linhas de aproximação foram desenvolvidas. Uma segue a hipótese de
Boussinesq (discutida em detalhe adiante) e a outra trabalha com equações de transporte do
próprio tensor de Reynolds.
Por resolver apenas as variáveis médias do escoamento, a modelagem clássica da
turbulência permite o uso de diversas simplificações para os problemas em estudo.
Simplificações como simetrias, escoamentos bidimensionais e em regime permanente são
aproximações freqüentemente utilizadas com modelos clássicos.
Também por essas características, os modelos de turbulência também ditos convencionais,
devem levar em consideração todas as escalas do movimento turbulento. Isso é a grande
dificuldade na elaboração desses modelos, principalmente pelo fato de que as grandes escalas da
turbulência são enormemente influenciadas pela geometria e pelas condições de contorno, o que
dificulta a obtenção de um modelo universal. Muitos modelos utilizados possuem parâmetros
que são ajustados empiricamente, porém para alguns escoamentos de referência como
turbulência isotrópica ou escoamentos turbulentos em canais. Em algumas situações, bastante
diferente daquelas para os quais foram ajustados, os modelos fornecem resultados bastante
imprecisos. Um dos exemplos mais conhecidos é a dificuldade do modelo k-ε em escoamentos
com separação de camada limite, o que prejudica, por exemplo, a previsão da sustentação ou do
arrasto de um perfil aerodinâmico a partir de um campo resolvido com esse modelo.
Outra grande dificuldade no uso desses modelos é de se obter resultados práticos quando as
propriedades de interesse na simulação são justamente aquelas ligadas às flutuações no campo de
escoamento. O que ocorre no estudo de vibrações induzidas pelo escoamento em um banco de
tubos, por exemplo.
3.3.
Simulação de grandes escalas
A simulação de grandes escalas - SGE (do inglês large-eddy simulation – LES) é uma
alternativa intermediária entre a simulação numérica direta e a modelagem clássica da
turbulência por médias de Reynolds. Ela foi introduzida pelo meteorologista Smagorinsky, 1963,
e desde então vem sendo estudada e aprimorada. Com a possibilidade de se tornar, num futuro
próximo, uma ferramenta prática de engenharia.
A simulação de grandes escalas consiste em resolver diretamente as maiores estruturas
presentes no escoamento (as maiores escalas) e modelar o efeito das menores estruturas
(menores escalas ou escalas sub-malha) naquelas resolvidas. Isso é feito através de uma
22
operação de filtragem (ou média espacial) aplicada sobre as equações de conservação (de massa
e quantidade de movimento) e uso de um modelo para os termos necessários para fechamento.
Esta metodologia permite o uso de malhas menos refinadas do que aquelas necessárias para
simulação numérica direta, sendo então mais viável.
A idéia da simulação de grandes escalas é de aproveitar as características das diferentes
escalas da turbulência já discutidas anteriormente. A menores escalas sendo mais universais,
permitem utilizar modelos de turbulência mais simples do que aqueles utilizados na modelagem
clássica para simular seus efeitos. As grandes escalas do escoamento, que são dependentes da
geometria e das condições de contorno, portanto peculiares a um problema específico, são
diretamente resolvidas. Lembrando-se que em muitas situações é o movimento das grandes
escalas que domina a dinâmica e dos escoamentos, esse tipo de simulação é bem mais realista
que a modelagem clássica, possibilitando a detecção de fenômenos que são “mascarados” na
modelagem por médias de Reynolds.
Boas revisões sobre a metodologia da simulação de grandes escalas podem ser encontradas
nos trabalhos de Rogallo e Moin, 1984, Lesieur et al., 1995 e Piomelli, 1999. Em linhas gerais
um bom texto (um pouco mais resumido) também é encontrado na referência Silva Freire et al.,
2002. Outros trabalhos citados ao longo deste ressaltam diferentes aspectos e cuidados especiais
a serem tomados na aplicação da simulação de grandes escalas.
A separação entre grandes escalar (resolvidas) e pequenas escalas (modeladas) é feita pelo
operador de filtragem; os movimentos de ordem superior ao comprimento característico do filtro
são as grandes escalas, e os de ordem inferior são as menores. Essa definição não é
necessariamente a definição utilizada no estudo teórico da turbulência, onde para altos valores do
número de Reynolds as maiores escalas são aquelas da ordem da geometria até a ordem da escala
inercial de Kolmogorov, e as menores são desta última até as escalas dissipativas. Por isso alguns
autores definem dois tipos de simulação de grandes escalas: a simulação de grandes escalas
física, onde a divisão de escalas é a mesma aplicada ao estudo teórico da turbulência; e a
simulação de grandes escalas numérica, onde essa condição não é cumprida [Pope, 2003]. Uma
classificação alternativa é apresentada por Ferziger, 1993. Ele classifica como captura de
estruturas coerentes (do inglês coherent structure capture – CSC) uma simulação nas bases da
simulação de grandes escalas, porém com uma resolução numérica suficiente apenas para
capturar os efeitos dos maiores vórtices. Com isso em casos onde os valores de interesse de um
dado escoamento são dependentes apenas das maiores escalas e das maiores flutuações de
pressão, bons resultados são obtidos com um baixo custo computacional. Entretanto, a captura de
23
estruturas coerentes falha em casos onde efeitos de arraste viscoso ou de estruturas menores que
a malha são predominantes.
As menores escalas do escoamento, que não são resolvidas na simulação de grandes
escalas, são também chamadas de escalas sub-malha (do inglês subgrid-scales – SGS). Seus
efeitos sobre o desenvolvimento do campo de escoamento resolvido são modelados através de
modelos especiais de fechamento [Findikakis e Street, 1982].
Como a simulação de grandes escalas tem o objetivo de simular grandes estruturas da
turbulência (como anteriormente apresentadas, tridimensionais) ela é tridimensional por
definição, alguma descrição bidimensional de um escoamento feita utilizando essa metodologia
só pode ser considerada em situações em que evidências físicas validam a hipótese de que as
grandes estruturas são predominantemente bidimensionais. Nesses casos as variáveis obtidas
representariam médias das variáveis reais na terceira dimensão [Findikakis e Street, 1982].
A qualidade de um resultado obtido numa simulação de grandes escalas é dependente de
diversos fatores como o método numérico, a resolução da discretização utilizada, o tamanho
característico do filtro, o modelo utilizado para as menores escalas e as condições de contorno
aplicadas. Cada um desses fatores tem uma influência sobre a resposta e diferentes combinações
destes também influenciam qualitativa e quantitativamente os resultados. Discussões sobre cada
alternativa serão feitas nas seções subseqüentes.
3.3.1.
Filtragem espacial para simulação de grandes escalas
Para a separação de escalas (em resolvidas e modeladas) em simulação de grandes escalas
de escoamentos turbulentos é aplicado um filtro sobre as equações de conservação [Silva Freire
et al., 2002]. As variáveis do problema são divididas em uma componente de grandes escalas (a
ser resolvida) e uma componente de pequenas escalas (chamadas escalas sub-malha, a serem
modeladas), assim:
i i i
u u u
p p p
′
= +
′
= +
(3.7)
Onde
i
u
e
p
são os valores filtrados (também chamados de valores médios espaciais) da
componente i do vetor velocidade e da pressão e
i
u
′
e
p
′
são os valores de flutuação da
componente
i
do vetor velocidade e da pressão.
24
A operação de filtragem é definida como uma convolução entre a variável de campo e a
função filtro [Leonard, 1974]. Os valores filtrados são obtidos a partir do operador filtro
mostrado na equação (3.8).
(
)
(
)
(
)
, ,
g t G g t d
′ ′ ′
= −
∫
x x x x x
(3.8)
Onde
g
é uma variável genérica,
g
seu valor filtrado e G a função filtro, que pode assumir
diversas formas, sendo mais comum o uso do filtro “caixa”. O filtro “caixa” é uma função que
apresenta um valor constante dentro do raio característico do filtro, como dado pela equação
(3.9) e mostrado na Figura 3.1a. Outros filtros freqüentemente utilizados são o filtro gaussiano,
definido pela equação (3.10) e mostrado também na Figura 3.1b, e o filtro passa alta (filtro de
corte) no espaço de Fourier (em números de onda) dado pela equação (3.11) [Leonard, 1974;
Silveira Neto et al., 1993; Piomelli, 1999] . Outros filtros são descritos nos artigos de Leonard,
1974, e Findikakis e Street, 1982.
( )
3
1
se
2
0 se
2
G
∆
′
− ≤
∆
′
− =
∆
′
− >
x x
x x
x x
(3.9)
Onde
∆
é a dimensão característica do filtro.
( )
2
2
3
9
2
6
G e
π
′
− −
∆
′
− =
∆
x x
x x (3.10)
1 se
( )
0 se
k
G k
k
π
π
≤
∆
=
>
∆
(3.11)
Onde
G
é a função filtro e k é o número de onda definidos no espaço de Fourier.
25
Figura 3.1. a) Forma da função filtro tipo “caixa”; b) Forma da função filtro gaussiana.
A operação de filtragem possui as seguintes propriedades:
0
g f g f
g f gf
g f
≠
′ ′
≠
′
≠
(3.12)
Essas propriedades a diferenciam da operação de média temporal da decomposição de
Reynolds. Além disso, se um filtro é uniforme (com comprimento característico
∆
constante,
dependente apenas de |
x
-
x
′|) é empregado, as operações de filtragem e derivação são
comutativas, ou seja:
i i
g g
x x
g g
t t
∂ ∂
=
∂ ∂
∂ ∂
=
∂ ∂
(3.13)
Quando um filtro não-uniforme é utilizado aparecem erros de comutação nas derivadas
espaciais. Piomelli, 1999, cita trabalhos que demonstram que esse erro é da ordem do quadrado
do tamanho característico do filtro, O(
∆
2
), e que quando são utilizados esquemas numéricos de
segunda ordem na discretização do problema o erro passa a ser da mesma ordem do erro de
truncamento. O tamanho característico do filtro não é necessariamente igual ao tamanho
característico da malha (h). Em alguns casos o tamanho característico do filtro é definido como a
raiz cúbica do produto dos comprimentos característicos da malha nas três dimensões, em outros,
G
3
1
∆
2
∆
−
2
∆
|x
-
x
′
′′
′
|
a)
G
|x
-
x
′
′′
′
|
b)
26
como em malhas bastante anisotrópicas, utiliza-se uma escala proporcional à norma euclidiana
da dimensão dos elementos [Findikakis e Street, 1982].
Uma das diferenças entre o filtro tipo caixa e o filtro gaussiano é que o filtro tipo caixa
possui contribuição de todas as escalas sub-malha, “absorvendo” seus movimentos. Por esse
motivo, não é possível obter um espectro completo da turbulência com o uso desse tipo de
função filtro. Para evitar esse tipo de problema muitos autores utilizam o filtro gaussiano, que
não elimina totalmente os movimentos das menores escalas [Findikakis e Street, 1982].
Aplicando-se a filtragem nas equações de conservação na forma incompressível, dadas
pelas equações (2.15) e (2.16), obtém-se as equações de conservação para as variáveis filtradas
(para simulação de grandes escalas) na forma incompressível, escritas como as equações (3.14)
para a conservação de massa filtrada e (3.15) para a conservação da quantidade de movimento
filtrada.
0
j
j
u
x
∂
=
∂
(3.14)
1
i i j
j i
i
j i j j i
u u
u p u u
f
t x x x x x
ν
ρ
∂
∂ ∂ ∂ ∂ ∂
+ = − + + +
∂ ∂ ∂ ∂ ∂ ∂
(3.15)
A equação da continuidade permanece inalterada em sua forma, passando apenas a
trabalhar com a velocidade filtrada. O campo filtrado de velocidades então deve satisfazer
também a conservação de massa. Já na equação de Navier-Stokes, o termo advectivo aparece
com o produto filtrado das velocidades, o que não pode ser trabalhado de forma direta.
Inserindo a definição de velocidade filtrada e flutuante dada pela equação (3.7) no produto
filtrado das velocidades chega-se à:
(
)
(
)
j i j i j i
j i j i i j j i
u u u u u u u u u u u u u u
′ ′ ′ ′ ′ ′
= + + = + + +
(3.16)
Agrupando-se os termos e inserindo os termos de Leonard chega-se à equação (3.17).
j i
j i ij ij ij
u u u u C L
τ
= + + +
(3.17)
27
onde
τ
ij
representa as tensões das escalas sub-malha (produto filtrado das flutuações de
velocidade), C
ij
os termos cruzados (com contribuições das variáveis filtradas e das variáveis
sub-malha) e L
ij
os termos de Leonard, e que são dados pelas expressões (3.18), (3.19) e (3.20)
respectivamente.
ij j i
u u
τ
′ ′
= (3.18)
j i
ij i j
C u u u u
′ ′
= + (3.19)
j i j i
ij
L u u u u
= − (3.20)
Inserindo a equação (3.17) na equação (3.15) chega-se à forma filtrada completa da
equação de Navier-Stokes para a forma incompressível.
1
i j i i j
ij ij ij
i
j i j j i
u u u p u u
C L f
t x x x x x
ν τ
ρ
∂ ∂ ∂ ∂ ∂ ∂
+ = − + + − − − +
∂ ∂ ∂ ∂ ∂ ∂
(3.21)
Os termos de tensões sub-malha e cruzados representam os efeitos das escalas sub-malha
sobre as escalas resolvidas e a interação entre as duas escalas. O termo de Leonard representam
exatamente a diferença entre o valor do produto filtrado das velocidades filtradas e o produto
simples das velocidades filtradas ,como mostrado na equação (3.20). Esse termo foi definido por
Leonard, 1974, para melhor representar a cascata de energia no escoamento turbulento, fazendo
com que os efeitos das menores escalas não fosse capturado apenas pelo modelo sub-malha
utilizado, diminuindo a dependência sobre esse modelo.
Os termos de Leonard e termos cruzados são aproximados por uma expansão em série de
Taylor a partir do centro do volume filtrado, obtendo-se a seguinte expressão [Findikakis e
Street, 1982; Silveira Neto et al., 1993]:
2
2
k
i j
ij ij
k k
u u
C L
x x
γ
∆ ∂ ∂
+
∂ ∂
≃
(3.22)
28
Onde
∆
k
é o tamanho característico do filtro na direção k e
γ
é um valor constante definido
pela função filtro utilizada; para um filtro tipo caixa, por exemplo, o valor de
γ
é 6. Leonard,
1974, apresentou em seu trabalho uma análise da contribuição desses termos na representação da
cascata de energia, chegando à conclusão de que eles não podiam ser desprezados. Porém em
análises posteriores e mais detalhadas, diversos autores, como Silveira Neto et al., 1993, e Petry
e Awruch, 1997, concluíram que na verdade esses termos apresentam valores desprezíveis.
No trabalho de Rogallo e Moin, 1984, é dito que a forma correta de medir-se a importância
dos termos de Leonard é medir-se a transferência de energia associada ao termo advectivo em
um campo previamente resolvido (como em uma simulação direta da turbulência). É citado
também que alguns estudos nesse sentido já haviam sido realizados e que existia sim
transferência de energia, porém em uma ordem de grandeza muito inferior àquela prevista por
Leonard, 1974.
Na prática, quando se utiliza um filtro tipo caixa (no domínio físico) ou um filtro passa de
corte (no domínio de Fourier) e esquemas de discretização de baixa ordem, a contribuição desses
termos é da ordem do erro numérico [Piomelli, 1999]. O que faz com que muitos autores nem
escrevam tais termos, ou considere seus efeitos como inseridos no termo de tensões sub-malha
[Lesieur et al., 1995; Piomelli, 1999].
O termo de tensões sub-malha sim recebe atenção especial. Ele representa diretamente a
contribuição das menores escalas para o escoamento resolvido, e não pode ser considerado
desprezível. A aproximação usual é feita através da hipótese de Boussinesq, mostrada na
equação (3.23).
2
3
ij
ij
ij kk T
S
δ
τ τ ν
− −
≃
(3.23)
Onde
δ
ij
é a função delta de Kroneker, ν
T
é a viscosidade turbulenta associada localmente
ao escoamento e
ij
S
é o tensor taxa de deformação de grandes escalas, calculado pela equação
(3.24). Essa é a primeira aproximação feita para contabilizar a influência das menores escalas no
escoamento médio como função do movimento das grandes escalas. O termo envolvendo as
tensões normais é relacionado com a energia cinética turbulenta, como mostrado pela equação
(3.25).
1
2
i j
ij
j i
u u
S
x x
∂ ∂
= +
∂ ∂
(3.24)
29
2
kk k k
u u k
τ
′ ′ ′
= = (3.25)
Onde
k
′
é a energia cinética turbulenta filtrada. O termo surge para considerar-se apenas a
parte anisotrópica do tensor de tensões sub-malha, inserindo a parte isotrópica junto com a
pressão estática do fluido (por ser uma tensão normal flutuante). O termo é necessário também
por se estar trabalhando (até este ponto do presente trabalho) com as equações de conservação na
forma incompressível. O que, na ausência desse termo adicional, geraria valores nulos de tensão
sub-malha quando i=j [Hinze, 1975].
Inserindo a aproximação (3.23) na equação (3.21) e, conforme justificado, desprezando os
termos cruzados e de Leonard, chega-se à:
( )
1
i j i i j
T
i
j i j j i
u u u P u u
f
t x x x x x
ν ν
ρ
∂ ∂ ∂ ∂ ∂ ∂
+ = − + + + +
∂ ∂ ∂ ∂ ∂ ∂
(3.26)
Onde a pressão foi redefinida como:
3
kk
P p
ρ
τ
≡ + (3.27)
Onde
P
é a pressão filtrada redefinida na equação de Navier-Stokes [Lesieur et al., 1995].
A única variável sem alguma equação de conservação ou constitutiva é a viscosidade
turbulenta, para ela é necessário algum tipo de equação de aproximação para fechamento do
problema. Apesar de representar os efeitos das escalas sub-malha no escoamento médio, uma
aproximação deve ser feita em função das escalas resolvidas do problema. Esses são os
diferentes modelos sub-malha utilizados na simulação de grandes escalas.
A aplicação da filtragem nas equações de Navier-Stokes fornece equações que são
semelhantes àquelas utilizadas na modelagem por médias de Reynolds. Porém a modelagem
clássica envolve um campo médio estacionário (ou com pequena variação no tempo) que varia
suavemente no espaço. Já a simulação de grandes escalas envolve um campo que é bastante
caótico no espaço e no tempo quando a malha é suficientemente pequena. Se o tamanho
característico da malha aumenta muito (e com ele o tamanho característico do filtro), o resultado
obtido com a simulação de grandes escalas de certa forma se aproxima daquele obtido usando-se
a modelagem clássica [Lesieur et al, 1995].
30
Trabalhando-se com as equações de conservação com hipótese de quase-
incompressibilidade, como apresentadas nas equações (2.13) e (2.14) obtém-se uma expressão
um pouco diferente. Principalmente pelo fato de que não é necessária a redefinição da pressão, já
que a parte relativa à energia cinética turbulenta na hipótese de Boussinesq não é necessária. A
forma da aproximação de Boussinesq utilizada para o termo de tensões sub-malha é vista na
equação (3.28) [Petry, 2002].
2
ij
ij T
S
τ ν
−
≃
(3.28)
Com a mesma definição do tensor taxa de cisalhamento mostrada anteriormente, efetuando
a operação de filtragem espacial nas equações de conservação para quase-incompressibilidade,
inserindo a essa aproximação para o tensor de tensões sub-malha e desprezando os termos
cruzados e de Leonard, chega-se à forma quase-incompressível das equações de conservação
para simulação de grandes escalas.
2
0
j
j
p u
c
t x
ρ
∂ ∂
+ =
∂ ∂
(3.29)
( )
1
i j i j i k
T
i
j i j i j i k
u u u p u u u
f
t x x x x x x x
λ
ν ν
ρ ρ
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
+ = − + + + + +
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
(3.30)
As equações (3.29) e (3.30) representam então as equações de conservação de massa e
quantidade de movimento, respectivamente, para o campo filtrado (de grandes escalas) de um
escoamento turbulento, tridimensional, transiente e isotérmico de um fluido newtoniano com
hipótese de quase-compressibilidade. Essas equações serão utilizadas no restante do trabalho.
Pode-se definir ainda uma nova viscosidade, chamada viscosidade efetiva, que é a soma da
viscosidade molecular do fluido com o valor de viscosidade turbulenta obtido para o
escoamento, essa viscosidade efetiva (
ν
ef
) é dada então por:
ef T
ν ν ν
= +
(3.31)
31
Na forma adimensional, obtida tanto pela adimensionalização das equações (3.29) e (3.30)
ou pela aplicação da operação de filtragem e aproximações nas equações (2.18) e (2.19), é
mostrada nas equações (3.32) e (3.33).
* *
* *
1
0
j
j
p u
t Mc x
∂ ∂
+ =
∂ ∂
(3.32)
* * * * * * *
*
* * * * * * * *
1 1 1 1
i j i j i k
i
j i j L T i j i k
u u u p u u u
f
t x Mc x x Re Re x x Re x x
λ
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
+ = − + + + + +
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
(3.33)
Onde Re
T
é o número de Reynolds turbulento, definido em função da viscosidade
turbulenta.
UL
T
T
Re
ν
= (3.34)
3.3.2.
Modelos sub-malha
O modelo sub-malha é utilizado para representar a influência das escalas não-resolvidas
(movimentos sub-malha) sobre as escalas resolvidas. Em modelos baseados no conceito de
viscosidade turbulenta (pela hipótese de Boussinesq), a função do modelo sub-malha é calcular
um valor para a viscosidade turbulenta que, quando inserido na equação do movimento filtrada,
cause uma dissipação de energia que simule a transferência de energia para as pequenas escalas
(a cascata de energia), e em casos especiais o efeito contrário, o fenômeno de backscatter.
O papel do modelo então não é de fornecer essas estatísticas para as grandes escalas, mas
prevenir que a omissão das pequenas escalas (modeladas) comprometa o cálculo do campo de
grandes escalas sobre o qual as estatísticas são tomadas.
A elaboração de um modelo sub-malha ideal ainda é um desafio para os estudiosos de
turbulência e, principalmente, para os numericistas que desenvolvem a simulação de grandes
escalas. As características desejadas para tal modelo é que ele seja capaz de prever corretamente
a dissipação geral de energia no escoamento, que seja capaz de se anular em escoamentos
laminares, que dependa fortemente das menores escalas resolvidas (mais do que do espectro
completo de movimentos) e que preveja precisamente a troca local de energia entre as escalas
resolvidas e as modeladas [Piomelli, 1999].
32
Nas sub-seções a seguir serão discutidos alguns dos principais modelos sub-malha.
a.
Modelo de Smagorinsky
O primeiro modelo sub-malha foi proposto por Smagorinsky, 1963, em seu trabalho
pioneiro, para simulação de grandes escalas de escoamentos atmosféricos quase-bidimensionais.
A partir da hipótese de equilíbrio para as menores escalas, que diz que as menores escalas
dissipam instantaneamente toda a energia que é transferida para elas, ele deduziu uma expressão
para a viscosidade turbulenta mostrada na equação (3.35).
(
)
2
T S
C S
ν
= ∆
(3.35)
Onde C
S
é a constante de Smagorinsky e
S
é o módulo do tensor taxa de cisalhamento
calculado pela equação (3.36).
2
ij ij
S S S
=
(3.36)
A constante de Smagorinsky foi calculada para turbulência isotrópica apresentando valores
entre 0,18 e 0,23. Em escoamentos com cisalhamento médio, próximo a contornos sólidos e em
escoamentos em transição esse valor deve ser reduzido, isso se deve ao fato de que o modelo de
Smagorinsky ser muito dissipativo.
Em escoamentos em transição, essa dissipação aumentada faz com que as primeiras
perturbações no escoamento sejam amortecidas e não causem a transição. Também em
escoamentos em transição e em escoamentos próximos a superfícies sólidas (escoamentos que
não satisfazem a condição de equilíbrio de energia para as menores escalas), ocorre o fenômeno
da transferência inversa de energia (o backscatter), fenômeno esse que não é captado pelo
modelo de Smagorinsky, que só permite transferência de energia das maiores para as menores
escalas. A capacidade de captar esse fenômeno é imprescindível para garantir-se a qualidade dos
resultados obtidos. Um estudo interessante sobre o fenômeno de transferência inversa de energia
no escoamento em um canal é apresentado em Piomelli et al., 1991.
Outra dificuldade do modelo é que ele não é capaz de anular as tensões turbulentas sub-
malha quando o escoamento é laminar [Findikakis e Street, 1982; Piomelli et al., 1993; Ferziger,
1993; Lesieur et al., 1995; Piomelli, 1999].
33
Para superar essas dificuldades muitos autores propuseram modificações no modelo de
Smagorinsky; funções de intermitência para “desligar” o modelo onde necessário e funções de
amortecimento tipo van Driest próximo a paredes são exemplos dessas modificações [Piomelli et
al., 1991; Ferziger, 1993]. Em Piomelli et al. 1988, são citados trabalhos que propuseram uma
decomposição do modelo de Smagorinsky em duas partes, uma para modelar os vórtices
isotrópicos e outra para modelar as não-homogeneidades locais nas paredes.
Apesar dessas dificuldades, muitos trabalhos foram realizados utilizando o modelo de
Smagorinsky com resultados de qualidade, trabalhando com escoamentos adequados às suas
limitações. Uma crítica feita sobre o modelo de Smagorinsky é de que os bons resultados obtidos
com sua aplicação refletem sua habilidade de estabilizar os cálculos realizados e não de
reproduzir com precisão os efeitos das escalas sub-malha, o que indicaria também que as maiores
escalas dos escoamentos são insensíveis aos detalhes dos movimentos sub-malha naqueles
escoamentos onde são obtidos resultados satisfatórios [Rogallo e Moin, 1984].
Um detalhe sobre o uso do modelo de Smagorinsky é a combinação dele com alguns tipos
de filtro, como mostrado no trabalho de Piomelli, 1999, a aplicação do modelo de Smagorinsky
só é adequada quando se tem o tamanho das escalas sub-malha bem definidos (ou seja, quando
se tem claramente a escala de corte das menores estruturas) como acontece quando se utiliza uma
função filtro tipo caixa. O modelo de Smagorinsky impede que se contabilize o efeito de todas as
escalas sobre o campo filtrado, que é uma das características do filtro gaussiano. Dessa forma a
combinação do modelo de Smagorinsky com um filtro gaussiano leva a resultados ruins, sendo
mais apropriada sua combinação com o filtro tipo caixa.
b.
Modelo dinâmico
Um dos modelos sub-malha que fornece os melhores resultados para escoamentos sob
diversas condições é o modelo dinâmico. O modelo dinâmico pode ser visto também como uma
estrutura de modelagem, não só como um modelo sub-malha específico, já que ele pode ser
aplicado sobre diversos outros modelos sub-malha. O modelo dinâmico baseia-se a idéia de que
as constantes do modelo sub-malha sobre o qual é utilizado podem ser avaliadas localmente no
escoamento a partir das menores escalas resolvidas. Desta forma as “constantes” se adaptam as
características locais (no espaço e no tempo) do escoamento, simulando com mais precisão os
efeitos das escalas modeladas.
O modelo dinâmico foi apresentado no trabalho de Germano et al., 1991, e mais tarde
Lilly, 1992, modificou a forma apresentada para o cálculo do coeficiente dinâmico. A proposta é
34
de se aplicar uma segunda filtragem sobre o escoamento resolvido e relacionar as tensões
turbulentas das menores escalas resolvidas como as tensões das escalas sub-malha.
A aplicação mais comum da modelagem dinâmica é sobre o modelo de Smagorinsky (para
o qual foi derivado). Nesse caso a constante de Smagorinsky C
s
é substituída por um valor
dinâmico C. Para a obtenção de uma expressão para o cálculo de C aplica-se uma segunda
filtragem (filtro teste) sobre a equação da conservação da quantidade de movimento na forma
(3.15). O tamanho característico do filtro teste
∆
é sempre tomado maior que o tamanho
característico da filtragem de grandes escalas, geralmente o dobro. A nova equação de
conservação com a segunda filtragem é escrita como:
ɶ
ɶ
ɶ
1
i i j
j i
i
j i j j i
u u
u p u u
f
t x x x x x
ν
ρ
∂
∂ ∂ ∂ ∂ ∂
+ = − + + +
∂ ∂ ∂ ∂ ∂ ∂
(3.37)
Define-se então um novo tensor de tensões sub-malha (no nível do filtro teste) como
mostrado na equação (3.38).
ɶ
ɶ
i j
ij i j
T u u u u
= −
(3.38)
Onde T
ij
é o tensor de tensões sub-malha. Define-se então um outro tensor que representa
as tensões causadas pelas menores escalas resolvidas da malha, essas escalas estão em uma faixa
de comprimento denominada “janela de teste”.
ɶ
ɶ
i j i j ij
ij ij
L u u u u T
τ
= − = −
ɶ
(3.39)
Aqui L
ij
é o tensor de tensões turbulentas das menores escalas resolvidas, escalas de
comprimento maior do que o primeiro filtro
∆
e menor do que o filtro teste
∆
, e
τ
ɶ
é o tensor
filtrado no nível teste de tensões sub-malha no nível da primeira filtragem.
Supondo que ambos os tensores de tensões turbulentas possam ser modelados pelo mesmo
funcional, o modelo de Smagorinsky, na forma da equação (3.35), chega-se as equações (3.40),
(3.41) e (3.42).
2
2
3
ij
ij
ij kk
C S S
δ
τ τ
− − ∆
≃
(3.40)
35
2
2
3
ij
ij
ij kk
T T C S S
δ
− − ∆
≃
(3.41)
2
3
ij
ij kk ij
L L CM
δ
− −
≃
(3.42)
com o tensor M
ij
calculado a partir das equações (3.39), (3.40) e (3.41) como:
2 2
2 2
ij ij
ij
M C S S C S S
= ∆ − ∆ (3.43)
Com a definição da equação (3.43), pode-se calcular o valor de C a partir da equação
(3.42). O problema que surge é que a equação (3.42) não fornece um único valor para C; isso foi
resolvido pela aplicação de uma aproximação de mínimos quadráticos para o erro total da
equação
2
2
3
ij
ij kk ij
Q L L CM
δ
= − +
(3.44)
onde Q é o erro quadrático da equação (3.42). O valor de C para o mínimo erro é obtido pela
operação:
0
Q
C
∂
=
∂
(3.45)
chegando-se a:
2
1
2
ij ij
ij
L M
C
M
= −
(3.46)
A viscosidade turbulenta, utilizada na equação de Navier-Stokes filtrada para simulação de
grandes escalas (equação (3.26), (3.30) ou (3.33)), é calculada então a partir da fórmula:
2
T
C S
ν
= ∆
(3.47)
36
com o valor de
C
obtido a partir da equação (3.46).
Utilizando esse modelo o valor da viscosidade turbulenta se anula para escoamentos
laminares e próximos a contornos sólidos (já que
L
ij
se anula nesses casos), o que é o
comportamento correto nesses casos. Além disso o modelo é capaz de reproduzir o fenômeno da
transferência inversa de energia (
backscatter
) o que o torna adequado a simular, por exemplo,
escoamentos em transição.
O único parâmetro definido fora do processo de cálculo é a razão entre o comprimento
característico do filtro teste (
∆
) e o comprimento característico do filtro para simulação de
grandes escalas (
∆
). O que se tem definido é que esse valor tem que ser maior que a unidade.
Testes demonstraram que o modelo não é muito sensível a esse número, mas um valor ótimo
conhecido é em torno de 2.
Um inconveniente do modelo dinâmico é o aparecimento de variações abruptas no espaço
e no tempo dos valores do coeficiente dinâmico (e portanto da viscosidade turbulenta). Isso pode
ser resolvido através de médias (no espaço ou no tempo) do valor do coeficiente. Em seu
trabalho, Germano et al., 1991, fez médias sobre todas as direções homogêneas dos casos
testados, resolvendo o problema. Uma média na escala de comprimento do filtro teste é
justificada pelo fato de que o mesmo coeficiente deve ser aplicado para as tensões turbulentas
sub-malha e para aquelas dentro do filtro teste [Piomelli, 1999].
Uma proposta de algoritmo para cálculo do coeficiente dinâmico e para a média espacial no nível
do filtro teste é apresentada em Petry e Awruch, 2004, onde o processo de cálculo e média
espacial é baseado no método de elementos finitos. No mesmo trabalho é empregado um
limitador para valores negativos de viscosidade turbulenta, sendo este limitador também usado
por Zang et al, 1993. O limitador faz com que a soma das viscosidades molecular e turbulenta (a
viscosidade efetiva) seja sempre maior ou igual a zero, como mostra a equação (3.48).
0
ef T
ν ν ν
= + ≥
(3.48)
Isso também evita que o processo numérico se torne instável, além do fato de que um valor
de viscosidade efetiva negativo não teria um sentido físico aceitável.
Como o modelo dinâmico foi aplicado sobre a base do modelo de Smagorinsky, deve-se
utilizar uma função filtro tipo caixa na filtragem das equações de conservação, sob pena de perda
de qualidade nos resultados se for utilizado um filtro tipo gaussiano, por exemplo.
37
Outros modelos sub-malha, como modelos mistos ou tipo função de estrutura também são
bastante utilizados. Detalhes sobre tais modelos podem ser encontrados nas revisões escritas por
Lesieur et al., 1995, e Piomelli, 1999.
38
4. MODELAGEM MATEMÁTICA E NUMÉRICA
O método numérico utilizado no trabalho para a aproximação das equações de conservação
é o método de elementos finitos. Inicialmente deduzido para trabalhar com problemas
estruturais, o método de elementos finitos tornou-se uma ferramenta completa para lidar com
diversos tipos de problemas de engenharia. A bibliografia sobre o método é bastante extensa e
detalhada. Alguns trabalhos indicados para estudo do método em geral e sua aplicação no
problema de dinâmica de fluidos computacional são aqueles de Hughes, 1987, Zienkiewicz e
Taylor, 2000, e Reddy e Gartling, 1994.
Através da aplicação do método de elementos finitos o sistema de equações diferenciais
formado pela equação da conservação de massa e pelas equações de conservação da quantidade
de movimento (duas no caso bidimensional e três no caso tridimensional) é transformado em um
sistema algébrico passível de solução. No presente trabalho as equações de conservação
utilizadas são aquelas na forma exibida pelas equações (2.18) e (2.19). Essas equações governam
o escoamento, tridimensional (ou bidimensional), isotérmico, transiente e turbulento, das grandes
escalas de um fluido newtoniano com hipótese de quase-incompressibilidade escritas na forma
indicial e adimensionalizadas. Para facilitar a leitura do trabalho essas equações são repetidas
abaixo.
* *
* *
1
0
j
j
p u
t Mc x
∂ ∂
+ =
∂ ∂
(4.1)
* * * * * * *
*
* * * * * * * *
1 1 1 1
i j i j i k
i
j i j L T i j i k
u u u p u u u
f
t x Mc x x Re Re x x Re x x
λ
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
+ = − + + + + +
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
(4.2)
Para que o sistema formado por essas equações seja completo e passível de solução são
necessárias também condições iniciais e condições de contorno, apresentadas nas equações (4.3)
e (4.4) como em Reddy e Gartling, 1994.
* *
0
* *
0
em 0
em 0
i i
u u t
p p t
= =
= =
(4.3)
39
⌢
*
*
*
*
*
sobre
sobre
i i
D
ij
i
i j N
u u
T n T
σ
ρ
= Γ
= = Γ
⌢
(4.4)
Onde
*
0
i
u
e
*
0
p
são valores prescritos de velocidade e pressão para o tempo inicial.
*
ij
σ
é a
tensão no fluido,
*
i
T
é a tração sobre o fluido na direção i e
*
i
u
⌢
e
⌢
*
i
T
são os valores da velocidade
e da tração no fluido nos contornos, todos na forma adimensional. Γ representa o contorno do
domínio do problema chamado de Ω; o contorno é dividido em duas partes, uma chamada Γ
D
onde é aplicada uma condição de valor prescrito das variáveis primárias, chamada de condição
de primeira espécie (condição de Dirichlet) e outra chamada Γ
N
onde é aplicada uma condição de
segunda espécie (condição de Neumann) envolvendo derivadas das variáveis primárias. Deve-se
ressaltar que não existe sobreposição entre as fronteiras Γ
D
e Γ
N
, ou seja, Γ
D
∩ Γ
N
= 0. A Figura
4.1 mostra um esquema do domínio e suas fronteiras.
Figura 4.1. Esquema do domínio Ω e das fronteiras Γ
D
e Γ
N
para o método de elementos
finitos.
A tração de grandes escalas no fluido é obtida a partir da equação (4.5).
( )
i j k
ij
i
j ij T ij j
j i k
p u u u
T n n
x x x
σ
λ
δ ν ν δ
ρ ρ ρ
∂ ∂ ∂
= = − + + + +
∂ ∂ ∂
(4.5)
Na forma adimensional, e equação (4.5) é dada por:
* * *
*
*
*
* * *
1 1 1
i j k
i
ij j ij ij j
L T j i k
p u u u
T n n
Mc Re Re x x Re x
λ
σ δ δ
∂ ∂ ∂
= = − + + + +
∂ ∂ ∂
(4.6)
Ω
Γ
N
Γ
D
40
4.1. Forma de resíduos ponderados das equações de conservação de massa e quantidade de
movimento
Para aplicar a aproximação de elementos finitos sobre alguma equação diferencial é
necessário aplicar o princípio de resíduos ponderados sobre essa equação. Utilizando o processo
descrito em Reddy, 1985, a forma de resíduos ponderados de uma equação é obtida seguindo-se
os passos mostrados a seguir.
Primeiro todos os termos não-nulos das equações são passados para um dos lados da
igualdade, multiplica-se então cada equação por uma função peso e integram-se as equações
resultantes sobre o domínio. Essas etapas aplicadas sobre as equações (4.1) e (4.2) levam a:
* *
* *
1
0
j
j
p u
q d q d
t Mc x
Ω Ω
∂ ∂
Ω+ Ω =
∂ ∂
∫ ∫
(4.7)
* * * *
* * *
* * *
*
* * * * *
1
1 1 1
0
i j i
i i i
j i
j i k
i i i
i
j L T i j i k
u u u p
w d w d w d
t x Mc x
u u u
w d w d w f d
x Re Re x x Re x x
λ
Ω Ω Ω
Ω Ω Ω
∂ ∂ ∂
Ω+ Ω+ Ω−
∂ ∂ ∂
∂ ∂ ∂ ∂ ∂
+ + + Ω− Ω− Ω =
∂ ∂ ∂ ∂ ∂
∫ ∫ ∫
∫ ∫ ∫
(4.8)
Onde q é a função peso associada à pressão e w
i
a função peso associada a componente i da
velocidade. Ambas as funções peso são arbitrárias, devendo apenas satisfazer às condições de
primeira espécie na forma homogênea.
Após isto se distribui a diferenciação entre as variáveis do problema e as funções peso
(usando as regras de derivação do produto de funções ou de integração por partes, regra da
cadeia e teorema de Gauss-Green). Chega-se então nas equações (4.9) e (4.10).
* *
* *
1
0
j
j
p u
q d q d
t Mc x
Ω Ω
∂ ∂
Ω+ Ω =
∂ ∂
∫ ∫
(4.9)
* * *
* *
* *
* * * *
* * * *
1 1 1
0
i j i
i
i i i i
i
j
j i k
i
ij ij
j L T i j k
u u u
w d w d w T d w f d
t x
w
p u u u
d
x Mc Re Re x x Re x
λ
δ δ
Ω Ω Γ Ω
Ω
∂ ∂
Ω+ Ω− Γ− Ω+
∂ ∂
∂
∂ ∂ ∂
+ − + + + + Ω =
∂ ∂ ∂ ∂
∫ ∫ ∫ ∫
∫
(4.10)
41
A equação da continuidade mantém a forma da equação (4.7). Na equação de Navier-
Stokes a distribuição da diferenciação, o uso do teorema de Gauss-Green e da definição da
condição de Neumann da equação (4.6), alteraram a equação passando algumas das derivadas
para a função peso e fazendo aparecer um termo de integração no contorno que não estava
presente na equação original. O fato da condição de contorno de segunda espécie ser utilizada na
derivação da forma variacional das equações de conservação faz com que esse tipo de condição
de contorno seja naturalmente satisfeita na resolução do problema. Assim a condição de
contorno de segunda espécie é também chamada em elementos finitos de condição natural. A
condição de primeira espécie é chamada de condição essencial, já que é necessário algum valor
pré-determinado das variáveis primárias para que o problema seja bem posto.
Observando as equações (4.10), (4.4) e (4.6) chega-se à conclusão de que a imposição de
condições de contorno na variável de pressão não é feita de forma trivial. Como a pressão não
aparece nas condições de contorno essenciais e sim nas naturais, prescrever uma condição de
contorno de pressão requer trabalhar-se com integrais de contorno na formulação do problema.
No trabalho de Gresho e Sani, 1987, a questão da imposição das condições de contorno de
pressão é discutida mais detalhadamente. Como a idéia do presente trabalho é o desenvolvimento
de um código a ser melhorado com futuros trabalhos, esta questão não irá ser abordada
diretamente e na primeira versão do código não serão utilizadas condições de contorno de
pressão prescrita.
O uso da formulação de quase-incompressibilidade iguala os requerimentos de derivação
das funções que descrevem os campos de velocidade e pressão através do aparecimento de um
termo de derivada do campo de pressões na equação da continuidade. Assim pode-se utilizar
elementos de mesma ordem para ambos os campos. Usualmente a condição de Ladyzhenskaya-
Babuska-Brezzi impede o uso deste tipo de elemento, a aplicação da formulação de quase-
incompressibilidade “burla” essa imposição. O termo de derivada no tempo da pressão também
faz com que não apareçam zeros na diagonal principal do sistema algébrico a ser formado, o que
facilita a solução numérica do mesmo.
Utilizando essa forma fraca das equações de conservação será possível transformar o
sistema de equações íntegro-diferenciais, dado pelas equações (4.9) e (4.10), em um sistema
algébrico de equações. Para tanto é feita uma aproximação de elementos finitos sobre as
equações citadas.
42
4.2. Aproximação de elementos finitos e formação do problema matricial
A aproximação de elementos finitos do problema consiste em discretizar o domínio do
problema em um número de elementos e definir funções para as funções peso e para as funções
aproximação. As funções aproximação (ou funções teste) são as funções utilizadas para
representar o comportamento das variáveis primárias dentro do domínio de cada elemento.
As funções peso, por serem arbitrárias, são responsáveis pela formação do sistema de
equações. Já que para cada função peso diferente escolhida se obtém um conjunto de equações
diferente. Podendo assim igualar o número de incógnitas do problema.
A discretização do domínio transforma o domínio Ω do problema em uma coleção de sub-
domínios Ω
e
de cada elemento, os quais não apresentam sobreposição, exceto na fronteira. Cada
elemento apresenta um número de nós e em cada nó são definidos valores das variáveis
primárias. O valor dessas variáveis no restante do domínio do elemento são dados pelas funções
de forma e pelos valores nodais das variáveis. As funções de forma têm esse nome justamente
por definirem a “forma” dos campos no interior do elemento, que pode ser linear, quadrática,
etc., dependendo da ordem das funções de forma utilizadas. O número de nós no interior do
elemento também está associado à ordem das funções de forma.
As funções de forma utilizadas no método de elementos finitos possuem “suporte local”,
ou seja, possuem valores não nulos apenas no interior do elemento a qual estão associadas. Isso
simplifica a montagem do sistema global de equações já que cada nó será influenciado apenas
pelos nós vizinhos a ele.
Uma forma de simplificar ainda mais a aplicação do método é utilizar as mesmas funções
de forma para as funções peso e para as funções aproximação. Esse esquema é conhecido como
método de Bubnov-Galerkin [Hughes, 1987]. As funções utilizadas no presente trabalho são
mostradas abaixo.
Para as funções peso, as aproximações são dadas por:
m
m
i
q
w
ψ
ψ
=
=
(4.11)
Onde
ψ
m
é a função de forma associada ao nó m do elemento, ou seja, existe uma função
de forma associada a cada um dos nós que compõem o elemento.
Para as funções teste as funções de forma são dadas pela combinação linear das funções de
forma definidas no interior do elemento, como mostrado na equação (4.12).
43
*
1
*
1
*
1
*
1
N
m m
m
N
m m
i
i
m
N
m m
i
i
m
N
m m
i
i
m
p p
u u
T T
f f
ψ
ψ
ψ
ψ
=
=
=
=
=
=
=
=
∑
∑
∑
∑
(4.12)
Aqui p
m
,
m
i
u
,
m
i
T
e
m
i
f
são os valores da pressão e da componente i da velocidade, da
tração no fluido e das forças de corpo adimensionais no nó m. N é o número de nós que compõe
o elemento. Apesar de poder gerar alguma confusão com os símbolos utilizados para pressão e
velocidade nas equações básicas, por simplicidade de escrita e de leitura serão utilizados esses
símbolos. No contexto atual a diferença das definições é claramente percebida.
As funções de forma são definidas como variáveis apenas no espaço (constantes no tempo)
e os valores nodais variáveis apenas no tempo (constantes no espaço), como mostrado na
equação (4.13).
(
)
( )
( )
*
*
*
m m
m m
m m
i i
p p t
u u t
ψ ψ
=
=
=
x
(4.13)
Inserindo as aproximações das equações (4.11) e (4.12) nas equações de conservação na
forma fraca, dadas pelas equações (4.9) e (4.10), obtém-se:
1 1
* *
1
0
e e
N N
n n n n
j
m m
n n
j
p u
d d
t Mc x
ψ ψ
ψ ψ
= =
Ω Ω
∂ ∂
Ω+ Ω =
∂ ∂
∑ ∑
∫ ∫
(4.14)
44
1 1 1
* *
1 1
1 1
*
*
1 1
*
1
1 1
e e
e e
e
N N N
n n n n l l
i i j
m m
n n l
j
N N
m n n m n n
i i
n n
N N
n n n n
k
n n
ij ij
k
m
N N
j
n n n n
j i
n n
L T i j
u u u
d d
t x
T d f d
p u
Mc Re x
x
u u
Re Re x x
λ
ψ ψ ψ
ψ ψ
ψ ψ ψ ψ
ψ ψ
δ δ
ψ
ψ ψ
= = =
Ω Ω
= =
Γ Ω
= =
Ω
= =
∂ ∂
Ω+ Ω+
∂ ∂
− Γ− Ω+
∂
− + +
∂
∂
+
∂
∂ ∂
+ + +
∂ ∂
∑ ∑ ∑
∫ ∫
∑ ∑
∫ ∫
∑ ∑
∫
∑ ∑
*
0
d
Ω =
(4.15)
onde Ω
e
e Γ
e
são o domínio e a fronteira de cada elemento. Como as funções de forma de um
elemento se anulam fora deste a integral é definida somente no interior do elemento ao qual é
associada. Utilizando as propriedades das funções de forma e dos valores nodais das variáveis
mostradas na equação (4.13) e evidenciando derivações e somatórios chega-se nas equações
(4.16) e (4.17).
* *
1
1
0
e e
n n
N
m n m n
j
n
j
p
d d u
t Mc x
ψ
ψ ψ ψ
=
Ω Ω
∂ ∂
Ω⋅ + Ω⋅ =
∂ ∂
∑
∫ ∫
(4.16)
* *
1
* * *
* * * *
1 1
e e e
e e e
e e
n
n l
N
m n m l n m n n
i
j i i
l
j
m m n
ij ij
m n n n n n
i k
j j k
m n m n
n n
j i
L T j i j j
u
d u d u d T
t x
d f d p d u
Mc x Re x x
d u d u
Re Re x x x x
λ
ψ ψ
ψ ψ ψ ψ ψ
δ δ
ψ ψ ψ
ψ ψ ψ
ψ ψ ψ ψ
=
Ω Ω Γ
Ω Ω Ω
Ω Ω
∂
∂
Ω⋅ + Ω ⋅ − Γ⋅ +
∂ ∂
∂ ∂ ∂
− Ω⋅ − Ω⋅ + Ω⋅ +
∂ ∂ ∂
∂ ∂ ∂ ∂
+ + Ω⋅ + Ω⋅
∂ ∂ ∂ ∂
∑
∫ ∫ ∫
∫ ∫ ∫
∫ ∫
1
0
N
n=
=
∑
(4.17)
Assim para cada função peso, ou melhor, para cada valor de m, se obtém um conjunto de
equações diferente para as n variáveis. Como tomamos m=n, forma-se um conjunto de n
equações e n incógnitas. Escrevendo na forma matricial tem-se:
(
)
⋅ + ⋅ = +
M U K U U F S
ɺ
(4.18)
45
onde M é a matriz de massa do problema, K(U) é a matriz de rigidez (matriz de coeficientes)
não-linear,
U
ɺ
é o vetor de derivadas no tempo das incógnitas,
U
o vetor de incógnitas,
F
o vetor
de carregamento (forças externas) e
S
o vetor de forças de tração sobre a fronteira do elemento.
No sistema acima ainda permanecem duas dificuldades: primeiro, ainda existe um termo de
derivadas dos valores das incógnitas em relação ao tempo e segundo, a matriz de rigidez do
problema é não-linear, o que impede a solução direta do problema. Um tratamento especial a
cada uma desses problemas será dado nas seções subseqüentes. As matrizes e vetores da equação
(4.18), para o caso tridimensional, são definidos abaixo.
0 0 0
0 0 0
0 0 0
0 0 0
x
y
z
P
=
M
M
M
M
M
(4.19)
e
mn mn m n
i P
M M d
ψ ψ
Ω
= = Ω
∫
(4.20)
,
0 0
x x
y y
z z
= =
F S
F S
F S
F S
(4.21)
1
1
e
e
N
m m n n
i i
n
N
m m n n
i i
n
F d f
S d T
ψ ψ
ψ ψ
=
Ω
=
Γ
= Ω⋅
= Γ⋅
∑
∫
∑
∫
(4.22)
{
}
{ }
{ }
{ }
m
m
m
m
u
v
w
p
=
U
ɺ
ɺ
ɺ
ɺ
ɺ
(4.23)
46
{
}
{ }
{ }
{ }
m
m
m
m
u
v
w
p
=
U
(4.24)
Onde
m
u
ɺ
,
m
v
ɺ
,
m
w
ɺ
e
m
p
ɺ
são os valores das derivadas no tempo das variáveis nos nós e
u
m
,
v
m
,
w
m
e
p
m
são os valores das variáveis nos nós do elemento. A matriz de coeficientes é na
verdade o somatório de outras matrizes, como mostrado na equação (4.25).
(
)
(
)
= + + + +
K U A U C V G D
(4.25)
onde
A
(
U
) é a matriz de termos advectivos, sendo a parte não-linear do problema,
C
é a matriz
de termos viscosos,
V
a matriz de termos viscosos volumétricos,
G
a matriz de termos do
gradiente de pressão e
D
a matriz de termos do divergente de velocidade da equação da
continuidade. Cada matriz é definida como mostrado abaixo, novamente para o caso
tridimensional.
( )
0 0 0
0 0 0
0 0 0
0 0 0 0
x
y
z
=
A
A
A U
A
(4.26)
*
1
e
n l
N
mn m l
i j
l
j
A u d
x
ψ ψ
ψ
=
Ω
∂
= Ω
∂
∑
∫
(4.27)
0 0
0 0
,
0 0
0 0 0 0 0 0 0 0
xx xy xz xx xy xz
yx yy yz yx yy yz
zx zy zz zx zy zz
= =
C C C V V V
C C C V V V
C V
C C C V V V
(4.28)
47
* * * *
* *
1 1
1
e e
e
m n m n
mn
ij ij
L T j i j j
m n
mn
ij
i j
C d d
Re Re x x x x
V d
Re x x
λ
ψ ψ ψ ψ
δ
ψ ψ
Ω Ω
Ω
∂ ∂ ∂ ∂
= + Ω+ Ω⋅
∂ ∂ ∂ ∂
∂ ∂
= Ω
∂ ∂
∫ ∫
∫
(4.29)
0 0 0 0
0 0 0
0 0 0 0
0 0 0
,
0 0 0 0
0 0 0
0
0 0 0 0
x
y
z
x y z
= =
G
G
G D
G
D D D
(4.30)
*
*
1
1
e
e
m
mn n
i
i
n
mn m
i
i
G d
Mc x
D d
Mc x
ψ
ψ
ψ
ψ
Ω
Ω
∂
= − Ω
∂
∂
= Ω
∂
∫
∫
(4.31)
O sistema mostrado na equação (4.18) junto com as condições de contorno e iniciais
dadas nas equações (4.3) e (4.4) aproximam as equações de balanço sobre um domínio
discretizado
Ω
.
4.3.
Discretização do avanço no tempo
Existem diversas formas de se tratar o problema no avanço no tempo disponíveis para
aplicação no método de elementos finitos. O avanço no tempo, no caso atual, é na verdade é a
forma utilizada para aproximar o vetor de derivadas temporais das variáveis primárias que
multiplica a matriz de massa na equação (4.18).
Usualmente para a simulação de grandes escalas de escoamentos turbulentos se opta por
esquemas totalmente explícitos [Petry, 2002; Popiolek et al., 2006], o que simplifica esse avanço
no tempo (e também o problema da não-linearidade) e diminui o requerimento de armazenagem
de valores. O problema associado a um esquema totalmente explícito é que o passo de tempo
passível de ser utilizado é muito pequeno. Algumas vezes pesando mais no custo computacional
do que a economia de memória do esquema. Para a formulação explícita utilizada no presente
trabalho, em problemas de advecção dominante, esse fator é bastante crítico, já que o passo de
tempo está limitado a valores abaixo do obtido através da equação (4.32), sob a pena de
divergência da solução caso se use um valor superior.
48
U
i
máximo
x
t
c
∆
∆ ≤
+
(4.32)
onde
∆t
máximo
é o máximo valor do passo de tempo para que a solução tenha convergência
garantida e
∆x
i
é a menor dimensão dos elementos da malha.
O esquema de avanço no tempo escolhido para o presente trabalho é a aproximação de
família
θ
[Hughes, 1987; Reddy, 1985]. Esse esquema pode, a partir da escolha do parâmetro
θ
,
reproduzir esquemas explícito, semi-implícito e totalmente implícito. O esquema é baseado na
aproximação de uma média ponderada da derivada no tempo das variáveis dependentes em dois
passos de tempo consecutivos por uma interpolação linear dos valores dessas variáveis nos dois
passos de tempo. Para aplicar essa aproximação, primeiro se reescreve a equação (4.18) para dois
passos de tempo subseqüentes, já que esta equação é válida para qualquer instante de tempo.
Chega-se então nas equações (4.33) e (4.34).
(
)
n n n n
n
⋅ + ⋅ = +
M U K U U F S
ɺ
(4.33)
(
)
1 1 1 1
1
n n n n
n
+ + + +
+
⋅ + ⋅ = +
M U K U U F S
ɺ
(4.34)
Onde o subscrito
n
e
n
+1 são indicadores do instante de tempo de cada equação. Entre os
instantes de tempo
n
e
n
+1 existe uma diferença
∆t
(passo de tempo). A aproximação da família
θ
aproxima as derivadas no tempo nos instantes
n
e
n
+1 por:
( )
1
1
1
n n
n n
t
θ θ
+
+
−
⋅ + − ⋅ =
∆
U U
U U
ɺ ɺ
(4.35)
Onde
θ
é o parâmetro de avanço no tempo.
Isolando em cada uma das equações (4.33) e (4.34) os vetores de derivada no tempo das
variáveis dependentes, tem-se:
(
)
(
)
1
n n n n n
−
= ⋅ + − ⋅
U M F S K U U
ɺ
(4.36)
(
)
(
)
1
1 1 1 1 1
n n n n n
−
+ + + + +
= ⋅ + − ⋅
U M F S K U U
ɺ
(4.37)
49
Inserindo as equações (4.36) e (4.37) na equação (4.35) chega-se a forma do problema
apresentada na equação .
(
)
(
)
( ) ( )
( )
1 1 1 1
1
1
n n n n
n n n n n n
t
t
θ
θ
+ + + +
+
∆ ⋅ + − ⋅ +
+∆ − + − ⋅ = ⋅ − ⋅
F S K U U
F S K U U M U M U
(4.38)
Agrupando-se os termos para cada instante de tempo.
(
)
(
)
(
)
(
)
(
)
( )
( )
( )
( )
1 1
1 1
1
1 1
n n n n
n n n n
t t
t t t t
θ θ
θ θ θ θ
+ +
+ +
+∆ ⋅ ⋅ ⋅ = −∆ ⋅ − ⋅ ⋅ +
+ ∆ ⋅ ⋅ +∆ ⋅ − ⋅ + ∆ ⋅ ⋅ +∆ ⋅ − ⋅
M K U U M K U U
F F S S
(4.39)
Considerando que o problema transiente irá sempre partir de um campo inicial dado e que
as condições de contorno são bem definidas para todos os passos de tempo, o lado direito da
equação (4.39) é sempre conhecido. Assim pode-se resolver o sistema algébrico para obter-se os
campos de pressão e velocidade para o passo de tempo
n
+1.
A escolha do valor de
θ
define o tipo de aproximação no tempo. Em princípio,
θ
pode
assumirqualquer valor entre 0 e 1. Porém classicamente são utilizados os valores mostrados na
Tabela 4.1).
Valor de
θ
Nome do esquema
θ
= 0
Esquema totalmente explícito (de Euler) ou esquema de diferença adiantada
θ
= ½
Esquema de Crank-Nicholson semi-implícito
θ
=
⅔
Esquema de Galerkin semi-implícito
θ
= 1
Esquema totalmente implícito ou esquema de diferença atrasada
Tabela 4.1. Valores do parâmetro de avanço no tempo.
Os dois esquemas semi-implícitos são teoricamente incondicionalmente estáveis para
problemas lineares, enquanto os outros dois são condicionalmente estáveis [Reddy, 1985]. A
condição para esta estabilidade implica numa limitação no valor máximo de passo de tempo que
pode ser utilizado. Como mostrado em Lange, 1992, tomando-se
θ
= ½ (esquema semi-
implcícito de Crank-Nicholson) obtém-se um esquema de segunda ordem para a aproximação no
tempo.
Para o esquema explícito o limite no valor do passo de tempo para convergência
condicional é dado pela equação (4.40).
50
U
i
máximo
x
t
∆
∆ ≤
(4.40)
Como pode ser percebido pela análise desta equação, em comparação com a equação
(4.32), dependendo da relação entre a velocidade característica do problema e da velocidade do
som no fluido o valor do passo de tempo fornecido pela equação (4.40) é consideravelmente
maior do que aquele obtido a partir da equação (4.32). Isto significa uma economia
computacional em termos de número de passos de tempo para se atingir os valores estatísticos
desejáveis em uma simulação de grandes escalas.
Com essa aproximação o sistema de equações diferenciais dado pela equação (4.18) é
aproximado por um sistema de equações algébricas dado pela equação (4.39).
O sistema de equações algébricas obtido ainda não pode ser resolvido diretamente, já que
ainda apresenta termos não-lineares.
4.4.
Método de Newton para sistemas de equações algébricas não-lineares
Para tratar o problema da não-linearidade do sistema de equações algébricas formado na
aproximação de elementos finitos das equações de conservação, é utilizado o método de Newton
[Reddy e Gartling, 1994]. Este método resolve um sistema algébrico não-linear através de um
processo iterativo, resolvendo a cada passo um sistema linear. O método de Newton é um
método de segunda ordem e fornece uma solução geralmente com um número pequeno de
iterações (convergência rápida) porém requer um campo inicial adequado (pequeno raio de
convergência).
Sendo uma extensão do método de Newton-Raphson para encontrar raízes de equações
não-lineares, o método é baseado em uma expansão em séries de Taylor truncada do resíduo.
Esse resíduo para o problema em questão é dado por:
(
)
(
)
(
)
(
)
(
)
( )
( )
( )
( )
1 1
1 1
( ) 1
1 1
n n n n
n n n n
t t
t t t t
θ θ
θ θ θ θ
+ +
+ +
= +∆ ⋅ ⋅ ⋅ − −∆ ⋅ − ⋅ ⋅ +
− ∆ ⋅ ⋅ +∆ ⋅ − ⋅ − ∆ ⋅ ⋅ + ∆ ⋅ − ⋅
R U M K U U M K U U
F F S S
(4.41)
Onde
R
(
U
) é o vetor de valores residuais da equação matricial aproximada do problema. O
valor do vetor de incógnitas
U
n+1
que anula esse resíduo é a solução do sistema não-linear.
Expandindo o resíduo em séries de Taylor ao redor da solução U
p
tem-se:
51
2
( ) ( ) 0
p
p
O
∂
+ ⋅∆ + ∆ =
∂
U
R
R U U U
U
(4.42)
onde U
p
é a solução da iteração
p
do processo iterativo do método de Newton,
∆U
é a variação
da solução entre as iterações
p
e
p+
1, ou seja,
∆U
=(
U
p+1
-
U
p
), e
O
(
∆U
2
) é o erro de truncamento
da ordem de
∆U
2
. Omitindo-se os termos de segunda ordem e superiores chega-se à equação
(4.43).
( ) ( ) ( )
1 1
( )
p
p p p p p p
+ +
∂
⋅ − ≡ ⋅ − = −
∂
U
R
U U J U U U R U
U
(4.43)
onde a matriz
J
é chamada de matriz jacobiana, ou matriz tangente. Essa matriz contém os
termos da derivada do resíduo em relação as variáveis do problema. Para o problema em questão
o sistema não-linear é aproximado pelo processo iterativo de se obter uma solução
U
p+1
para o
sistema linear dado pela equação, para o qual o valor de
∆U
seja menor do que um critério de
convergência.
(
)
(
)
1
1 1 1 1
( , )
p p p p
n n n n n
+
+ + + +
⋅ − = −
J U U U R U U
(4.44)
onde os subscritos
p
e
p
+1 se referem ao número da iteração no método de Newton e os
subscritos
n
e
n
+1 se referem ao passo de tempo no processo transiente. A matriz jacobiana é
obtida a partir da expressão abaixo.
(
)
(
)
(
)
(
)
1 1 1
p p p
n n n
t
θ
+ + +
′
= + ∆ ⋅ +J U M K U A U
(4.45)
Onde
′
A
é a matriz da derivada da matriz de termos advectivos em relação as variáveis do
problema, dada pelas equações (4.46) e (4.47).
( )
1
0
0
0
0 0 0 0
xx xy xz
yx yy yz
p
n
zx zy zz
+
′ ′ ′
′ ′ ′
′
=
′ ′ ′
A A A
A A A
A U
A A A
(4.46)
52
*
1
e
n l
N
mn m l
ij i
l
j
A d u
x
ψ ψ
ψ
=
Ω
∂
′
= Ω
∂
∑
∫
(4.47)
A solução transiente do problema é então obtida a partir dos seguintes passos:
1.
A partir das condições iniciais e de contorno monta-se o primeiro passo do sistema linear
dado pela equação (4.44), com
1
p
n n
+
=
U U
;
2.
Resolve-se o sistema e calcula-se
1
1 1 1
p p p
n n n
+
+ + +
= ∆ +
U U U
;
3.
O valor de
1
1
p
n
+
+
U
obtido passa a ser
1
p
n
+
U
e monta-se novamente o sistema linear;
4.
Resolve-se o sistema novamente e recalcula-se
1
p
n
+
∆U
, se este valor estiver abaixo do critério
de convergência a solução é considerada convergida para esse passo de tempo, se o valor
obtido estiver acima do critério de convergência volta-se ao passo 3;
5.
O valor obtido de
1
1
p
n
+
+
U
passa a ser
n
U
e o processo é reiniciado a partir do passo 1 para o
novo passo de tempo.
4.5.
Oscilações no campo de pressão e diagonalização seletiva da matriz de massa
Um dos problemas enfrentados quando se utiliza o método de elementos finitos para
aproximar as equações de conservação para simulação de escoamentos é a compatibilização dos
campos de pressão e velocidade e a instabilidade do campo de pressão. No presente trabalho a
primeira dificuldade já foi superada pela aplicação da formulação de quase-incompressibilidade,
a segunda, porém, ainda requer algum tratamento especial.
Uma das formas de instabilidade do campo de pressão, que ocorre mesmo quando se
utilizam elementos de baixa ordem para o campo de pressão (respeitando a condição de
Ladyzhenskaya-Babuska-Brezzi) é o efeito “tabuleiro de xadrez” (conhecido pela expressão em
inglês
checkerboard
, que significa “tabuleiro de damas”) que recebe este nome pela aparência do
campo de pressão obtido com a instabilidade. O problema desses modos de pressão instáveis é
discutido em Gresho e Sani, 1999, onde algumas soluções também são apontadas. Selvam, 1997,
utiliza em seu trabalho um esquema de solução alternada dos campos de velocidade e pressão em
quatro passos, conseguindo assim evitar a instabilização da pressão.
Em seu trabalho, Lee et al., 1979 sugerem algumas técnicas de suavização dos campos de
pressão (e vorticidade) entre cada passo de tempo ou entre cada iteração do método de solução
53
do sistema não-linear. Diferentes formas de suavização são testadas. Todas com resultados
próximos porém com variados custos computacionais.
Uma outra alternativa para a estabilização do sistema é o uso de diagonalização seletiva da
matriz de massa. Esse método é utilizado nos trabalhos anteriores do grupo [Petry, 2002] e é
discutido e utilizado no trabalho de Kawahara e Hirano, 1983. A diagonalização consiste em
concentrar os valores obtidos na matriz de massa do problema na diagonal principal da mesma,
através de algum tipo de média, conservando a massa do elemento. Essa técnica também é
chamada de amortecimento da matriz de massa (já que esta é a matriz mais característica do
problema transiente). O uso de alguma forma de diagonalização da matriz de massa pode trazer
consigo uma certa perda na precisão numérica da solução, porém é conveniente pelo ganho em
eficiência e estabilidade do esquema de solução, o que compensa a teórica perda numérica.
Algumas vezes a diagonalização da matriz de massa pode até aumentar a precisão numérica do
sistema. Isso ocorre em casos onde a não-utilização desse tipo de subterfúgio faz com que
apareçam oscilações numéricas (instabilidades como o efeito de tabuleiro de xadrez). Tais
oscilações são amortecidas pelo uso de esquemas de diagonalização da matriz de massa,
melhorando assim a precisão do método.
Diferentes alternativas de diagonalização são utilizadas; algumas são mostradas no
trabalho de Zienkiewicz e Taylor, 2000. A alternativa aplicada em trabalhos anteriores do grupo
vem do trabalho de Kawahara e Hirano, 1983, em que a parte da matriz de massa referente às
variáveis de velocidade é totalmente diagonalizada e aquela referente às variáveis de pressão são
diagonalizadas de forma seletiva, o que leva ao aparecimento de um novo parâmetro arbitrário
para o problema.
Durante a elaboração do código utilizado no presente trabalho, foram testadas duas formas
de diagonalização da matriz de massa, aquela igual à citada de Kawahara e Hirano, 1983, e uma
nova proposta de diagonalização completa apenas da parte da matriz de massa referente aos
termos de pressão.
A diagonalização seletiva proposta por Kawahara e Hirano, 1983, se dá pelo uso de duas
matrizes de massa adicionais, uma consistente, uma totalmente diagonalizada e uma
diagonalizada seletivamente, como mostradas nas equações (4.48) e (4.50).
0 0 0
0 0 0
0 0 0
0 0 0
x
y
z
P
=
M
M
M
M
M
(4.48)
54
1
e
N
mn mn m n
i P mn
n
M M d
ψ ψ δ
=
Ω
= = Ω⋅
∑
∫
(4.49)
0 0 0
0 0 0
0 0 0
0 0 0
x
y
z
P
=
M
M
M
M
M
ɶ
ɶ
(4.50)
( )
1
1
e e
N
mn m n m n
P mn
n
M e d e d
ψ ψ δ ψ ψ
=
Ω Ω
= ⋅ Ω⋅ + − ⋅ Ω
∑
∫ ∫
ɶ
(4.51)
onde as
M
é a matriz de massa totalmente diagonalizada e
M
ɶ
a matriz de massa diagonalizada
seletivamente. As sub-matrizes
i
M
(diagonais) e
P
M
ɶ
(diagonalizada seletivamente) são
calculadas como mostrado nas equações (4.49) e (4.51). Na equação (4.51)
e
é o parâmetro de
diagonalização seletiva que pode variar entre os valores
0 1
e
≤ ≤
(geralmente adota-se e=0,7 ou
e=0,8).
Utilizando esta primeira forma de diagonalização seletiva, o sistema de equações a ser
resolvido em cada iteração do método de Newton se mantém o mesmo.
(
)
(
)
1
1 1 1 1
( , )
p p p p
n n n n n
+
+ + + +
⋅ − = −
J U U U R U U
(4.52)
Porém com as definições de
R
e
J
alteradas para:
(
)
(
)
( ) ( )
(
)
( )
( )
( )
( )
1 1 1
1 1
( ) 1
1 1
p p
n n n n n
n n n n
t t
t t t t
θ θ
θ θ θ θ
+ + +
+ +
= +∆ ⋅ ⋅ ⋅ − −∆ ⋅ − ⋅ ⋅ +
− ∆ ⋅ ⋅ +∆ ⋅ − ⋅ − ∆ ⋅ ⋅ + ∆ ⋅ − ⋅
R U M K U U M K U U
F F S S
ɶ
(4.53)
(
)
(
)
(
)
(
)
1 1 1
p p p
n n n
t
θ
+ + +
′
= + ∆ ⋅ +
J U M K U A U (4.54)
A segunda forma de diagonalização da matriz de massa, proposta no trabalho, altera
apenas as linhas da matriz de massa referentes às variáveis de pressão. Assim o sistema se
mantém inalterado para as variáveis de velocidade, mantendo sua precisão, e as variáveis de
55
pressão são amortecidas como desejado. A matriz de massa diagonalizada na pressão
M
⌢
fica
então na forma:
0 0 0
0 0 0
0 0 0
0 0 0
x
y
z
P
=
M
M
M
M
M
⌢
⌢
(4.55)
Com os termos da matriz de massa diagonalizada na pressão
P
M
⌢
calculados como
mostrado pela equação (4.56).
1
e
N
mn m n
P mn
n
M d
ψ ψ δ
=
Ω
= Ω⋅
∑
∫
⌢
(4.56)
O sistema de equações lineares a ser resolvido a cada iteração do método de Newton, dado
pela equação (4.44), permanece inalterado. Porém, assim como na primeira forma de
diagonalização seletiva apresentada as definições de R e J são alteradas para:
(
)
(
)
( ) ( )
(
)
( )
( )
( )
( )
1 1 1
1 1
( ) 1
1 1
p p
n n n n n
n n n n
t t
t t t t
θ θ
θ θ θ θ
+ + +
+ +
= +∆ ⋅ ⋅ ⋅ − −∆ ⋅ − ⋅ ⋅ +
− ∆ ⋅ ⋅ +∆ ⋅ − ⋅ − ∆ ⋅ ⋅ + ∆ ⋅ − ⋅
R U M K U U M K U U
F F S S
⌢
(4.57)
(
)
(
)
(
)
(
)
1 1 1
p p p
n n n
t
θ
+ + +
′
= + ∆ ⋅ +
J U M K U A U (4.58)
As equações (4.57) e (4.58) mostram que a forma de diagonalização proposta além de
preservar as equações das variáveis de velocidade ainda mantém a forma original da matriz de
massa multiplicando os valores do passo de tempo anterior, conservando ainda mais a precisão e
coerência dos resultados.
4.6.
Método de Picard para sistemas de equações algébricas não-lineares
Uma alternativa para solução do problema da não-linearidade do sistema de equações
mostrado na equação (4.39) é o método de Picard [Reddy e Gartling, 1994], ou método da
substituição progressiva. Esse método, diferente do método de Newton, é um método de primeira
56
ordem (com convergência e precisão um pouco menores), mas que apresenta um maior raio de
convergência quando comparado ao método de Newton.
O raio de convergência é importante no sentido de que se o campo inicial fornecido para o
início do processo iterativo estiver muito longe da solução real, um método com um pequeno
raio de convergência pode não conseguir resolver o problema. Assim o método de Picard, que
apresenta um raio de convergência maior que o método de Newton, pode ser utilizado em casos
onde o segundo falhe nos passos iniciais. Assim um resultado intermediário obtido com o
método de Picard pode ser passado para o método de Newton como um campo inicial mais
próximo da solução real, facilitando a convergência do método de Newton.
Para o método de Picard o processo de busca da solução se dá da seguinte maneira: Monta-
se o sistema dado pela equação (4.39) utilizando o campo inicial fornecido como primeira
aproximação para o campo solução,
1
p
n n
+
=
U U
, na forma mostrada na equação (4.59).
(
)
(
)
( ) ( )
(
)
( )
( )
( )
( )
1
1 1
1 1
1
1 1
p p
n n n n
n n n n
t t
t t t t
θ θ
θ θ θ θ
+
+ +
+ +
+∆ ⋅ ⋅ ⋅ = −∆ ⋅ − ⋅ ⋅ +
+ ∆ ⋅ ⋅ +∆ ⋅ − ⋅ + ∆ ⋅ ⋅ + ∆ ⋅ − ⋅
M K U U M K U U
F F S S
(4.59)
Como pode ser observado, na equação (4.59) a única incógnita é o valor de
1
1
p
n
+
+
U
.
Resolvendo-se então o sistema linear se obtém um novo valor para
1
1
p
n
+
+
U
. Este valor é comparado
com o valor da iteração anterior p,
1
p
n
+
U
, e caso tenha sido atingido o critério de convergência se
passa para o próximo passo de tempo. Caso não tenha sido atingida a convergência se utiliza o
valor obtido de
1
1
p
n
+
+
U
como
1
p
n
+
U
e se reinicia o processo.
Para o método de Picard também foi implementada a diagonalização seletiva da matriz de
massa, o que simplesmente reescreve a equação (4.59) na equação (4.60) para a diagonalização
seletiva de Kawahara e Hirano, 1983, ou na equação (4.61) para a diagonalização proposta no
trabalho.
(
)
(
)
( ) ( )
(
)
( )
( )
( )
( )
1
1 1
1 1
1
1 1
p p
n n n n
n n n n
t t
t t t t
θ θ
θ θ θ θ
+
+ +
+ +
+∆ ⋅ ⋅ ⋅ = −∆ ⋅ − ⋅ ⋅ +
+ ∆ ⋅ ⋅ +∆ ⋅ − ⋅ + ∆ ⋅ ⋅ + ∆ ⋅ − ⋅
M K U U M K U U
F F S S
ɶ
(4.60)
(
)
(
)
( ) ( )
(
)
( )
( )
( )
( )
1
1 1
1 1
1
1 1
p p
n n n n
n n n n
t t
t t t t
θ θ
θ θ θ θ
+
+ +
+ +
+∆ ⋅ ⋅ ⋅ = −∆ ⋅ − ⋅ ⋅ +
+ ∆ ⋅ ⋅ +∆ ⋅ − ⋅ + ∆ ⋅ ⋅ + ∆ ⋅ − ⋅
M K U U M K U U
F F S S
⌢
(4.61)
57
Onde a definição das matrizes de massa totalmente diagonal
M
, diagonalizada
seletivamente
M
ɶ
e diagonalizada na pressão
M
⌢
continuam as mesmas dadas nas equações
(4.48), (4.50) e (4.55) respectivamente.
4.7.
Método explícito simplificado e explícito diagonalizado
Outra possibilidade de avanço no tempo sem o problema da não-linearidade do sistema é o
uso de
θ
=0. Que caracteriza o método totalmente explícito. Nesse método de avanço no tempo o
valor dos campos de velocidade e pressão no passo de tempo futuro são dependentes apenas dos
valores desses mesmos campos no passo de tempo anterior, eliminando a não-linearidade.
Reescrevendo o problema dado pela equação (4.39) para
θ
=0 tem-se:
(
)
(
)
(
)
1
n n n n n
t t t
+
⋅ = −∆ ⋅ ⋅ + ∆ ⋅ ∆ ⋅
M U M K U U F S
(4.62)
Como visto na equação (4.62) o método explícito de avanço no tempo simplifica a
montagem do sistema de equações. Utilizando a diagonalização seletiva pode-se também
eliminar a necessidade de um método iterativo, ou de um método direto mais caro, para a
solução do sistema linear, que passa a ser simplesmente uma operação algébrica.
Aplicando sobre a equação (4.62) as duas formas de diagonalização seletiva utilizadas no
trabalho são obtidas as equações (4.63) e (4.64).
(
)
(
)
(
)
1
n n n n n
t t t
+
⋅ = −∆ ⋅ ⋅ + ∆ ⋅ ∆ ⋅
M U M K U U F S
ɶ
(4.63)
(
)
(
)
(
)
1
n n n n n
t t t
+
⋅ = −∆ ⋅ ⋅ + ∆ ⋅ ∆ ⋅
M U M K U U F S
⌢
(4.64)
onde na equação (4.63) é utilizada a diagonalização seletiva de Kawahara e Hirano, 1983, e na
equação (4.64) é utilizada a matriz de massa diagonalizada na pressão, proposta no trabalho.
Na equação (4.63) percebe-se que cada incógnita é multiplicada apenas por um valor da
matriz de massa totalmente diagonalizada que aparece do lado esquerdo. Desta forma não é
necessária nenhuma técnica para solução do sistema linear, sendo somente preciso dividir o lado
direito da equação pelo valor da diagonal principal da matriz de massa
M
na mesma linha.
58
A principal limitação dos métodos totalmente explícitos é a valor máximo do passo de
tempo permitido para o problema para o qual a convergência é assegurada. Utilizando as
formulações discutidas nessa seção, a literatura e trabalhos anteriores do grupo mostram que o
limite para o passo de tempo é dado pela equação (4.32).
A idéia de ter esses métodos todos implementados no código é, de acordo com a
necessidade, utilizar um esquema com convergência mais segura para casos mais problemáticos,
ou esquemas de avanço mais rápido do tempo em casos que permitam esse tipo de tratamento.
4.8.
Solução do sistema linear
Em todos os métodos discutidos acima, exceto no método explícito com a diagonalização
seletiva, aparece em algum estagio da aproximação um sistema de equações algébricas lineares
que precisa ser resolvido. Inúmeras técnicas são utilizadas para a resolução desse tipo de sistema
e, dentre essas, duas foram escolhidas para serem implementadas no método. São eles o método
da sobre-relaxação sucessiva (conhecido pela sigla em inglês SOR, successive over-relaxation) e
o método dos gradientes conjugados. A maioria dos livros sobre métodos dinâmica dos fluidos
computacionais traz algum capítulo sobre métodos de solução de sistemas lineares; os métodos
utilizados no presente trabalho são bem discutidos nos trabalhos de Ferziger e Peri
ć
, 2002, e no
livro de ferramentas numéricas de Press et al., 1992.
Ambos os métodos selecionados são métodos iterativos de solução. Métodos iterativos
apresentam, geralmente, um custo computacional menor do que os métodos diretos, sendo
bastante atraentes para uso nas aplicações de dinâmica dos fluidos computacional.
59
5.
ASPECTOS COMPUTACIONAIS
Como já apresentado, o objetivo principal do trabalho é elaborar um código numérico que
será utilizado em trabalhos futuros do grupo. Cada parte em separado desse código base pode ser
motivo de um novo estudo para sua otimização e implementação de novas técnicas e modelos.
Com essa finalidade esse código base deve apresentar aspectos computacionais que
permitam a remodelagem e novas implementações de forma fácil.
Nesse sentido boa parte da atenção desse trabalho de criação do novo código foi voltada
para sua estruturação de forma adequada e ao uso de técnicas computacionais de alto
desempenho. Nas próximas seções serão mostrados alguns itens aplicados nesse trabalho de
estruturação e desenvolvimento.
5.1.
Estruturação do novo código
As principais características do novo código desenvolvido são: estruturas de programação
e dados; armazenagem compacta de matrizes; e utilização de processamento paralelo.
A estrutura de programação é tal que possibilita a implementação de novos modelos de
turbulência, técnicas de montagem do sistema de equações e métodos de solução de forma
simplificada. Pra cada etapa do processamento existe uma coleção de sub-rotinas separadas por
módulos (em um arquivo distinto). Assim uma alteração em uma parte específica do programa
pode ser feita editando-se apenas a parte correspondente no código, preservando-se entradas e
saídas de dados como especificado. Essa estrutura também simplifica a verificação de cada
alteração feita no programa. Todos os dados lidos e escritos (de arquivos externos) estão listados
no código e programa gera saídas em seu formato padrão e nos formatos lidos pelos programas
TecPlot (através de uma sub-rotina interna) e ParaView (através de um programa conversor
externo).
A estrutura completa do código pode ser visualizada através do diagrama na figura abaixo:
60
Figura 5.1. Esquema da estrutura do programa dividido por módulos.
Detalhes sobre a armazenagem compacta de matrizes esparsas e a implementação do
processamento paralelo serão discutidos nas próximas seções.
5.2.
Armazenagem compacta de matrizes
De acordo com o modelo numérico apresentado no capítulo 4, o programa desenvolvido
pode trabalhar com esquemas de avanço no tempo semi-implícitos e explícitos. Como mostrado
na formulação desses esquemas semi-implícitos, a armazenagem das matrizes globais completas
do sistema de elementos finitos é necessária para sua aplicação.
O problema da armazenagem completa das matrizes formadas por métodos numéricos
como diferenças, volumes e elementos finitos é o tamanho dessas matrizes. No caso de
elementos finitos as matrizes na forma completa irão conter n×n posições, onde n é o número
total de variáveis (número de nós multiplicado pelo número de variáveis por nó) do sistema. Isso
inviabiliza completamente o uso da armazenagem completa de matrizes. Porém uma das
características dessas matrizes é que elas são bastante esparsas (chegando a apresentar casos com
apenas 5% de posições não-nulas). Pode-se então fazer uso desta esparsidade e armazenar apenas
as posições não-nulas das matrizes, reduzindo drasticamente a requisição de memória.
61
Diversas técnicas para armazenagem compacta de matrizes foram desenvolvidas, dentre
elas as técnicas Skyline [Hughes, 1987] e armazenagem indexada por linhas [Press et al., 1992].
Essas técnicas consistem em formas de se armazenar as matrizes de forma reduzida (eliminando
o máximo possível de posições nulas), de uma maneira que seja possível a reconstrução da
matriz completa ou a identificação das posições de elementos da matriz completa na
armazenagem compacta.
Dentre as técnicas analisadas, a escolhida para ser implementada no código foi a
armazenagem compacta de matrizes esparsas indexadas por linhas. Esta escolha se baseou em
dois principais fatores: o método é extremamente eficiente, armazenando apenas a diagonal
principal da matriz e os valores não-nulos fora da diagonal; e o sistema de armazenagem permite
a identificação fácil e rápida das posições dos valores não-nulos na matriz completa, dados os
vetores de armazenagem. Esse segundo fator é de extrema importância para um bom
desempenho na manipulação das matrizes de forma compacta e na solução do sistema linear
formado.
Para exemplificar o uso deste método será reproduzido aqui o exemplo apresentado em
[Press et al., 1992]. Dada uma matriz A de tamanho N×N, por exemplo:
3 0 1 0 0
0 4 0 0 0
0 7 5 9 0
0 0 0 0 2
0 0 0 6 5
=
A (5.1)
Se utilizam as seguintes regras para a geração dos vetores para armazenagem compacta:
6.
As primeiras N posições do vetor de valores sa armazenam os valores da diagonal principal da
matriz A em ordem;
7.
Cada uma das N primeiras posições do vetor de índices ija armazena o índice no vetor sa que
contém o primeiro valor não-nulo de A na linha N. Caso a linha não contenha valores fora da
diagonal este valor é igual último índice armazenado de sa mais 1;
8.
A posição 1 do vetor ija é igual a N+2 (pode ser lido posteriormente para determinar N);
9.
A posição N+1 do vetor ija é igual à posição do último elemento fora da diagonal armazenado
no vetor sa mais 1 (pode ser lido para obter o número de elementos fora da diagonal da
matriz A ou o tamanho dos vetores sa e ija). A posição N+1 em sa não é utilizada e pode ser
definida arbitrariamente;
62
10.
As posições N+2 em sa contêm os valores não-nulos fora da diagonal da matriz A ordenados
por linhas;
11.
As posições N+2 em ija contêm os números das colunas na matriz A onde os valores fora da
diagonal armazenados em sa correspondentes são localizados.
Aplicando esses passos na matriz A utilizada como exemplo obtemos os seguintes vetores
de índices ija e de valores sa:
(5.2)
A regra número 1 impõe que os valores da diagonal principal da matriz devem ser
armazenados, mesmo se forem nulos. Isso caracterizaria uma certa ineficiência do método.
Porém, como no método numérico utilizado não aparecem valores nulos na diagonal principal
das matrizes, isso não gera nenhuma armazenagem desnecessária. Tais valores nulos na diagonal
principal são inclusive indesejados e, como já discutido, o uso da aproximação de quase-
incompressibilidade evita seu aparecimento.
O exemplo mostrado não apresenta uma grande redução na quantidade de informações
armazenadas, mas para grandes matrizes esparsas o método é extremamente eficiente.
Como se pode perceber, os valores da diagonal principal da matriz A são facilmente
acessados no vetor de valores. Outra facilidade é que os valores fora da diagonal para uma
mesma linha estão todos em um único bloco e ordenados, com o número da coluna que ocupam
na matriz completa na posição correspondente no vetor ija. Estes dois pontos simplificam a
implementação de algoritmos de solução do sistema linear. Como aqueles utilizados nesse
trabalho.
Na implementação dessa técnica no código desenvolvido foram criados algoritmos de
busca para definir as posições de cada valor calculado diretamente na matriz compacta. Assim,
após o cálculo das matrizes de elemento, o sistema de equações global já é montado diretamente
na estrutura compacta. Esses algoritmos de localização ainda necessitam de melhoras no seu
desempenho. Porém, acompanhando o tempo de processamento de cada etapa do programa
verificou-se que essa etapa não é a mais cara computacionalmente, viabilizando seu uso.
A técnica de armazenagem de matrizes compactas indexadas por linhas é bastante
apropriada para o uso com elementos finitos. Principalmente pelo fato de que todas as matrizes a
serem armazenadas possuem termos não-nulos nas mesmas posições. E essas posições podem ser
obtidas numa etapa de pré-processamento, já que são apenas dependentes da malha utilizada
índice k 1 2 3 4 5 6 7 8 9 10 11
ija(k) 7 8 8 10 11 12 3 2 4 5 4
sa(k) 3 4 5 0 5 x 1 7 9 2 6
63
(conectividades). O vetor ija pode então ser calculado antes do processamento do problema e
armazenado uma única vez.
Essa semelhança entre as matrizes simplifica e muito uma série de outras operações entre
as matrizes compactas armazenadas. Por exemplo, a soma de duas matrizes compactas com
vetores de índices idênticos é apenas a soma dos vetores. Operações como produto das matrizes
compactas por um escalar são bastante fáceis e rápidas.
5.3.
Implementação de processamento paralelo utilizando a técnica OpenMP
A parte mais pesada do código é o cálculo das matrizes de elemento. Boa parte dessas
matrizes deve ser recalculada a cada iteração (dentro de cada passo de tempo) para gerar um
novo sistema global a ser resolvido (o sistema global depende do esquema de solução escolhido,
como mostrado no capítulo 4). Para esta etapa foi então implementado o processamento em
paralelo utilizando a técnica OpenMP. Como descrito no livro de Chandra, 2001, essa técnica de
paralelização é bastante simples e eficiente, porém é limitada ao uso em máquinas de memória
compartilhada (onde todos os processadores têm acesso à mesma memória RAM).
A aplicação dessa técnica consiste no uso de diretivas escritas dentro do código que serão
interpretadas pelo compilador (que deve ter suporte à OpenMP) quais tarefas deverão ser
distribuídas e como será feita esta distribuição. O interessante dessa técnica é que quando se
executa o programa, ele automaticamente verifica a disponibilidade de processadores e distribui
as tarefas de forma dinâmica, ajustando o particionamento ao número de processadores e a carga
de cada processador.
No código desenvolvido, o laço que calcula as matrizes de elemento varre todos os
elementos da malha, calcula as matrizes e soma o resultado no sistema global. Como as matrizes
de elemento são independentes entre si e a montagem do sistema global só envolve soma de
valores em posições específicas, essas operações (cálculo das matrizes de elemento) são
independentes de ordem. A paralelização do laço de cálculo das matrizes de elemento divide o
número de elementos entre os processadores disponíveis. Cada processador fica responsável pelo
cálculo das matrizes de elemento de uma faixa dos elementos da malha, e o resultado de cada
matriz de elemento é passado para o sistema global. Nota-se então que a paralelização
implementada dessa forma não apresenta sobreposição de trabalho entre os processadores,
melhorando a eficiência do método.
A Figura 5.2 mostra um exemplo de trecho do código com as diretivas de compilação para
paralelização com OpenMP
64
Figura 5.2. Trecho do código mostrando as diretivas de compilação in-line para o uso de
paralelização.
Como mostrado na figura, com apenas duas linhas adicionais de diretivas de compilação o
processamento paralelo já estava implementado. Devido à sintaxe das linhas de diretiva de
compilação o código pode ser compilado utilizando compiladores que não têm suporte à técnica
OpenMP sem problema algum, já que as linhas referentes ao seu uso serão interpretadas como
comentários.
Fora o acréscimo das diretivas no interior do código, na compilação devem ser utilizados
argumentos para o compilador que ativam o suporte ao OpenMP e procuram no código as
diretivas.
O processamento dos casos apresentados nesse trabalho foi feito em duas máquinas, uma
máquina Prism da SGI e uma máquina Altix também da SGI. Ambas as estações possuíam 8
processadores com compartilhamento de memória. Não foi necessária qualquer alteração no
código para a compilação nas diferentes máquinas com diferentes sistemas operacionais (SuSe 9
no caso da Prism e RedHat no caso da Altix).
65
6.
RESULTADOS
Dois casos de referência, com resultados disponíveis na literatura foram escolhidos para
validação do código computacional elaborado. Os dois casos são: o caso do escoamento dentro
de uma cavidade confinada e o escoamento sobre um degrau.
O problema da cavidade é um caso de referência utilizado na validação de muitos trabalhos
publicados, como por exemplo, nos artigos recentes de Elias et al., 2004, Pereira e Campos Silva,
2005, e Costa et al., 2005. No presente trabalho o caso da cavidade foi simulado para as
condições bidimensional laminar e tridimensional turbulenta.
6.1.
Escoamento em uma cavidade bidimensional
O problema do escoamento no interior de uma cavidade consiste em um fluido confinado
dentro de uma cavidade (quadrada no caso) que se movimenta devido a uma condição de
velocidade constante forçada em uma das paredes.
Para o caso de uma cavidade quadrada de lado igual a 1m, resultados bastante confiáveis
foram publicados por Ghia et al., 1982, obtidos através de uma solução de diferenças finitas
utilizando métodos multi-grid. Esses resultados vêm sendo utilizados como referência e desde
então são reproduzidos com sucesso em outros trabalhos, como Schreiber e Keller, 1983, e Petry,
2002. A Figura 6.1 mostra o domínio, as condições de contorno e condições iniciais utilizadas
para o problema.
66
Figura 6.1. Esquema com descrição do problema da cavidade bidimensional, valores de
condições de contorno, iniciais e propriedades do fluido.
Na figura u e v são as velocidades impostas nas paredes nas direções x e y,
respectivamente, U é a velocidade de referência para o caso e os valores u, v e p definidos no
tempo t = 0s são as condições iniciais utilizadas para o problema. As propriedades do fluido
massa específica (
ρ
), viscosidade volumétrica (
λ
) e velocidade do som (c) são definidas, a
viscosidade absoluta (
µ
) é deixada como variável para que se possa variar o número de Reynolds
através da mudança de seu valor. O número de Reynolds para este caso é mostrado na equação
(6.1).
UL
cavidade
Re
ρ
µ
=
(6.1)
onde L é o comprimento característico da cavidade, no caso L = 1m.
O problema foi resolvido para valores de Re
cavidade
de 100, 400 e, 1000. Todos em regime
laminar e com resultados publicados por Ghia et al., 1982, e outros citados. Os valores da
velocidade horizontal (u) são comparados ao longo da linha vertical no centro da cavidade
(posição x = 0,5m) e os valores da velocidade vertical (v) são comparados ao longo da linha
horizontal no centro da cavidade (posição y = 0,5m). Para testar a sensibilidade do código ao
refino de malha foram utilizadas quatro malhas diferentes para cada situação. O número de
x
v = 0
u=U = 100m/s
u
=
v
=
0
u
=
v
=
0
u
=
v
=
0
u (t = 0) = v (t = 0) = 0
p (t = 0) = 0
y
ρ
= 1 kg/m
3
c = 350 m/s
µ
= variável
λ = 0 Pa⋅s
x=
1
m
y=
1
m
67
elementos em cada direção para cada uma das malhas que foram adotadas são os seguintes:
25x25 para a mais grosseira, 50x50, 100x100 e 200x200 elementos para as malhas mais
refinadas. Todas as malhas eram compostas por elementos quadriláteros (bidimensionais). Como
o programa resolve os escoamentos de forma transiente, a solução apresentada aqui é aquela
obtida com tempo de processamento suficiente para que se atingisse uma solução com o
escoamento desenvolvido no tempo. O passo de tempo utilizado é dependente da malha e era
calculado utilizando-se a seguinte expressão:
U
min
x
t
∆
∆ =
, (6.2)
onde
∆
x
min
é o tamanho do menor elemento da malha. Para visualizar a diferença na resolução
das malhas a Figura 6.2 mostra a malha mais grosseira e a malha mais refinada utilizadas no
trabalho.
Figura 6.2. Malha mais grosseira (25x25 elementos) e malha mais refinada (200x200
elementos) utilizadas no caso da cavidade para comparação.
As simulações foram realizadas utilizando o sistema semi-implícito de avanço no tempo
(com
θ
= 0,5) e utilizando o método de Newton com diagonalização seletiva apenas nos termos
de pressão das matrizes do problema (metodologia proposta) utilizando e=0,7.
Os gráficos com as comparações dos resultados para Re
cavidade
= 100 são mostrados nas
Figuras 6.3 a 6.4. Os valores de posição e velocidade mostrados nesses gráficos, assim como em
todos os outros para o caso, são adimensionais. A adimensionalização é feita utilizando-se o
68
comprimento de referência da cavidade (L = 1m) e a velocidade de referência (U = 100m/s)
segundo as expressões abaixo:
* ; *
L L
* ; *
U U
x y
x y
u v
u v
= =
= =
(6.3)
Figura 6.3. Perfil de velocidade u* (adimensional) ao longo da linha vertical central da
cavidade (em x = 0.5m); comparação com outros resultados publicados. Re
cavidade
= 100.
69
Figura 6.4. Perfil de velocidade v* (adimensional) ao longo da linha vertical central da
cavidade (em y = 0.5m); comparação com outros resultados publicados. Re
cavidade
= 100.
Nestes primeiros resultados exibidos percebe-se a excelente concordância dos perfis
obtidos com aqueles publicados por Ghia et al., 1982. No perfil de velocidade v* se percebe um
pequeno desvio em relação aos valores de referência para a malha mais grosseira, o que não se
repete nas outras.
Exceto para a malha de 25x25 elementos, os resultados de todos os casos rodados com o
código desenvolvido foram coincidentes. Isso mostra que para esse caso (Re
cavidade
= 100) a
malha de 50x50 elementos já fornece um resultado independente de malha.
Nas Figuras 6.5 a 6.8 são mostrados os resultados obtidos para os valores de Re
cavidade
de
400 e 1000.
70
Figura 6.5. Perfil de velocidade u* (adimensional) ao longo da linha vertical central da
cavidade (em x = 0.5m); comparação com outros resultados publicados. Re
cavidade
= 400.
Figura 6.6. Perfil de velocidade v* (adimensional) ao longo da linha vertical central da
cavidade (em y = 0.5m); comparação com outros resultados publicados. Re
cavidade
= 400.
71
Figura 6.7. Perfil de velocidade u* (adimensional) ao longo da linha vertical central da
cavidade (em x = 0.5m); comparação com outros resultados publicados. Re
cavidade
= 1000.
Figura 6.8. Perfil de velocidade v* (adimensional) ao longo da linha vertical central da
cavidade (em y = 0.5m); comparação com outros resultados publicados. Re
cavidade
= 1000.
Esses resultados apresentam uma boa concordância com os resultados de referência. Os
valores obtidos para as malhas de 50x50, 100x100 e 200x200 elementos coincidiram com os
resultados de referência. A concordância dos resultados obtidos com aqueles de Ghia et al.,
1982, foi semelhante àquela de Schreiber & Keller, 1983, que utilizou técnicas numéricas de alta
72
ordem (com um método tipo Newton adaptativo para o sistema não-linear e um procedimento de
continuação para varrer uma ampla faixa de números de Reynolds) e melhor do que a
apresentada por Costa et al., 2005, que utilizou uma técnica estabilizada de elementos finitos
(com esquema SUPG) para obtenção dos resultados.
As Figuras 6.9 e 6.10 mostram os campos de velocidade (módulo) e de pressão obtidos na
solução a Re
cavidade
= 1000.
Figura 6.9. Campo de velocidade (módulo) para Re
cavidade
= 1000.
Figura 6.10. Isolinhas de pressão para Re
cavidade
= 1000.
73
No gráfico de isolinhas de pressão observam-se os picos de pressão negativa e positiva
nos dois cantos da cavidade. Pelo formato das linhas observa-se também que o cálculo da
pressão está implementado de forma bastante estável.
Figura 6.11. Linhas de corrente obtidas sobre o campo resolvido para Re
cavidade
= 1000.
Figura 6.12. Linhas de corrente publicadas por Elias et al., 2004, para Re
cavidade
= 1000.
Pela comparação das Figuras 6.11 e 6.12 percebe-se que as linhas de corrente publicadas
em ambos os trabalhos são coerentes. No trabalho de Elias et al., 2004, foi utilizado um esquema
de Newton inexato associado a uma formulação de elementos finitos estabilizados SUPG/PSPG
74
para a resolução das equações. A malha utilizada era composta por 3200 elementos triangulares
(totalizando 1681 nós).
A seguir são apresentados resultados da análise tridimensional e turbulenta do problema
da cavidade.
6.2. Escoamento turbulento em uma cavidade tridimensional
O problema da cavidade tridimensional foi estudado experimentalmente por Prasad et al.,
1989, onde perfis de velocidade e valores médios de flutuações foram medidos através de
anemometria por LASER-Doppler. Todas as medições foram realizadas sobre o escoamento em
cavidades de seção quadrada de lado 1m e com três profundidades diferentes (0,25m, 0,5m e
1m). Os valores do número de Reynolds utilizado foram entre 3200 e 10000, sendo que o
número de Reynolds para este caso dói definido da mesma forma que para a cavidade
bidimensional, pela equação (6.1).
Mais tarde em um novo trabalho, publicado por Zang et al., 1993, o grupo apresentou um
conjunto de resultados de simulações, utilizando simulação de grandes escalas com modelo sub-
malha dinâmico, com o intuito de reproduzir os dados experimentais.
Dos casos apresentados dois foram escolhidos para validar o uso do código desenvolvido
em uma simulação tridimensional turbulenta. Os dois casos foram o caso a Re
cavidade
= 3200 com
razão de aspecto (lado/altura/profundidade) da cavidade igual a 1:1:1 e o caso a
Re
cavidade
= 10000 com razão de aspecto igual a 1:1:0,5.
6.2.1. Escoamento em uma cavidade tridimensional a Re
cavidade
= 3200
A Figura 6.13 mostra um esquema do domínio do problema, assim como condições de
contorno e propriedades do fluido utilizadas.
75
Figura 6.13. Descrição do problema da cavidade 3D com domínio, condições de contorno e
propriedades do fluido.
Como se pode observar o problema é semelhante ao caso da cavidade bidimensional,
exceto pelo acréscimo da componente w da velocidade (para a terceira dimensão). Nas paredes
laterais e na inferior condições de velocidade nula são impostas. Na parede superior é utilizada
uma condição de impermeabilidade e não-deslizamento, ou seja, o fluido acompanha a
velocidade da parede. Um valor constante é atribuído para a viscosidade absoluta de forma a
obter um valor fixo de Re
cavidade
de 3200. Como condição inicial para o problema, foi rodado um
caso bidimensional com o mesmo número de Reynolds e com a mesma malha (no plano x-y),
esse resultado (para um plano) foi copiado para todos os planos da discretização na terceira
dimensão e este campo foi utilizado como condição inicial do problema. Este procedimento tem
duas vantagens: primeiro acelera o desenvolvimento do escoamento (por não partir de um campo
nulo, mas sim um campo já com as principais características do escoamento); e segundo mostra
que o código desenvolvido consegue mostrar as características tridimensionais dos escoamentos
mesmo partindo de condições “bidimensionais”.
Todas as comparações de resultados com os valores publicados são feitos a partir de
medições no plano central da cavidade ilustrado na Figura 6.13.
A malha utilizada para resolução do problema foi semelhante àquela utilizada por Zang et
al., 1993. Era composta por 32x32x32 elementos hexaédricos, não-uniformes, com os menores
parede superior
u = U = 100m/s
v
=
w
= 0
paredes laterais
e inferior
u = v = w = 0
ρ
= 1 kg/m
3
c = 350 m/s
µ
= 0,03125 Pa⋅s
λ = 0 Pa⋅s
x
y
z
plano central
(z=0,5m)
B
=
1m
H
=
1m
D
=
1m
76
elementos junto à parede com uma razão de tamanho ∆x
min
/B (ou ∆y
min
/B) igual a 1,0⋅10
-2
. A
Figura 6.14 mostra a malha utilizada.
Figura 6.14. Vista frontal e tridimensional da malha para o caso da cavidade 3D,
Re
cavidade
= 3200
Novamente as simulações foram realizadas utilizando o sistema semi-implícito de avanço
no tempo (com
θ
= 0,5) e utilizando o método de Newton com diagonalização seletiva apenas
nos termos de pressão das matrizes do problema (metodologia proposta). O passo de tempo para
o caso foi também calculado pela equação (6.2), seu valor foi de 1,0⋅10
-4
s.
Um tempo de simulação suficiente para que se atingisse um escoamento médio permanente
foi utilizado. Quando atingido esse comportamento, as médias foram tomadas por um tempo
igual a três tempos característicos do problema, sendo o tempo característico considerado
calculado conforme a equação (6.4). O valor obtido foi de t
c
= 0.01s.
L
U
c
t
=
(6.4)
onde L é o comprimento de referência da cavidade sendo que, para o caso, L=B=D=1m.
Os primeiros resultados comparados foram as velocidades médias nas direções x e y
(componentes u e v) nas linhas médias do plano central da cavidade, de forma semelhante ao
caso bidimensional. Os resultados de velocidade média são mostrados nas Figuras 6.15 e 6.16.
77
Figura 6.15. Perfil de velocidade u* (adimensional) ao longo da linha vertical central da
cavidade (em x = 0,5) no plano central; comparação com outros resultados publicados.
Re
cavidade
= 3200.
Figura 6.16. Perfil de velocidade v* (adimensional) ao longo da linha horizontal central da
cavidade (em y = 0,5) no plano central; comparação com outros resultados publicados.
Re
cavidade
= 3200.
Analisando os resultados mostrados nas Figuras 6.15 e 6.16 observa-se uma boa
concordância dos resultados obtidos com o código desenvolvido com os resultados publicados.
78
Nos gráficos pode-se perceber a falta de resolução da malha (caracterizada por “pontas” nas
linhas dos perfis de velocidades). Essa falta de resolução é ainda mais evidente quando se faz a
comparação dos perfis de flutuação RMS, mostrados nas Figuras 6.17 e 6.18. As flutuações
RMS são definidas conforme mostrado nas equações (6.5) e (6.6).
2
RMS
10
*=
U
u
u
′
′
(6.5)
2
RMS
10
*=
U
v
v
′
′
(6.6)
onde
u
′
e
v
′
são os valores das flutuações definidas na equação (3.6).
Figura 6.17. Perfil de velocidade u’
RMS
* (adimensional) ao longo da linha vertical central
da cavidade (em x = 0,5) no plano central; comparação com outros resultados publicados.
Re
cavidade
= 3200.
79
Figura 6.18. Perfil de velocidade v’
RMS
* (adimensional) ao longo da linha vertical central
da cavidade (em y = 0,5) no plano central; comparação com outros resultados publicados.
Re
cavidade
= 3200.
Como foi dito anteriormente, na análise das flutuações médias (RMS) fica mais evidente a
falta de resolução de malha para obter-se uma solução adequada do problema. Essa falta de
resolução de malha também é evidente nos resultados publicados por Zang et al., 1993. Esta falta
de resolução, associada a diferença entre os modelos sub-malha utilizados no trabalho atual e no
publicado por Zang et al., 1993, são os motivos das discrepâncias entre os dois resultados.
Mesmo assim os gráficos mostram uma excelente concordância qualitativa, e uma boa
concordância quantitativa quando comparados aos resultados experimentais.
A Figura 6.19 mostra o campo de velocidade média (magnitude do vetor de velocidade) no
plano central da cavidade. A Figura 6.20 mostra isolinhas de pressão média também no mesmo
plano central e a Figura 6.21 mostra algumas linhas de corrente exibindo a característica de
tridimensionalidade do escoamento.
80
Figura 6.19. Campo de velocidade média (magnitude do vetor) no plano central da
cavidade (z = 0,5m). Re
cavidade
= 3200.
Figura 6.20. Isolinhas de pressão estática média no plano central da cavidade, com escala
de valores reduzida entre -1,0⋅10
-4
Pa e 1,0⋅10
-4
Pa. Re
cavidade
= 3200.
81
Figura 6.21. Linha de corrente partindo do ponto central da cavidade colorida pelo módulo
da velocidade média (local). Re
cavidade
= 3200.
As Figuras 6.19 e 6.20 mostram claramente que uma instabilidade numérica apareceu nos
cantos superiores da cavidade. Esta instabilidade numérica tem algumas origens bem definidas,
conforme discutido no trabalho e na literatura. Uma das fontes dessa instabilidade é a condição
de contorno imposta no problema. Esse fenômeno é bem conhecido e tende a aparecer quando se
usa uma condição de cavidade totalmente fechada. Nessa condição, nos cantos superiores da
cavidade, a velocidade imposta pelas condições de contorno varia de zero ao valor imposto na
parede superior em um único elemento, situação bastante difícil de estabilizar. Outra fonte da
instabilidade é a baixa resolução utilizada nessa malha. Pode-se também pensar numa
instabilidade no esquema de avanço no tempo implementado, porém são necessários mais alguns
testes e comparações para chegar-se a uma conclusão.
Ainda com o intuito de exibir a tridimensionalidade do escoamento no interior da cavidade
a Figura 6.22 mostra isosuperfícies de pressão, nas quais percebe-se tanto a tridimensionalidade
do escoamento quanto as instabilidades no campo de pressão.
82
Figura 6.22. Isosuperfícies de pressão dentro da cavidade. Pressão estática com escala
reduzida entre -1,0⋅10
-4
Pa e 1,0⋅10
-4
Pa. Re
cavidade
= 3200.
6.2.2. Escoamento em uma cavidade tridimensional a Re
cavidade
= 10000
Para o caso do escoamento da cavidade bidimensional, foram feitas três alterações no
esquema apresentado na seção anterior. Primeiro o domínio mostrado na Figura 6.13 é reduzido
à metade na direção z, ou seja, a razão de aspecto lado/altura/profundidade passa a ser de 1:1:0,5.
Os valores de B, D e H conforme definidos na Figura 6.13 passam a ser 1m, 1m e 0,5m
respectivamente. O plano central da cavidade, que no caso anterior estava localizado em z =
0,5m, passa para a cota z = 0,25m.
A segunda alteração é no valor da viscosidade absoluta do fluido. Esse valor é alterado de
forma a obter um número de Reynolds, conforme definido na equação (6.1) igual a 10000. O
valor da
viscosidade absoluta para tanto é então de
µ
= 0,01 Pa⋅s.
A última alteração é na malha utilizada. Novamente foi utilizada uma malha semelhante àquela
apresentada no artigo de Zang et al., 1993. A malha para este caso com Re
cavidade
= 10000 é composta por
64x64x32 elementos
hexaédricos,
não-uniformes,
com os menores elementos junto à parede com
uma razão de tamanho ∆x
min
/B (ou ∆y
min
/B) igual a 5,0⋅10
-3
. A Figura 6.23 mostra a malha
utilizada.
83
Figura 6.23. Isosuperfícies de pressão dentro da cavidade. Pressão estática com escala
reduzida entre -1,0⋅10
-4
Pa e 1,0⋅10
-4
Pa. Re
cavidade
= 10000.
Assim como no caso a Re
cavidade
= 3200 as simulações foram feitas com o sistema semi-
implícito de avanço no tempo (com
θ
= 0,5) e utilizando o método de Newton com
diagonalização seletiva apenas nos termos de pressão das matrizes do problema. O passo de
tempo utilizado foi de 5,0⋅10
-5
s, obtido pela equação (6.2). Aqui também foi simulado um tempo
longo o suficiente para que se atingisse um comportamento médio permanente. Após atingir-se
esse estágio as médias foram tomadas por um tempo igual a três tempos característicos do
problema, com o tempo característico definido novamente pela (6.4) de valor t
c
= 0.01s.
Os resultados dos perfis médios de velocidade são mostrados nas Figuras 6.24 e 6.25.
84
Figura 6.24. Perfil de velocidade u* (adimensional) ao longo da linha vertical central da
cavidade (em x = 0,5) no plano central; comparação com outros resultados publicados.
Re
cavidade
= 10000.
Figura 6.25. Perfil de velocidade v* (adimensional) ao longo da linha horizontal central da
cavidade (em y = 0,5) no plano central; comparação com outros resultados publicados.
Re
cavidade
= 10000.
85
Figura 6.26. Perfil de velocidade u’
RMS
* (adimensional) ao longo da linha vertical central
da cavidade (em x = 0,5) no plano central; comparação com outros resultados publicados.
Re
cavidade
= 10000.
Figura 6.27. Perfil de velocidade v’
RMS
* (adimensional) ao longo da linha vertical central
da cavidade (em y = 0,5) no plano central; comparação com outros resultados publicados.
Re
cavidade
= 10000.
Para este caso um pequeno desvio entre os resultados obtidos com o código desenvolvido e
os resultados de Zang et al., 1993, e de Prasad et al., 1989, é percebido nos perfis de velocidade
86
média. Qualitativamente os perfis ficaram bastante próximos dos resultados de simulação
publicados.
Nos perfis de flutuações médias novamente percebe-se uma boa tendência qualitativa nos
valores obtidos e uma concordância com os resultados experimentais, principalmente nas zonas
onde a malha é mais refinada (junto às paredes). As discrepâncias apresentadas pelo resultado
obtido no presente trabalho novamente podem ser justificadas pela baixa resolução da malha
utilizada e pelo tipo do modelo sub-malha utilizado. Como explicado no texto, a principal
característica do modelo dinâmico é adaptação local às características da turbulência no
escoamento. Assim, para uma mesma malha, espera-se sempre que os resultados obtidos com
esse modelo sub-malha sejam superiores àqueles obtidos com o modelo de Smagorinsky. A boa
concordância, em termos qualitativos e de ordens de grandeza, é extremamente satisfatória dadas
as condições o problema.
Ressalta-se novamente ainda a incerteza quanto à estabilidade no esquema de avanço no
tempo implementado no código desenvolvido e a instabilidade causada pela condição de
contorno do problema, ambos já discutidos e que podem ser motivo de futuros trabalhos.
As Figuras 6.28 a 6.31 mostram o campo de velocidade média (módulo) e linhas de presão
constante (ambos no plano central da cavidade), uma linha de corrente partindo do ponto central
da cavidade (mostrando a tridimensionalidade do escoamento) e superfícies de pressão constante
dentro da cavidade.
Figura 6.28. Campo de velocidade média (magnitude do vetor) no plano central da
cavidade (z = 0,5m). Re
cavidade
= 10000.
87
Figura 6.29. Isolinhas de pressão estática média no plano central da cavidade, com escala
de valores reduzida entre -2,0⋅10
-4
Pa e 2,0⋅10
-4
Pa. Re
cavidade
= 10000.
Figura 6.30. Linha de corrente partindo do ponto central da cavidade colorida pelo módulo
da velocidade média (local). Re
cavidade
= 10000.
88
Figura 6.31. Isosuperfícies de pressão dentro da cavidade. Pressão estática com escala
reduzida entre -2,0⋅10
-4
Pa e 2,0⋅10
-4
Pa. Re
cavidade
= 10000.
Das Figuras 6.30 e 6.31 fica evidente o comportamento tridimensional do problema. Da
figura 6.31 percebe-se também um comportamento simétrico (na direção z) deste campo médio
de pressões.
Para ambos os casos tridimensionais da cavidade apresentados aqui, o tempo de
processamento ainda era muito alto (apesar dos métodos implementados para melhorar o
desempenho numérico). Cada simulação levava em torno de 15 a 20 dias na estação Altix
disponibilizada para fornecer um resultado com tempo de simulação adequado.
Os principais motivos para esta demora foram o pequeno passo de tempo utilizado nas
simulações e a forma com que o código calcula as matrizes de elemento e monta as matrizes
globais (rotinas onde maior parte do tempo de processamento era consumido). O primeiro não
pode ser alterado, por se trata de um parâmetro físico dos escoamentos em questão, o segundo
pode ser resolvido alterando a forma com que essas matrizes são calculadas, melhorando o seu
desempenho. Essa tarefa é simplificada pela estrutura do código, já discutida.
6.3. Escoamento em um degrau bidimensional
O escoamento sobre um degrau é também um caso de referência utilizado em inúmeros
trabalhos. Este escoamento apresenta características bastante interessantes como descolamento e
recolamento de camada limite, formação de uma camada de mistura e, quando em regime
89
turbulento, formação de estruturas já bastante conhecidas devido à diversas investigações
(numéricas e experimentais) disponíveis na literatura.
Ainda para a validação do código desenvolvido neste trabalho, o caso do escoamento sobre
o degrau, em regimes laminar e turbulento, foi simulado. As simulações foram todas realizadas
sobre um mesmo domínio e uma mesma malha bidimensional. No momento a discussão sobre
simulação bidimensional de grandes escalas é deixada de lado, já que isto já foi discutido
anteriormente.
Um detalhado estudo do caso, para uma ampla faixa de números de Reynolds, foi
publicado por Armaly et al., 1983. Nesse artigo diversas medições experimentais, utilizando a
técnica de LASER-Doppler, são apresentados. Perfis de velocidade e características do
escoamento como pontos de descolamento e recolamento de camada limite foram publicados. No
mesmo trabalho um estudo numérico bidimensional utilizando o método de diferenças finitas
também foi realizado.
A Figura 6.32 mostra esquematicamente o domínio, as condições de contorno e as
propriedades do fluido utilizadas.
H = 1m
h = 0,94m
l = 1m L = 30m
x
R
paredes:
u = v = 0
ρ
= 1 kg/m3
c = 350 m/s
µ
= variável
λ
=
0
perfil parabólico
recirculação
principal
Figura 6.32. Diagrama esquemático do problema do degrau, com dimensões, condições de
contorno e propriedades do fluido.
Na figura acima H é a altura do canal de entrada, h a altura do degrau, l o comprimento da
seção de entrada, L o comprimento do canal após o degrau e x
R
o a distância de recolamento da
recirculação principal.
As condições de contorno são de não-deslizamento nas paredes, um perfil parabólico
imposto na entrada do canal, com velocidade média de 100m/s (visando reproduzir a entrada da
seção experimental do caso publicado) e tração nula na saída do canal.
Para o caso laminar foram simulados 8 casos com números de Reynolds entre 100 e 800
com passo de 100. O objetivo era de reproduzir a curva de ponto de recolamento da recirculação
90
principal publicado por Armaly et al., 1983. Para o primeiro caso foi utilizada uma condição
inicial de velocidades e pressão nula. Para os casos subseqüentes foi utilizado o campo resultante
da simulação com número de Reynolds mais baixo como condição inicial. O número de
Reynolds para o caso do degrau é definido como mostrado na equação (6.7).
U2H
degrau
Re
ρ
µ
= (6.7)
onde U é a velocidade média do perfil parabólico na entrada, e H é a altura do canal de entrada.
A malha utilizada em todos os casos foi a mesma e era composta por 10380 elementos
quadriláteros com tamanho uniforme na região do degrau e um pouco mais grosseira mais
próximo da saída do canal. A região próxima ao canal é mostrada na Figura 6.33.
Figura 6.33. Malha utilizada no caso do degrau bidimensional, região próxima ao degrau.
As simulações foram conduzidas utilizando o sistema semi-implícito de avanço no tempo
(com
θ
= 0,5) e com o método de Newton com diagonalização seletiva apenas nos termos de
pressão das matrizes do problema. O passo de tempo utilizado foi de 1,0⋅10
-4
s, um pouco abaixo
do obtido pela equação (6.2), por motivo de garantir-se a estabilidade para automatização das
rodadas (feita através de scripts). Aqui também foi simulado um tempo longo o suficiente para
que se atingisse um comportamento médio permanente.
O ponto de recolamento da recirculação principal é definido aqui como o ponto em que a
velocidade horizontal é nula a partir dos valores medidos na primeira linha de nós acima da
parede inferior. A Figura 6.34 mostra um campo de velocidade (instantâneo) e a Figura 6.35
mostra o respectivo campo de pressão estática obtidos para o caso com Re
degrau
= 800.
91
Figura 6.34. Campo de velocidades obtido para Re
degrau
= 800.
Figura 6.35. Campo de pressão estática obtido para Re
degrau
= 800.
Com todos os resultados do ponto de recolamento foi montado o gráfico mostrado na
Figura 6.36 com a comparação dos valores obtidos com os publicados por Armaly et al., 1983.
Figura 6.36. Comparação dos pontos de recolamento obtidos com valores publicados por
Armaly et al., 1893.
92
Nas figuras dos campos de velocidades e pressão se percebe a característica oscilatória dos
campos na região de saída do canal e também a recuperação de pressão estática devido à
expansão.
No gráfico a posição do ponto de recolamento é expressa em termos de x
R
/S, onde S é a
razão entre a altura do canal de entrada (H) e a altura do degrau (h). Os valores obtidos
apresentaram uma boa concordância com os valores experimentais e numéricos publicados.
Conforme esperado, para valores de Re
degrau
maiores que 400 os valores previstos pela simulação
começam a se distanciar do valor experimental. Isto acontece porque já a partir desses valores,
efeitos de tridimensionalidade do escoamento começam a ser consideráveis, efeitos estes não
capturados por esta simulação.
Os resultados obtidos com o código desenvolvido foram um pouco melhores do que os
publicados por Armaly et al., 1983, para a região de maiores valores de Reynolds. Isto se deve ao
método de resolução utilizado, mais robusto e preciso do que o método simples de diferenças
finitas utilizado pelo autor.
93
7. CONCLUSÃO
A principal conclusão do trabalho é que o código de dinâmica dos fluidos computacional,
para simulação de grandes escalas de escoamentos bi- e tridimensionais, transientes,
incompressíveis, isotérmicos e turbulentos, que servirá como base de futuros trabalhos do grupo,
foi desenvolvido com sucesso. Os resultados das simulações rodadas com o novo código foram
bastante satisfatórios e validaram a metodologia implementada.
Apesar disso algumas dificuldades se evidenciaram no uso deste código. A principal delas
foi o alto tempo de processamento necessário para casos tridimensionais. Como dito isto é
devido à forma com que são calculadas as matrizes de elemento e montadas as matrizes globais
do problema, o que pode ser melhorado em trabalhos futuros. As instabilidades que surgiram no
caso das simulações a altos valores do número de Reynolds também foram preocupantes, e
parecem estar relacionadas ao esquema de avanço no tempo escolhido, o que pode também ser
estudado adiante.
Mesmo com essas questões de desempenho das rotinas utilizadas, os métodos
implementados para melhorar o desempenho do programa foram bastante eficientes. A redução
de memória necessária utilizando o método de armazenagem compacta de matrizes indexadas
por linhas viabilizou o uso de esquemas semi-implícitos de avanço no tempo e se mostrou muito
mais eficiente do que os outros métodos de armazenagem compacta. Além disso esta
armazenagem tem um desempenho excelente com os métodos de solução de sistemas lineares
implementados.
O processamento em paralelo utilizando a técnica OpenMP também se mostrou bastante
eficaz. Nas estações disponibilizadas pela SGI a distribuição dinâmica de processamento, de
acordo com a disponibilidade de cada processador foi acompanhada e sempre manteve o máximo
número de processos paralelos possíveis. Apesar da técnica ser limitada a máquinas de memória
compartilhada, ela pode ser extendida para uso com máquinas de memória distribuída através de
ferramentas de criação de uma “máquina paralela de memória compartilhada virtual” ou pode
ainda ser utilizada em conjunto com técnicas de processamento entre máquinas de memória
distribuída.
Os seguintes temas são dados como sugestões de próximos trabalhos: melhorias na
estrutura de cálculo das matrizes de elemento e montagem das matrizes globais; estudo e
implementação de outros esquemas de avanço no tempo para aumentar a precisão e a
estabilidade do método; implementação do modelo dinâmico e validação e ampliação da
94
biblioteca de elementos para de forma a permitir que o programa utilize, por exemplo, elementos
triangulares e tetraédricos facilitando a geração de malha em geometrias mais complexas.
Outros modelos de turbulência para as escalas sub-malha da simulação de grandes escalas
também podem ser implementados no código.
95
8. REFERÊNCIAS BIBLIOGRÁFICAS
Armaly, B.F., Durst, F., Pereira, J.C.F., Schönung, B., 1983
.
“Experimental and
Theoretical Investigation of Backward-Facing Step Flow”,
Journal of Fluid Mechanics
, v. 127,
pp. 473-496.
Costa, G.K., Lyra, P.R.M. e Oliveira Lira, C.A.B., 2005. “Numerical simulation of two
dimensional compressible and incompressible flows”,
Journal of the Brazilian Society of
Mechanical Sciences and Engineering
, v. XXVII (4), pp. 372-380.
Elias, R.N., Coutinho, A.L.G.A. e Martins, M.A.D., 2004. “Inexact Newton-type methods
for non-linear problems arising from the SUPG/PSPG solution of steady incompressible Navier
Stokes Equations”,
Journal of the Brazilian Society of Mechanical Sciences and
Engineering
, v. XXVI (3), pp. 330-339.
Ferziger, J., 1993. “Simulation of Complex Turbulent Flows: Recent Advances and
Prospects in Wind Engineering”,
Journal of Wind Engineering and Industrial
Aerodynamics
, vol. 46 & 47, pp. 195-212.
Ferziger, J.H. e Perić, M., 2002. “
Computational methods for fluid dynamics
”, 3ª
edição, Springer.
Findikakis, A.N. e Street, R.L., 1982. “Mathematical Description of Turbulent Flows”,
Journal of Hydraulics Division
, vol. 108, no. HY8, 17265, pp. 887-903.
Findikakis, A.N. e Street, R.L., 1982. “Finite Element Simulation of Stratified Turbulent
Flows”,
Journal of Hydraulics Division
, vol. 108, no. HY8, 17266, pp. 904-920.
Fox, R.W. e McDonald, A.T., 2001. “
Introdução à Mecânica dos Fluidos
”, LTC Editora,
Rio de Janeiro.
Germano, M. Piomelli, U., Moin, P. Cabot, W.H., 1991. “A dynamic sub-grid-scale eddy
viscosity model”,
Physics of Fluids,
A3 (7), 1760-1765.
96
Ghia, U., Ghia, K. N., Shin, C. T., 1982. “High-Re Solutions for Incompressible Flow
Using the Navier-Stokes Equations and a Multi-Grid”.
Journal of Computational Physics
, v.
48, pp. 387-411.
Gresho, P.M. e Sani, R.L., 1987. “On pressure boundary conditions for the incompressible
Navier-Stokes equations”, Finite Elements in Fluids, v. 7, pp. 123 – 157.
Gresho, P.M. and Sani, R.L., 1999. “Incompressible
Flow and the Finite Element
Method
”, v. 1&2, John Wiley and Sons, Ltd.
Hinze, J.O., 1975. “
Turbulence
”, McGraw-Hill, New York.
Hughes, T. J. R., 1987.
“The Finite Element Method”
, Prentice-Hall, New Jersey.
Hussain, A.K.M.F., 1983. “Coherent Structures: Reality and Myth”,
Physics of Fluids
,
vol. 26(10), pp. 2816-2850.
Kawahara, M. e Hirano, H., 1983. “A Finite Element Method for High Reynolds Number
Viscous Fluid Flow Using Two Step Explicit Scheme”,
International Journal For Numerical
Methods in Fluids
, vol. 3, pp. 137-163.
Lange, C.F., 1992. “Simulação de escoamentos incompressíveis não-isotérmicos pelo
método de elementos finitos com função de penalidade”,
Dissertação de Mestrado
, Programa
de Pós-Graduação em Engenharia Mecânica (PROMEC), Universidade Federal do Rio Grande
do Sul (UFRGS), Porto Alegre.
Lee, R.L., Gresho, P.M. e Sani, R.L., 1979. “Smoothing techniques for certain primitive
variable solutions of the Navier-Stokes equations”, International Journal for Numerical Methods
in Engineering, v. 14, pp. 1785-1804.
Leonard, A., 1974. “Energy Cascade in Large-Eddy Simulations of Tubulent Fluid Flows”,
Advances in Geophysics
, vol. 18A, pp. 237-248.
97
Lesieur, M., Comte, P. e Métais, O., 1995. “Numerical Simulations of Coherent Vortices in
Turbulence”,
Applied Mechanics Reviews
, vol. 48(3), pp. 121-149.
Lilly, D.K., 1992. “A proposed moddification of the Germano subgrid-scale closure
method”,
Physics of Fluids,
A4
(3), 633-635.
Nallasamy, M. e Prasad, K.K., 1977. “On cavity flow at high Reynolds numbers”,
Journal
of Fluid Mechanics
, v. 79 (2), pp. 391-414.
Pereira, V.D. e Campos Silva, J. B., 2005. “Simulations if incompressible fluid flows by a
least squares finite element method”,
Journal of the Brazilian Society of Mechanical Sciences
and Engineering
, v. XXVII (3), pp. 274-282.
Petry, A.P. e Awruch, A., 1997. “Simulação de Escoamentos Turbulentos Pelo Método de
Elementos Finitos Através da Simulação Direta de Grandes Vórtices”,
ENIEF
.
Petry, A.P., 2002. “Análise Numérica de Escoamentos Turbulentos Tridimensionais
Empregando o Método de Elementos Finitos e Simulação de Grandes Escalas”,
Tese de
Doutorado
, Programa de Pós-Graduação em Engenharia Mecânica (PROMEC), Universidade
Federal do Rio Grande do Sul (UFRGS), Porto Alegre.
Petry, A.P. e Awruch, A.M., 2004. “Large eddy simulation of three-dimensional turbulent
flows by the finite element method”,
IV Escola de Primavera de Transição e Turbulência
,
Porto Alegre.
Piomelli, U., Moin, P. e Ferziger, J.H., 1988. “Model Consistency in Large Eddy
Simulation of Turbulent Channel Flows”,
Physics of Fluids
, vol. 31(7), pp. 1884-1891.
Piomelli, U., Ferziger, J., Moin, P. e Kim, J., 1989. “New Approximate Boundary
Conditions for Large Eddy Simulations of Wall-Bounded Flows”,
Physics of Fluids
A, vol. 1(6),
pp. 1061-1068.
Piomelli, U., Cabot, W.H., Moin, P. e Lee, S., 1991. “Subgrid-scale Backscatter in
Turbulent and Transitional flows”,
Physics of Fluids
A, vol. 3(7), pp. 1766-1771.
98
Piomelli, U., 1993. “High Reynolds number calculations using the dynamic subgrid-scale
stress model”,
Physics of Fluids
A, vol. 5(6), pp. 1484-1490.
Piomelli, U., 1999. “Large-Eddy Simulation: Achievements and Challenges”,
Progress in
Aerospace Sciences
, vol. 35, pp. 335-362.
Pope, S.B., 2003. “Ten Questions Concerning the Large-Eddy Simulation of Turbulent
Flows”,
New Journal of Physics
, November.
Popiolek, T. L., Awruch, A.M. e Teixeira, P.R.F., 2006. “Finite element analysis of
laminar and turbulent flows using LES and sub-grid scale models”,
Applied Mathematical
Modelling
, v. 30, pp 177-199.
Prasad. A. K., Koseff, J. R., 1989. “Reynolds Number and End-wall Effects on a Liddriven
Cavity Flow”.
Physics of Fluids,
A, v. 1, n. 2, pp 208-218.
Press, S.A., Vetterling, W.T. e Flannery, B.P., 1992. “
Numerical Recipes in Fortran
77: the Art of Scientific Computing
”, Cambridge, Cambridge.
Reddy, J.N., 1985. “
An introduction to the finite element method
”, McGraw-Hill, New
York.
Reddy, J.N. e Gartling, D.K., 1994. “
The Finite Element Method in Heat Transfer and
Fluid Dynamics
”, CRC Press, Boca Raton.
Rogallo, R.S. e Moin P., 1984. “Numerical Simulation of Turbulent Flows”,
Annual
Reviews on Fluid Mechanics
, vol. 16, pp. 99-137.
Schlichting, H., 1979. “
Boundary-layer Theory
”, McGraw-Hill, New York.
Schreiber, R. e Keller, H.B., 1983. “Driven cavity flows by efficient numerical
techniques”,
Journal of Computational Physics
, v. 49, pp. 310-333.
99
Selvam, R.P., 1997. “Finite element modeling of flow around a circular cylinder using
LES”, Journal of Wind Engineering and Industrial Aerodynamics, v. 67&68, pp. 129-139.
Silva Freire, A.P., Menut, P.P.M e Su, J., 2002. “
Turbulência
”, Associação Brasileira de
Ciências Mecânicas (ABCM), Rio de Janeiro.
Silveira Neto, A., Grand, D., Métais, O. e Lesieur, M., 1993. “A Numerical Investigation
of the Coherent Vortices in Turbulence Behind a Backward-Facing Step”, Journal of Fluid
Mechanics, vol. 256, pp. 1-25.
Smagorinsky, J., 1963. “General circulation experiments with the primitive equations. I.
The basic experiment”, Monthly Weather Reviews, v. 91, pp. 99-164.
Sonntag, R.E., Borgnakke, C. e Van Wylen, G.J., 1998. “
Fundamentos da
Termodinâmica
”, Edgard Blücher Ltda., São Paulo.
Tennekes, H. e Lumley, J.L., 1994. “
A First Course in Turbulence
”, MIT Press, Londres.
White, F.M., 1999. “
Mecânica dos Fluidos
”, McGraw-Hill, Rio de Janeiro.
Wilcox, D.C., 1994. “
Turbulence Modeling for CFD
”, DCW Industries Inc. La Cañada,
Califórnia.
Zang, Y., Street, R. L., Koseff, J. R., 1993. “A Dynamic Mixed Subgrid-Scale Model and
its Application to Turbulent Recirculating Flows”.
Physics of Fluids,
A, v. 5, n. 12, pp. 3186-
3196.
Zienkiewicz, O.C. e Taylor, R.L., 2000. “
The Finite Element Method - Volume 3: Fluid
Dynamics
”, Butterworth-Heinemann, Oxford.
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
