Embed Size (px)
Citation preview
UNIVERSIDADE FEDERAL FLUMINENSECTC - Centro TecnológicoTCE - Escola de EngenhariaTEM - Departamento de Engenharia Mecânica
PROJETO DE GRADUAÇÃO II'
&
$
%
Título do Projeto:
ANÁLISE DA INFLUÊNCIA DE MODELOS DETURBULÊNCIA NA FLUIDODINÂMICA DE UM SELO
TIPO HOLE-PATTERN
'
&
$
%
Autor(es):
ALEX DA CUNHA MARROIG
LUCAS DE ABREU NIEDNER NUNES
'
&
$
%
Orientador(es):
JOÃO FELIPE MITRE DE ARAUJO, Ph.D.
,
Data: 9 de Agosto de 2016
ALEX DA CUNHA MARROIG
LUCAS DE ABREU NIEDNER NUNES
ANÁLISE DA INFLUÊNCIA DE MODELOS DETURBULÊNCIA NA FLUIDODINÂMICA DE UM SELO TIPO
HOLE-PATTERN
Trabalho de Conclusão de Curso apresentado aoCurso de Engenharia Mecânica da Universidade Federal Flu-minense, como requisito parcial para obtenção do grau deEngenheiro Mecânico.
Orientador(es):
JOÃO FELIPE MITRE DE ARAUJO, Ph.D.
,
Niterói
9 de Agosto de 2016
Ficha Catalográfica elaborada pela Biblioteca da Escola de Engenharia e Instituto de Computação da UFF
M361 Marroig, Alex da Cunha
Análise da influência de modelos de turbulência na
fluidodinâmica de um selo tipo Hole-Pattern / Alex da Cunha
Marroig, Lucas de Abreu Niedner Nunes. – Niterói, RJ : [s.n.], 2016.
129 f.
Trabalho (Conclusão de Curso) – Departamento de Engenharia
Mecânica, Universidade Federal Fluminense, 2016.
Orientador: João Felipe Mitre de Araujo.
1. Dinâmica dos fluidos (Engenharia). 2. Selo mecânico. 3.
Fluidodinâmica computacional. I. Nunes, Lucas de Abreu Niedner.
Título.
CDD 620.1064
UNIVERSIDADE FEDERAL FLUMINENSECTC - Centro TecnológicoTCE - Escola de EngenhariaTEM - Departamento de Engenharia Mecânica
PROJETO DE GRADUAÇÃO II
AVALIAÇÃO FINAL DO TRABALHO
Título do Trabalho:ANÁLISE DA INFLUÊNCIA DE MODELOS DE TURBULÊNCIA NA
FLUIDODINÂMICA DE UM SELO TIPO HOLE-PATTERN
Parecer do Professor Orientador da Disciplina:
− Grau Final recebido pelos Relatórios de Acompanhamento:
− Grau atribuído ao grupo nos Seminários de Progresso:
Parecer do Professor(es) Orientador(es):
Nome e Assinatura do Professor(es) Orientador(es):
Prof.: João Felipe Mitre de Araujo. Assinatura:
Prof.: . Assinatura:
Parecer Conclusivo da Banca Examinadora do Trabalho:
Projeto Aprovado Sem Restrições
Projeto Aprovado Com Restrições
Prazo concedido para cumprimento das exigências:
Discriminação das exigências e/ou observações adicionais:
UNIVERSIDADE FEDERAL FLUMINENSECTC - Centro TecnológicoTCE - Escola de EngenhariaTEM - Departamento de Engenharia Mecânica
PROJETO DE GRADUAÇÃO II
AVALIAÇÃO FINAL DO TRABALHO(continuação)
Aluno: Alex da Cunha Marroig. Grau:
Aluno: Lucas de Abreu Niedner Nunes. Grau:
Composição da Banca Examinadora:
Prof.: João Felipe Mitre de Araujo, Ph.D. Assinatura:
Prof.: Hugo Alvarenga Oliveira, Ph.D. Assinatura:
Prof.: Roger Matsumoto Moreira, Ph.D. Assinatura:
Local e Data de Defesa do Trabalho:
Departamento de Engenharia Mecânica, 09/08/2016
DEDICATÓRIA
A todas as pessoas que superam as adversidades com um sorriso no rosto e trabalham
todos os dias para tornar nosso mundo um lugar melhor.
6
AGRADECIMENTOS
Primeiramente agradeço à Deus, pois Ele possibilitou o acontecimento de tudo que foi
feito até hoje.
À minha mãe e meu pai pelo apoio em todos os momentos da minha vida e pelas sábias
orientações que sobre as coisas mais importantes da vida e que não há dever maior quanto o
dever de ser feliz. Obrigado pela compreensão, pela paciência nas horas complicadas, pelo
amor, e por todos os valores que me ensinaram que levarei para sempre.
Agradeço aos meus amigos pela eterna e inigualável amizade que me acompanharam
durante todos estes anos. Aos meus colegas de curso, aos amigos que fiz no PET-MEC,
obrigado por fazer meu dias na faculdade mais felizes. Agradeço principalmente ao meu
amigo Lucas Niedner, que me acompanhou desde início da minha trajetória na faculdade.
Obrigado por me proporcionar sua companhia, sempre bem humorada, de caráter único e por
todos momentos de superação, alegria e vitória.
Agradeço aos mestres que passaram não só apenas conhecimento, mas sabedoria para
lidar com todos percalços desta vida. Obrigado a Universidade Federal Fluminense, por ter
me proporcionado oportunidades únicas de amadurecimento e desenvolvimento. Agradeço
especialmente ao Professor João Felipe Mitre e a Professora Fabiana R. Leta pelos conselho,
orientações e por ter acredito na minha capacidade para me desenvolver nessa trajetória.
Alex da Cunha Marroig
7
8
Agradeço primeiramente aos meus pais, Adilson e Cecile, por terem sempre feito questão
de me proporcionar educação da melhor qualidade possível, dentro de casa e fora dela, o
que foi primordial para que eu me tornasse a pessoa que sou hoje e para que chegasse onde
cheguei.
Agradeço a meu irmão Paulo e demais familiares por sempre acreditarem em mim e me
incentivarem a ser uma pessoa cada vez melhor.
Agradeço a todos os amigos que fiz durante minha trajetória na faculdade e por todos os
momentos de diversão e descontração que me proporcionaram, momentos estes que ajudaram
bastante a enfrentar todas as dificuldades e pressões da graduação com alegria e confiança.
Aos meus colegas de turma, aos amigos da Equipe Buffalo com quem passei incontáveis
fins de semana trabalhando e a todas as outras grandes amizades que construí. Agradeço em
especial ao meu amigo Alex Marroig pela amizade que começou desde o primeiro dia de trote
e se estendeu ao longo de todos esses anos, obrigado por confiar em mim e fazer este trabalho
comigo.
Agradeço também aos bons professores que tive nesta universidade, por muitas vezes
terem aberto meus olhos e minha mente para novos conhecimentos e por exercerem uma
profissão tão nobre a qual, infelizmente, não é dada a devida importância em nosso país.
Agradecimento em especial ao professor orientador João Felipe Mitre por ter confiado em
nossa capacidade para realizar este trabalho e por ter sempre se colocado à disposição de nos
ajudar.
Por fim, agradeço à Universidade Federal Fluminense por ter me acolhido tão bem, tendo
se tornado minha segunda casa (durante alguns anos posso dizer até que se tornou minha
primeira casa) e por ter me proporcionado os melhores anos da minha vida até agora.
Lucas de Abreu Niedner Nunes
RESUMO
As máquinas de fluxo estão presentes em boa parte da indústria de produção de energia
por apresentarem inúmeras aplicações, fornecerem grande potência e por serem eficientes. A
importância destas máquinas neste ramo da indústria tem motivado diversos estudos na área
de dinâmica dos fluidos tanto para as máquinas em si quanto para seus vários componentes.
Com os recentes avanços na área de modelagem computacional e simulação numérica, muitos
destes estudos têm sido realizados com auxílio da fluidodinâmica computacional, ou CFD (do
inglês Computational Fluid Dynamics). O presente trabalho se propõe a analisar os efeitos de
quatro diferentes modelos de turbulência (κ−ε, κ−ω, SST e RNG κ−ε) na fluidodinâmica
de um escoamento através de uma seção angular de um selo mecânico tipo hole-pattern,
um dos principais componentes de uma turbomáquina. Para tal análise, este escoamento é
simulado para cada modelo de turbulência utilizando o software comercial ANSYS CFX®,
a partir do qual são obtidos dados de distribuição de pressão, tensões cisalhantes e vazão
mássica. Os dados obtidos de distribuição de pressão e tensões cisalhantes são comparados
entre os modelos enquanto os dados obtidos para vazão mássica são comparados a um valor
experimental. Os resultados deste trabalho demonstram que os quatro modelos analisados
possuem praticamente a mesma eficiência na determinação da distribuição de pressão e das
tensões cisalhantes, tanto para uma malha mais grossa quanto para uma malha mais refinada,
porém na determinação da vazão mássica o modelo SST apresentou valores melhores quando
comparados ao valor experimental.
Palavras-Chave: Modelo de Turbulência. Fluidodinâmica Computacional. Selo Hole-pattern.
rotodinâmica
9
ABSTRACT
The turbomachines are very relevant to the energy production industry due to their big
number of applications, for being a great power source and having good efficiency. The
importance of these machines in this sector of the industry have been motivating several
studies in the area of fluid dynamics for both the machines themselves and their components.
With the recent improvements in the area of computational modeling and numerical simulation,
most of these studies are being done with the aid of computational fluid dynamics. The purpose
of this study is to analyse the effects of four different turbulence models (κ−ε, κ−ω, SST e
RNG κ−ε) on the fluid dynamics of a flow through an angular section of a hole-pattern gas
seal, one of the main components of a turbomachine. For this analysis, the flow is simulated for
each turbulence model using the comercial software ANSYS CFX®, from which is obtained
data of pressure distribution, shear stresses and mass flow. The obtained data for pressure
distribution and shear stress are compared among the models while the data for mass flow
is compared to an experimental value. The results of this study show that the four analysed
models have similar efficiency when determining pressure distribution and shear stresses, for
both coarse and fine meshes, however, when determining the mass flow rate the SST model
showed better results when compared to the experimental value.
Key-Words: Turbulence models. Computer Fluid Dynamics. Hole-Pattern Seal. Rotordynamic
10
LISTA DE FIGURAS
2.1 Corte de um Compressor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Corte de um selo mecânico tipo labirinto . . . . . . . . . . . . . . . . . . . . . 5
2.3 Modelo simplificado de um selo mecânico com indicações dos componentes . 6
2.4 Modelo simplificado dos selos planos cilíndricos, cônicos e escalonados, da
esquerda para a direita respectivamente . . . . . . . . . . . . . . . . . . . . . . 7
2.5 Distribuição de pressão num selo devido ao efeito lomakin . . . . . . . . . . . 7
2.6 Abordagem física do caso em que existe excentricidade entre o eixo do rotor
e do estator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.7 Simplificação do caso transiente com a mudança de referencial . . . . . . . . . 11
2.8 Possíveis disposições das lâminas do selo labirinto . . . . . . . . . . . . . . . . 13
2.9 Configuração do selo colméia . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.10 Validade da hipótese do contínuo . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.11 Balanço de massa do elemento de um fluido . . . . . . . . . . . . . . . . . . . 23
2.12 Componentes das tensões na superfície de controle . . . . . . . . . . . . . . . 26
2.13 Componentes das tensões na direção x . . . . . . . . . . . . . . . . . . . . . . . 27
2.14 Fluxograma de modelagem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.15 Abordagens das simulações de escoamentos turbulentos disponíveis . . . . . . 33
2.16 Medição de velocidade de um escoamento turbulento . . . . . . . . . . . . . . 35
2.17 Distribuições típicas de velocidade e tensões cisalhantes próximas à parede . 39
2.18 Verificação experimental das leis das camadas interna, intermediária e externa 40
2.19 Uso de funções de parede . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.20 Representação do cálculo usando volumes de controle . . . . . . . . . . . . . . 53
2.21 Oscilações numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
2.22 Gráfico sobre a suavização dos altos gradientes . . . . . . . . . . . . . . . . . . 55
2.23 Subdivisão de um domínio em elementos de malha . . . . . . . . . . . . . . . 58
2.24 Malha refinada nas paredes para captação de efeitos específicos . . . . . . . . 59
2.25 Malha híbrida na cauda de um aerofólio . . . . . . . . . . . . . . . . . . . . . . 61
11
12
2.26 Exemplos dos principais parâmetros de malha . . . . . . . . . . . . . . . . . . 62
2.27 Tipos usuais de malha estruturada . . . . . . . . . . . . . . . . . . . . . . . . . 63
2.28 Extrapolação usada no método Multigrid . . . . . . . . . . . . . . . . . . . . . 65
3.1 Geometria do domínio de estudo e trecho de 6.2º graus da geometria . . . . . 69
3.2 Demonstração do escoamento no domínio analisado . . . . . . . . . . . . . . . 69
3.3 Visualização da malha próxima à parede - detalhe 1 . . . . . . . . . . . . . . . 75
3.4 Visualização da malha nos centros dos cilindros - detalhe 2 . . . . . . . . . . . 75
3.5 Configuração do Regime Estacionário . . . . . . . . . . . . . . . . . . . . . . . 77
3.6 Configuração do Fluido de trabalho . . . . . . . . . . . . . . . . . . . . . . . . 77
3.7 Configuração do modelo de turbulência e energia . . . . . . . . . . . . . . . . 78
3.8 Configuração da Condições de Contorno na Entrada . . . . . . . . . . . . . . . 79
3.9 Configuração da condições de contorno na saída . . . . . . . . . . . . . . . . . 79
3.10 Configuração da região do estator . . . . . . . . . . . . . . . . . . . . . . . . . 80
3.11 Configuração das condições de contorno no rotor . . . . . . . . . . . . . . . . 80
3.12 Configuração das condições de contorno da interface lateral da geometria . . 81
3.13 Configuração do controle de processamento - Solver Control . . . . . . . . . . 82
3.14 Monitoração dos parâmetros físicos do solver . . . . . . . . . . . . . . . . . . . 83
3.15 Planos transversais na geometria . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.1 Residuais do modelo κ−ε com equação da energia completa . . . . . . . . . . 86
4.2 Variação total de temperatura do modelo κ−ε com equação da energia completa 87
4.3 Visualização das linhas de correntes no modelo κ−ε . . . . . . . . . . . . . . 88
4.4 Visualização Y-Plus em ambas malhas M1 e M2 . . . . . . . . . . . . . . . . . 88
4.5 Erro numérico apresentado pelo modelo κ−ε na região de entrada na malha M1 89
4.6 Visualização dos vetores de velocidade na parte lateral do selo - resultado
obtido a partir do modelo SST na malha M2 . . . . . . . . . . . . . . . . . . . 90
4.7 Visualização dos vetores de velocidade na parte lateral do selo - resultado
obtido a partir do modelo SST na malha M2 . . . . . . . . . . . . . . . . . . . 90
4.8 Residuais obtidos utilizando o modelo κ−ε em ambas malhas M1 e M2 . . . 92
13
4.9 Residuais obtidos utilizando o modelo κ−ω em ambas malhas M1 e M2 . . . 93
4.10 Residuais obtidos utilizando o modelo RNG κ−ε em ambas malhas M1 e M2 94
4.11 Residuais obtidos utilizando o modelo SST em ambas malhas M1 e M2 . . . 95
4.12 Visualização das tensões cisalhantes na parede do rotor do modelo κ−ε em
ambas malhas M1 e M2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.13 Visualização das tensões cisalhantes na parede do rotor do modelo κ−ω em
ambas malhas M1 e M2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.14 Visualização das tensões cisalhantes na parede do rotor do modelo RNG κ−εem ambas malhas M1 e M2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.15 Visualização das tensões cisalhantes na parede do rotor do modelo SST em
ambas malhas M1 e M2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.16 Visualização das escalas de parede do modelo SST em ambas malhas M1 e M2 98
4.17 Visualização do uso do primeiro blend SST em ambas malhas M1 e M2 . . . 98
4.18 Visualização do uso do segundo blend SST em ambas malhas M1 e M2 . . . 99
4.19 Gráfico das pressões totais média obtidas nos planos transversais para a malha
M1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4.20 Gráfico das pressões totais média obtidas nos planos transversais para a malha
M2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
LISTA DE TABELAS
2.1 Balanço de fluxos no volume de controle elementar 2.11 . . . . . . . . . . . . 25
2.2 Valores das propriedades φ,Γφ,Sφ . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.1 Geometria do Selo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.2 Configurações das Condições de Contorno . . . . . . . . . . . . . . . . . . . . 70
3.3 Discretização do domínio computacional . . . . . . . . . . . . . . . . . . . . . 74
4.1 Tabela dos parâmetros da resolução dos modelos de turbulência para malha M1 91
4.2 Tabela dos parâmetros da resolução dos modelos de turbulência para malha M2 91
4.3 Tabela das vazões obtidas em cada simulação para malha M1 . . . . . . . . . 100
4.4 Tabela das vazões obtidas em cada simulação para malha M2 . . . . . . . . . 100
4.5 Tensões Cisalhantes no rotor para a malha M1 . . . . . . . . . . . . . . . . . . 101
4.6 Tensões Cisalhantes no rotor para a malha M2 . . . . . . . . . . . . . . . . . . 101
6.1 Configuração dos planos transversais para análise de pressão média . . . . . . 105
14
NOMENCLATURA
x, y, z coordenadas cartesianas
u, v, w componentes cartesianas da velocidade
y+ distância normal em coordenadas de parede
E energia
κ energia cinética turbulenta
L comprimento
Cv calor específico à volume constante
Cp calor específico à pressão constante
M matriz massa
K matriz rigidez ou rigidez direta
k rigidez cruzada
C matriz de amortecimento
x, x, x vetores deslocamento, velocidade e aceleração
Símbolos Gregos
ρ massa específica
ν viscosidade cinemática
ϑ constante de linearização
δ distancia infinitesimal
µ viscosidade dinâmica
ε taxa de dissipação turbulenta
θ temperatura
τ tensão cisalhante
Ω velocidade do eixo central em torno do estator
Subscritos
i , j índices de coordenadas cartesianas
15
16
p pressão constante
up upstream
i p integration point
t tempo
Números adimensionais
Re número de Reynolds
Pe número de Peclet
K u número de knudsen
Pr número de Prandtl
Sc número de Schmidt
Siglas
CFD Computational Fluid Dynamics
DNS Direct Numerical Simulation
LES Large Eddy Simulation
RANS Reynolds Averaged Navier-Stokes
SST Shear Stress Tensor
RNG Re-Normalisation Group
RSM Reynolds Stress Model
MEF Método de Elementos Finitos
MDF Método de Diferenças Finitas
MVF Método de Volumes Finitos
EbFVM Element-based Control Volume Methods
CDS Central Differencing Scheme
UDS Upwind Differencing Scheme
AMG Algebric Multigrid
SUMÁRIO
1. Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Motivação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Objetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Agradecimentos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2. Revisão Bibliográfica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1 Turbomáquinas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Tipos de Selo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2.1 Selo Plano ou Liso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2.2 Selo de Anéis Flutuantes . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.3 Selo de Contato . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.4 Selo Labirinto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.5 Selo Honeycomb . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.6 Selo Hole-pattern . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3 Hipótese do Contínuo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4 Números Admensionais Relevantes . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4.1 Número de Reynolds . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4.2 Número de Prandtl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4.3 Número de Schmidt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.4.4 Número de Péclet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.5 Pequenas Escalas de Kolmogorov . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.6 Princípios de Mecânica dos Fluidos . . . . . . . . . . . . . . . . . . . . . . . . 22
2.6.1 Lei de Conservação da Massa . . . . . . . . . . . . . . . . . . . . . . . 23
2.6.2 Conservação da Quantidade de Movimento . . . . . . . . . . . . . . . 25
2.6.3 Conservação da Energia . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.7 Turbulência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.7.1 Tratamento da Turbulência . . . . . . . . . . . . . . . . . . . . . . . . . 32
17
18
2.7.2 Lei da Parede . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.7.3 Funções de Parede . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.8 Modelos de Turbulência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.8.1 Modelos de Equação Zero . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.8.2 Modelos de Uma equação . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.8.3 Modelos de Duas Equações . . . . . . . . . . . . . . . . . . . . . . . . 43
2.9 Modelagem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.9.1 Métodos de Discretização na Análise Numérica Computacional . . . 49
2.9.2 Malha . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
2.9.3 Tipos de Malha . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
2.9.4 Classificação de Elementos . . . . . . . . . . . . . . . . . . . . . . . . 61
2.9.5 Estratégias para Uso de Malha Estruturada . . . . . . . . . . . . . . . . 63
2.9.6 Aspecto Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
2.9.7 Análise de Erros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3. Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.1 Geometria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.2 Condições de Contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.3 Pré-Análise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.4 Discretização do Domínio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.5 Pré-Processamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.5.1 Condições de Contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.5.2 Controle do Processamento . . . . . . . . . . . . . . . . . . . . . . . . 81
4. Resultados e Discussão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.0.1 Verificação e Validação Preliminar . . . . . . . . . . . . . . . . . . . . 85
4.0.2 Resultados dos Modelos de Turbulência . . . . . . . . . . . . . . . . . 87
4.0.3 Análise dos Residuais . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.0.4 Tensões Cisalhantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
4.0.5 Comparações entre os Modelos de Turbulência . . . . . . . . . . . . . 99
19
5. Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.1 Sugestões para Trabalhos Futuros . . . . . . . . . . . . . . . . . . . . . . . . . 103
6. Apêndice A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
7. Referências Bibliográficas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
1 – Introdução
As turbomáquinas ou máquinas de fluxo estão presentes em diversos setores da engenharia,
sendo componentes críticos por fornecer grande potência e ter boa eficiência. Além de serem
fundamentais para qualquer planta de processo, podem ser usadas em diferentes aplicações
industriais que envolvem a troca de energia mecânica entre um sistema mecânico e um fluido
sendo então encontradas nas indústrias militar, aeronáutica, aeroespacial, automotiva, naval e
de geração de energia com alta eficiência como por exemplo, hidrelétricas e termoelétricas.
Dentre os componentes essenciais para o funcionamento das turbomáquinas, pode-se destacar
o selo mecânico, cuja função é garantir uma vedação apropriada entre o eixo do rotor e
a carcaça, influenciando também nas características rotodinâmicas da máquina [1]. O selo
mecânico desempenha um papel fundamental para a determinação destes coeficientes, atuando
não só para a vedação do mesmo, quanto para a redução de vibrações que causam redução da
vida útil do maquinário e da eficiência.
O estudo do impacto do selo mecânico na vibração de turbomáquinas foi iniciado por
Lomakin [2] em 1958, no qual analisou-se o desbalanceamento causado por um fluxo radial
em uma passagem anelar deste selo e posteriormente foram elaborados trabalhos onde foram
considerados o deslocamento do eixo não ocorrendo ao redor de um ponto de equilíbrio [3],
os efeitos da pressão em selos longos considerando pequenos movimentos ao redor da posição
centrada [4], a turbulência para o desenvolvimento de um modelo analítico, a partir do método
das pertubações baseado nas equações de Hirs Lubrication Equation para determinação de
coeficientes de inércia, amortecimento e rigidez, sendo conhecido como método bulk-flow [5].
A determinação destes coeficientes e a melhor compreensão dos fenômenos de turbulência
associados a este fluxo radial vêm sendo estudada e aprimorada ao longo dos anos tornando-se
fundamental para o desenvolvimento de maquinários mais eficientes e potentes.
1
2
1.1 Motivação
A competição entre as empresas pelo mercado têm levado a busca por melhorias em seus
produtos, na qual faz-se necessário o uso de ferramentas computacionais para se otimizar
o tempo e obter equipamentos mais eficientes. Assim, o estudo dos modelos de turbulência
oferece a possibilidade da obtenção de resultados numéricos mais precisos, com um custo
computacional menor, tornando os processos de desenvolvimento e de criação muito mais
rentáveis e viáveis.
Os trabalhos anteriores não consideraram os diversos modelos de turbulência que poderiam
fornecer resultados mais precisos, assim, a motivação deste trabalho é a escolha do modelo de
turbulência mais adequado para obtenção de melhores resultados para a análise dos efeitos
rotodinâmicos influenciados pelo selo mecânico em questão.
1.2 Objetivo
O objetivo deste trabalho é compreender a influência do fenômeno da turbulência por
meio dos modelos de turbulência nas características de rotodinâmica originados a partir de
um escoamento através de uma folga em um selo mecânico tipo Hole-Pattern por meio do
uso da fluidodinâmica computacional.
Assim, é desenvolvido um modelo matemático baseado em fluidodinâmica computacional
(Computer Fluid Dynamics, CFD) que possa capturar de uma maneira suficientemente precisa,
os efeitos de turbulência para o problema de interesse.
É utilizado para a modelagem da geometria referente ao espaço entre rotor e o estator,
um programa de Computer Aided Design, o Design Modeler®. Para a geração de malhas e
domínio da solução do problema, os programas comerciais Pointwise® 17R4, para a geração
de malhas, e ANSYS CFX®, para a resolução e visualização dos resultados, são utilizados.
O modelo matemático usado se baseia no conceito da média de Reynolds (Reynolds-
Averaged Navier-Stokes, RANS) adicionado de um modelo de turbulência de duas equações
para capturar os efeitos da tensão de cisalhamento relativo às forças viscosas e à vazão obtida
no selo tipo Hole-Pattern.
3
1.3 Agradecimentos
Gostaríamos de agradecer à PETROBRAS por ceder a licença de uso dos programas
computacionais no âmbito do projeto "‘Desenvolvimento de metodologias e ferramentas
numéricas para a obtenção de coeficientes dinâmicos de selos internos de compressores
centrífugos"’ (processo 2012/00251-0). Também agradecemos a este projeto por nos ter
fornecido as malhas utilizadas neste estudo.
2 – Revisão Bibliográfica
2.1 Turbomáquinas
As turbomáquinas podem ser classificadas simplesmente por fornecer energia como
turbinas ou consumir energia como compressores e bombas [6]. A Figura 2.1 mostra o corte
de um compressor.
Fig. 2.1: Corte de um Compressor
Fonte:Retirado de Gonçalves [6]
Para se otimizar uma equipamento mecânico, deve-se analisar os componentes críticos do
mesmo em busca de uma solução mais inteligente que vise o aperfeiçoamento do equipamento.
Assim, dentre os diversos componentes das turbomáquinas há o selo de mecânico, mostrado
na Figura 2.2, cujas características são objeto de estudo deste trabalho.
Sua principal função é garantir uma vedação apropriada entre o eixo do rotor e a carcaça,
influenciando também nas características rotodinâmicas da máquina [1].
4
5
Fig. 2.2: Corte de um selo mecânico tipo labirinto
Fonte: Retirado do site www.eagleburgmann.com [7]
O selo mecânico pode ser utilizado em equipamentos rotativos como bombas, com-
pressores, misturadores e turbomáquinas em geral. Assim como na aplicação em outros
equipamentos, os selos mecânicos são essenciais para o funcionamento de máquinas de fluxo
de grande e médio porte, pois atuam minimizando o fluxo indesejável (vazamento) do fluido
de trabalho que flui da zona de maior pressão para zona de menor pressão, sendo emitido para
o meio externo (atmosfera) [8].
Como o fluido bombeado por estes equipamentos rotativos tende a ocupar todos os espaços
possíveis e escapar para a atmosfera através das aberturas, inclusive pelo eixo, a utilização de
juntas planas, anéis de seção “o” e outros retentores é recomendada para vedações estáticas,
contudo os selos mecânicos são em geral aplicados no eixo, pois seu uso é indicado para
os casos na qual os retentores convencionais, como gaxetas, não podem ser aplicados em
situações onde há uma condição de alta pressão, temperatura e velocidade.
Estes selos causam um impacto significante nas características dinâmicas de compressores
e turbinas. As forças desenvolvidas nas vedações são aproximadamente proporcionais ao
gradiente de pressão (∇P ) entre os selos e a densidade do fluido no interior dos mesmos. Por
causa dessa pendência da densidade, a vedação tem um grande impacto na rotodinâmica de
turbinas a vapor e compressores de alta pressão onde a densidade é maior que nas turbinas a
gás [1].
A primeira publicação sobre a influência dos selos na vibração de rotores afirma que a
6
diferença de geometria destes selos proporciona um efeito assimétrico na diferença de pressão
causada por uma força radial gerada.
Os selos apresentam diferentes disposições de montagem e modelos. Quanto ao tipo de
montagem, os mais usuais são selos simples e selos duplos [6].
Fig. 2.3: Modelo simplificado de um selo mecânico com indicações dos componentes
Fonte: Retirado de Gonçalves [6]
Como mostrado na Figura 2.3, o selo mecânico mais simples, em geral, é montado dentro
da caixa de vedação e possui um único conjunto de sede. O selo fica imerso no líquido
existente na caixa de vedação, que mantém os componentes do selo sobre compressão e
penetra na interface dos anéis primários de forma a lubrificá-los.
A montagem de selo mecânico duplo é utilizada quando o selo simples não é suficiente
para atingir a redução de vazamento desejada. Esta montagem pode ser feita com selos
opostos, em série ou frontais [6].
Quantos aos modelos, os principais são descritos na seção a seguir.
2.2 Tipos de Selo
2.2.1 Selo Plano ou Liso
O selo plano ou liso pode ser cilíndrico, cônico ou escalonado (Figura 2.4). Seu arranjo é
similar aos mancais mas possui uma folga maior para evitar o contato com a parte girante.
7
Fig. 2.4: Modelo simplificado dos selos planos cilíndricos, cônicos e escalonados, da esquerdapara a direita respectivamente
Fonte: Tiwari [9]
É o tipo mais simples de selo e sua modelagem matemática costuma ser usada como
base para o cálculo dos coeficientes rotodinâmicos mencionados anteriormente. A seguir é
demonstrado o procedimento analítico para a determinação destes coeficientes.
Rotodinâmica
Considerando a existência de excentricidade entre o rotor e o estator, folgas diferentes
surgem em cada lado do conjunto. No lado com maior folga haverá mais fluido passando,
logo, maior velocidade e menor pressão. No lado com menor folga acontece o inverso, menos
fluido passando e uma pressão maior que gera uma força que tende a fazer o anel se mover
em volta do centro, ou seja, a diferença na distribuição de pressão causada pelo deslocamento
radial produz um efeito de rigidez radial chamado “Efeito Lomakin” [6].
Fig. 2.5: Distribuição de pressão num selo devido ao efeito lomakin
Fonte: Retirado de Maurice [10]
Seguindo o exemplo da Figura 2.5 sem contabilizar o efeito da rotação do eixo, tem-se
8
que o vetor de força radial ~f é perpendicular à excentricidade (e) e sua magnitude é expressa
a seguir:
f =−∫ L
0
∫ 2π
0p(θ, z)Rsi nθdθd z (2.1)
Esta força radial ~f no rotor é de natureza restaurativa e aumenta com a redução do
espaçamento deste anel, gerando também como efeito adicional um aumento na rigidez
induzida do fluido, o que induz a redução da vibração geral da bomba e da deflexão do eixo e,
consequentemente, leva ao aumento da vida da bomba.
O amortecimento não previne diretamente a deflexão do eixo mas minimiza a resposta do
rotor à excitação das forças do mesmo modo que um amortecedor de carro atua. Ao reduzir-se
a folga do selo, aumenta-se o efeito amortecedor o que resulta numa maior estabilidade para
o rotor [11]. A redução da folga aumenta também a rigidez do eixo. Esta rigidez adicional
é derivada da força de correção que ocorre quando o rotor é excêntrico, o que leva a um
diferencial de velocidade e pressão que causa a força de correção.
É interessante determinar equações que definem como esta força se comporta em função
dos coeficientes rotodinâmicos. Para isso é feita uma abordagem física do problema da
excentricidade no qual o sistema é analisado como um sistema massa-mola-amortecedor. Para
isso se faz a seguinte aproximação linear:
M x =−K x −C x +Fext (2.2)
Onde M é a matriz massa, K é a matriz rigidez e C é a matriz de amortecimento. Os termos
x, x e x representam os vetores deslocamento, velocidade e aceleração respectivamente.
Como pode ser visto na Figura 2.6, a velocidade do movimento rotatório do rotor em
torno do próprio eixo é determinada por ω e a velocidade do movimento de seu eixo central
em torno do centro geométrico do estator é dada por Ω. O vetor r determina o comprimento
da excentricidade.
9
Fig. 2.6: Abordagem física do caso em que existe excentricidade entre o eixo do rotor e doestator
Fonte: Retirado de Yan [12] e editado pelos autores
Para este problema, o vetor x pode ser escrito em coordenadas polares, levando em
consideração sua variação no tempo:
x =
x
y
=
r cos(Ωt )
r sen(Ωt )
(2.3)
Analogamente, o vetor velocidade x pode ser escrito da forma:
x = d x
d t=
−Ωr sen(Ωt )
Ωr cos(Ωt )
(2.4)
E o vetor aceleração fica da forma:
x = d 2x
d t 2=
−Ω2r cos(Ωt )
−Ω2r sen(Ωt )
(2.5)
A partir disto, devem-se fazer as devidas simplificações físicas para facilitar a resolução
do problema. Inicialmente, considera-se que a excentricidade e as deformações são muito
10
pequenas, o que permite considerar que as matrizes M , K e C são constantes:
M = ∂F
∂x= const ante (2.6)
K = ∂F
∂x= const ante (2.7)
C = ∂F
∂x= const ante (2.8)
Para estas matrizes também é considerada a hipótese de isotropia, ou seja, as propriedades
são iguais em todas as direções. Isto permite que a matriz K por exemplo seja escrita da
seguinte forma:
K =
Kxx Kx y
Ky x Ky y
=
K k
−k K
(2.9)
Onde Kxx = Ky y = K e Kx y = Ky x = k. O coeficiente K (maiúsculo) é denominado rigidez
direta e o coeficiente k (minúsculo) é denominado rigidez cruzada.
A partir da mesma hipótese, as demais matrizes podem ser simplificadas da mesma forma:
M =
Mxx Mx y
My x My y
=
M m
−m M
(2.10)
C =
Cxx Cx y
Cy x Cy y
=
C c
−c C
(2.11)
Onde C (maiúsculo) é o coeficiente de amortecimento direto e c (minúsculo) é o coeficiente
de amortecimento cruzado. Para gases turbulentos, que é o caso estudado neste trabalho, a
matriz massa pode ser ainda mais simplificada de forma que m = 0, levando a:
M =
M 0
0 M
(2.12)
Onde coeficiente M representa a massa de fluido no volume de controle analisado.
Feitas estas simplificações, ainda existe o problema da transiência. Para simplificá-lo,
muda-se o referencial do problema, fazendo com que o movimento seja visto a partir do rotor,
ou seja, o rotor fica parado enquanto o estator gira em volta dele. Desta forma a situação
11
simplificada fica como mostrada na Figura 2.7:
Fig. 2.7: Simplificação do caso transiente com a mudança de referencial
Fonte: Retirado de Yan [12] e editado pelos autores
É visto na Figura 2.7 que foi escolhida a posição do rotor referente ao t = 0 de forma
a simplificar os vetores deslocamento, velocidade e aceleração, uma vez que cos(Ωt ) = 1 e
sen(Ωt ) = 0. Desta forma, estes vetores ficam escritos como:
x =
r
0
(2.13)
x =
0
Ωr
(2.14)
x =
−Ω2r
0
(2.15)
A partir de todas estas simplificações, a equação 2.2 pode ser escrita da seguinte forma:
M 0
0 M
−Ω2r
0
+
K k
−k K
r
0
+
C c
−c C
0
Ωr
=
Fx
Fy
(2.16)
Resolvendo a equação acima, chega-se a:
−MΩ2r
0
+
K r
−kr
+
cΩr
CΩr
=
Fx
Fy
(2.17)
A partir da qual se obtém o sistema de equações para o cálculo da força nas direções x e
12
y :
Fx
r= K + cΩ−MΩ2 (2.18)
Fy
r=−k +CΩ (2.19)
Para a determinação dos coeficientes rotodinâmicos, deve-se simular o problema para
diferentes valores de rotação (Ω). Para cada valor de rotação, irá se obter uma distribuição
de pressão que pode ser integrada pela Equação 2.1 para se obter uma força. Os valores de
rotação e força são então plotados em um gráfico F x Ω. Para obtenção dos coeficientes da
equação de Fx , deve-se ligar os pontos seguindo uma equação de segundo grau enquanto que
para obtenção dos coeficientes da equação de Fy os pontos podem ser aproximados segundo
uma reta.
A partir das equações obtidas para as curvas, determinam-se os coeficientes rotodinâmicos
do conjunto rotor-selo.
2.2.2 Selo de Anéis Flutuantes
O selo de anéis flutuantes é geometricamente semelhante ao selo plano, porém neste caso
o espaço livre entre o selo e o eixo do rotor é preenchido com óleo com o objetivo de diminuir
os vazamentos. O anel de vedação orbita e vibra com o rotor, atenuando desbalanceamentos
[6].
2.2.3 Selo de Contato
No selo de contato não há folga projetada, portanto neste tipo de selo surgem altas
temperaturas devido ao atrito, o que reduz a vida útil do equipamento. Seu uso é recomendado
em máquinas de fluxo com baixa velocidade de rotação ou máquinas em que o fluido de
trabalho possa operar como fluido refrigerante [6].
13
2.2.4 Selo Labirinto
Selos labirinto, também chamados de anéis labirinto, são um tipo de selo mecânico feito
geralmente por uma união permanente entre um anel interno e um anel externo. Além da
união destes anéis, o principal componente do selo são as lâminas que formam o labirinto e
devido a esta geometria, este selo possui perfis muito parecidos entre si o que ajuda a manter
as folgas entre os dois anéis quase idênticas. A Figura 2.9 mostra alguns possíveis arranjos
de lâminas para o selo labirinto, pode-se ver que em alguns casos o próprio rotor também
apresenta lâminas que se encaixam com as do selo, dificultando ainda mais a passagem de
fluido.
Fig. 2.8: Possíveis disposições das lâminas do selo labirinto
Fonte: Retirado de Gonçalves [6]
Os lados das faces opostas deste tipo de selo não se encostam, o que faz reduzir bastante o
atrito já que entre o rotor e o estator há uma pequena folga. Em muitos casos, esta folga chega
até a 0,0762mm [13].
Para cada lâmina pela qual o fluido passa ocorre uma queda de pressão, o que dificulta a
passagem de fluido para as partes protegidas pela vedação. Quanto mais complexo o labirinto
de lâminas, maior será a vedação garantida pelo selo [6].
Este selo é excelente para reter óleo de lubrificação e bloquear contaminação por partículas
e misturas de fluidos. Porém uma de suas desvantagens é que os selos labirintos em geral são
mais caros e requerem adequação da geometria do equipamento (retrofit) para a sua instalação
14
apropriada.
2.2.5 Selo Honeycomb
Derivados da tecnologia aeronáutica, os selos honeycomb apresentam furos com o padrão
de cavidade de prisma hexagonal (formato colméia ou honeycomb) ao longo de sua superfície.
Fig. 2.9: Configuração do selo colméia
Fonte: Retirado de Massini [14]
De forma semelhante ao selo labirinto, estes furos geram regiões onde há queda de pressão
no fluido, dificultando sua passagem. Este padrão começou a ser amplamente utilizado a partir
dos anos 60, onde substituiu os antigos selos labirinto que eram consumidos pelo desgaste
e possuíam uma capacidade de vedação pior que os selos padrão colméia para o mesmo
tamanho de folga e também apresentavam uma melhora nas características rotodinâmicas [15].
No final da década de 1980, Childs [16] realizou a primeira pesquisa envolvendo selos colméia
analisando como a diminuição vazamento de vazão mássica por ser obtida melhorando a
estabilidade rotodinâmica. Após uma extensa investigação, na qual 7 tipos de selos colméias
foram analisados, os resultados mostraram que os selos colméias eram melhores em termos
de estabilidade rotodinâmica do que os selos labirintos quando não usado nenhum adicional
para a redução de redomoinho para o caso do selo labirinto [14].
Assim, diversos estudos foram realizados e evidenciaram que o selo colméia é referência
em situações que exigem estabilidade rotodinâmica com uma boa eficiência de vedação em
condições de alta pressão, temperatura e velocidade [14].
15
2.2.6 Selo Hole-pattern
Os selos padrão hole-pattern, que são o foco desse trabalho, começaram a ser utilizados
na mesma época que os selos padrão colméia e apresentam características rotodinâmicas
semelhantes. O conjunto de regiões cilíndricas é usado para diminuir o vazamento interno,
tanto o selo hole-pattern quanto o selo colméia combinam bom desempenho de vedação e
baixo custo.
−
Fx
Fy
=
K k
−k K
X
Y
+
C c
−c C
X
Y
+M
X
Y
(2.20)
Na qual K,C são os termos diretos e k e c são os termos cruzados de rigidez e amorte-
cimento respectivamente. O termo M é o relacionado ao termo direto de adição de massa.
A maioria dos trabalhos que desenvolvem modelos rotodinâmicos para selos mecânicos
geralmente desprezamo termo M e o termo c por não representar nenhum impacto direto no
amortecimento geral.
Os selos de parede rugosas, como o selo tipo hole-pattern ou labirinto acoplados no estator,
tendem a reduzir o vazamento por meio da redução da velocidade circunferencial do selo,
também chamado relacionado ao swirl, e aumentam o coeficiente de amortecimento cruzado
k. Assim, a substituição destes selos por um selo amortecedor faz com que esta velocidade
aumente e, consequentemente,a sua instabilidade rotodinâmica.
O efeito dos selos no coeficiente de amortecimento cruzado k, pode ser mais bem entendi-
mento se reescrevemos as forças radiais e tangenciais atuantes no rotor em um movimento
síncrono a uma velocidade angular constante e amplitude e:
−
Fn
Ft
=
K +cω−Mω2
Cω−k
e (2.21)
O trabalho de Childs e Kim [17] mostra que o coeficiente de amortecimento efetivo
Ce f f = C−kω é maior em selos amortecedores por causa da redução do coeficiente k.
Um selo mecânico típico é caraterizado por pequenas folgas entre o rotor e o estator com-
parado com o diâmetro destes. Geralmente, a parte do selo mecânico acoplada a superfície do
16
estator é caracterizada por diversos tipos de cavidades, na qual produzem expansão/contração
do escoamento nas direções tanto axiais quanto circunferenciais e dessa forma reduzem o
vazamento e os redemoinhos. A superfície destes selos chega a ter milhares de cavidades em
três dimensões, cujos comprimentos são da ordem da folga do selo. Portanto, estas geometrias
discretas afetam toda física envolvendo o escoamento nesta folga, tornando o estudo analítico
destes aspectos inviáveis [18].
Embora o uso de selos hole-pattern seja comum em compressores centrífugos, não existem
ainda ferramentas de ampla utilização que permitam calcular os coeficientes rotodinâmicos
de um modo mais eficiente [19]. Atualmente, os dois métodos mais difundidos para o cálculo
desses coeficientes são o método bulk-flow e o método de dinâmica dos fluidos computacio-
nal(CFD) [19].
Tradicionalmente, o estudo analítico do escoamento nos selos amortecedores se baseiam
no modelo bulk flow, cuja teoria se baseia na adaptação do Blasius type friction para o fator
de atrito, f como mostrado na Equação 2.22:
f = τw12ρU 2
= n (Re)m = n
(ρU (2H)
µ
)m
(2.22)
Onde U é a velocidade bulk-flow relativa à parede do estator ou rotor, e ρ é a massa
específica do fluido em questão e τw é a tensão cisalhante na parede. O número de Reynolds
é baseado no diâmetro hidráulico do problema, que neste caso é duas vezes a folga do selo.
Os termos n e m são coeficientes empíricos obtidos experimentalmente [18].
A primeira análise completa baseada no modelo bulk-flow foi realizado por Nelson [20],
na qual foi utilizado o modelo para um selo anelar afunilado operando em regime turbulento.
Sua análise foi baseada nas equações de ordem zero e primeira ordem de pertubação do
modelo bulk-flow.
Os resultados obtidos por Neslon [20] tiveram boa precisão usando a teoria de non-pre
rotation (sem pré-rotação), porém paro o caso não estático inicial (pre rotated case) a precisão
não é satisfatória.
Assim, trabalhos posteriores modificaram o modelo de forma a atingir boas precisões em
diferentes condições físicas e geometrias.
17
Como estes selos mecânicos geralmente funcionam em escoamento de alto número de
Reynolds onde as forças de inércia são dominantes sobre as forças viscosas. O modelo bulk-
flow aliado ao uso de ferramentas computacionais, garante uma boa relação entre a precisão e
o custo computacional [21].
No trabalho de Migliorini [8] este modelo bulk-flow é utilizado para a determinação
dos coeficientes rotodinâmicos de um selo hole-pattern. A análise é feita considerando a
hipótese de concentricidade, o que impede o cálculo destes coeficientes pelo método analítico
descrito na seção do Selo Simples, mas permite que apenas uma seção do selo seja analisada
utilizando-se da condição de simetria circunferencial. A análise de apenas uma seção do
domínio leva a um custo computacional muito menor. Combinando o modelo analítico e o
computacional podem-se determinar os coeficientes rotodinâmicos, através da modelagem
com este modelo.
2.3 Hipótese do Contínuo
A hipótese do contínuo consiste em uma simplificação que considera a natureza discreta
das moléculas dos fluidos com uma visão do contínuo. Assim, para que qualquer propriedade
local do fluido permaneça inalterada independente do tamanho da amostra, o fluido deve ser
tratado como um meio contínuo, que é válido a partir de um determinado tamanho garantindo
assim uma teoria de movimento para os fluidos (mecânica dos fluidos). A figura 2.10 mostra
o domínio de validade da hipótese do contínuo a partir da variação da massa específica ρ e da
variação de volume.
18
Fig. 2.10: Validade da hipótese do contínuo
Fonte: Retirado de Azevedo [22]
Pode-se definir um número adimensional dado pela razão entre o comprimento livre
médio entre as partículas e uma escala de comprimento do escoamento, chamado número de
Knudsen [23].
Este número adimensional, é também utilizado para determinar se a formulação da
mecânica do contínuo ou a formulação da mecânica estatística deve ser usada.
K u = λ
l(2.23)
Se o número de Knudsen é próximo ou maior que um, o caminho médio livre de uma mo-
lécula é comparável à escala de comprimento do problema, e a consideração de continuidade
da mecânica dos fluidos não é mais uma boa aproximação. Neste caso a mecânica estatística
deve ser usada.
2.4 Números Admensionais Relevantes
2.4.1 Número de Reynolds
O número adimensional de Reynolds é dado pela razão entre as forças de inércia e as
forças viscosas atuando sobre um fluido. Sua principal função é determinar se um dado
escoamento é laminar ou turbulento, de acordo com seu valor e com o tipo de fluido a ser
19
analisado. A água, por exemplo, possui comportamento turbulento para valores de Reynolds
maiores que 2300.
A principio, as equações de Navier Stokes descrevem ambos os regimes de escoamento,
laminar e turbulento, sem a necessidade de informação adicional. Porém, os escoamentos
turbulentos com número de Reynolds elevado abrangem uma ampla faixa de comprimento,
tempo e velocidade.
O número de Reynolds é definido como:
Re = U L
ν= f or ças de i nér ci a
f or ças vi scosas= convecção de momentum
di f usão de momentum(2.24)
Onde U é a velocidade média do fluido, L é o comprimento característico e ν é a viscosi-
dade cinemática.
2.4.2 Número de Prandtl
O número adimensional de Prandtl é dado pela razão entre a difusividade de momento e a
difusividade térmica de um fluido. Ele é responsável por determinar a espessura relativa das
camadas limites térmica e de momento em problemas de transferência de calor. É definido
como:
Pr = CPµ
k= t axa de di f usão vi scosa
t axa de di f usão tér mi ca(2.25)
Onde CP é o calor específico, µ é a viscosidade dinâmica e k é a condutividade térmica.
Como a equação não depende de um comprimento característico, o número de Prandtl
está relacionado somente ao tipo de fluido e seu estado. Em função disso este número é
frequentemente encontrado em tabelas de propriedades de fluidos.
A taxa de difusão viscosa está ligada a difusão de calor por convecção enquanto a difusão
térmica representa a difusão de calor por condução. Para metais líquidos como o mercúrio,
que é conhecido por apresentar alta condutividade térmica, o número de Prandtl apresenta um
valor baixo enquanto que para alguns óleos hidrocarbonetos (óleo de motor), a eficiência na
transferência de energia por convecção é maior o que gera um número de Prandtl alto.
20
Para ar e gases em geral o número de Prantl apresenta valores próximos de 1, indicando
que não há predominância de convecção sobre condução ou vice-versa.
2.4.3 Número de Schmidt
O número adimensional de Schmidt é considerado o análogo do número de Prandtl
para transferência de massa. É definido como a razão entre a difusividade de momento e a
difusividade de massa e segue a seguinte equação:
Sc = µ
ρD= t axa de di f usão vi scosa
t axa de di f usão de massa(2.26)
Onde µ é a viscosidade dinâmica do fluido, ρ é a massa específica e D é a difusividade de
massa.
Para gases, assim como o número de Prandtl, o número de Schmidt é da ordem de uma
unidade.
2.4.4 Número de Péclet
O número de Péclet é um número adimensional relevante no estudo dos fenômenos de
transporte no contínuo. É definido como a razão entre a taxa de advecção de uma grandeza
física e a taxa de difusão desta mesma grandeza seguindo um gradiente apropriado. Este
número segue a seguinte definição:
Pe = t axa de tr anspor te por ad vecção
t axa de tr anspor te por di f usão(2.27)
O número de Péclet apresenta variantes para transferência de calor e para transferência de
massa, tendo relação com os números de Prandtl e Schmidt respectivamente.
Para transferência de calor, é definido como:
Pe = RePr (2.28)
21
Para transferência de massa, é definido como:
Pe = ReSc (2.29)
Para valores altos do número de Péclet, recorrentes em aplicações de engenharia, as
variáveis no escoamento tendem a se tornar unidirecionais, pois o transporte por convecção
é predominante e ocorre em uma direção específica, sendo essa a direção principal em que
se variam as propriedades. Isso permite que modelos computacionais mais simples sejam
adotados [24].
2.5 Pequenas Escalas de Kolmogorov
As escalas de Kolmogorov são as escalas de dimensão, tempo e velocidade dos menores
vórtices encontrados em um escoamento turbulento, a partir dos quais a energia cinética
turbulenta é dissipada.
Estas escalas são utilizadas pelo método de simulação direta numérica para captar todos
os efeitos da turbulência. Pode se determinar a escala de dimensão pela teoria de Kolmogorov
através de três grandezas básicas:
1. Viscosidade Cinemática
ν=(µ
ρ
)(2.30)
Onde ν é viscosidade cinemática, µ é a viscosidade dinâmica e ρ é a massa específica
do fluido.
2. Velocidade do Escoamento Livre
3. Dimensão Característica do objeto
De acordo com Kolmogorov, as menores escalas de comprimento, tempo e velocidade
22
podem ser calculadas a partir de [25]:
lk =(ν3
ε
) 14
(2.31)
uk = (νε)14 (2.32)
tk =(νε
) 12 (2.33)
As escalas acima 2.31, 2.32 e 2.33 podem ser obtidas através da análise dimensional e
elas só dependem da taxa de dissipação viscosa (ε) e a viscosidade cinemática (ν), parâmetros
importantes para o escoamento em pequenas escalas.
Nestas escalas de Kolmogorov, o número de Reynolds é escrito na seguinte forma:
Rek =(
uk lk
ν
)(2.34)
Substituindo as Equações 2.31 e 2.33 na Equação 2.34, tem-se:
Rek = (1) (2.35)
Observa-se que o Reynolds é unitário nas pequenas escalas de Kolmogorov, mostrando
que os efeitos viscosos passam a dominar os efeitos de inércia. Assim, turbilhões menores de
lk são dissipados pela viscosidade e não se desenvolvem [25].
2.6 Princípios de Mecânica dos Fluidos
Os princípios físicos mais básicos que envolvem os problemas de mecânica de fluidos se
baseiam nas seguintes descrições matemáticas das leis de conservação da física:
• Considerando o espaço a ser analisado, a massa do fluido é conservada.
• A taxa de variação da quantidade de movimento linear deve ser igual à soma das forças
atuantes sobre a partícula de fluido.
• A Taxa de variação de energia é igual à soma da taxa de adição de calor e a taxa de
23
trabalho realizado pela partícula de fluido.
O comportamento dos fluidos são analisados através das variáveis de campo de veloci-
dade e campos de pressão que juntos definem o comportamento do fluido. Assim, para um
escoamento em três dimensões, tem-se quatro incógnitas u= u(u, v, w) e P pertencentes aos
campos de velocidade e pressão respectivamente.
As equações que juntas são utilizadas para resolução de tais campos são o balanço de
massa e as equações de momento nas três direções, totalizando quatro equações para 4
incógnitas, constituindo-se assim de um sistema com solução.
Alguma vezes a pressão pode variar com a temperatura, nesse caso utilizam-se as equações
de energia para a definição do campo de temperatura e consequente campo de pressão.
2.6.1 Lei de Conservação da Massa
A conservação de massa é um dos princípios universais da mecânica dos fluidos. É
apresentada a seguir a dedução deste princípio, supondo válida a hipótese do contínuo e sua
discretização em volumes de controle infinitesimais de um fluido.
Para o cálculo da lei de conservação de massa para três dimensões deve-se primeiramente
analisar o balanço de massa de um elemento do fluido como mostrado na Figura 2.11 abaixo:
Fig. 2.11: Balanço de massa do elemento de um fluido
Fonte: Retirado de Malalasekera [26]
24
A taxa que representa o aumento de massa para o elemento do fluido é dada por:
∂
∂t
(ρδxδyδz
)= ∂ρ
∂tδxδyδz (2.36)
A partir da figura 2.11, tem-se que o balanço de massa das faces do elemento de fluido é
dada por:
(ρu − ∂(ρu)
∂x
1
2δx
)δyδz −
(ρu + ∂(ρu)
∂x
1
2δx
)δyδz
+(ρv − ∂(ρv)
∂y
1
2δy
)δxδz −
(ρv + ∂(ρv)
∂y
1
2δy
)δxδz
+(ρw − ∂(ρw)
∂z
1
2δz
)δxδy −
(ρw + ∂(ρw)
∂z
1
2δz
)δxδy
=−(∂(ρu)
∂x+ ∂(ρv)
∂y+ ∂(ρw)
∂z
)δxδyδz
(2.37)
Somando-se a Equação 2.37 acima com 2.36 chega-se a:
−(∂(ρu)
∂x+ ∂(ρv)
∂y+ ∂(ρw)
∂z− ∂ρ
∂t
)δxδyδz (2.38)
Dividindo-se os termos acima por δxδyδz, tem-se que a equação de continuidade fica da
seguinte forma:
∂ρ
∂t+ ∂(ρu)
∂x+ ∂(ρv)
∂y+ ∂(ρw)
∂z= 0 (2.39)
Ou em notação vetorial:
∂ρ
∂t+∇· (ρu
)= 0 (2.40)
A equação da continuidade para regimes permanentes (não transientes) e fluidos incom-
pressíveis fica na forma:∂u
∂x+ ∂v
∂y+ ∂w
∂z= 0 (2.41)
25
2.6.2 Conservação da Quantidade de Movimento
A quantidade de movimento linear ou momentum é uma propriedade conservativa que
pode ser definida a partir da Segunda Lei de Newton. Se as vizinhanças exercem uma força
resultante sobre o sistema, esta lei estabelece que a massa no sistema começa a se acelerar:
F = ma= mdudx
= d
dt(mu) (2.42)
Em mecânica dos fluidos, esta lei é chamada de equação da conservação da quantidade
de movimento linear e pode ser obtida a partir da mesma dedução usando o mesmo volume
de controle elementar da Figura 2.11, podendo então se escrever a seguinte relação para a
quantidade de movimento linear:
∑F= ∂
∂t
(∫V.C
uρ dV ol
)+∑
(mi ui )saíd a −∑
(mi ui )entr ad a (2.43)
Como o volume elementar é muito pequeno, então a integral de volume pode ser aproxi-
mada para um termo diferencial:
∂
∂t
(uρ dV ol
)≈
∂
∂t
(ρu dx dy dz
)(2.44)
Avaliando os fluxos de quantidade de movimento que ocorrem nas seis faces do volume
elementar, obtém-se a Tabela 2.1:
Tab. 2.1: Balanço de fluxos no volume de controle elementar 2.11
Face Fluxo de entrada de Momentum Fluxo de saída de Momentum
x ρuu dy dz
(ρuu+ (
∂ρuu∂x
) dx
)dy dz
y ρvu dx dz
(ρvu+ (
∂ρvu∂y
) dy
)dx dz
z ρwu dx dy
(ρwu+ (
∂ρwu∂z
) dy
)dx dy
26
Manipulando os termos chega-se a:
= ∂
∂t
(ρu
)+ ∂
∂x
(ρuu
)+ ∂
∂y
(ρvu
)+ ∂
∂z
(ρwu
)= u
[∂ρ
∂t+∇· (ρu
)]︸ ︷︷ ︸
Conti nui d ade
+ρ(∂u∂t
+u∂u∂x
+ v∂u∂y
+w∂u∂z
)(2.45)
A equação da continuidade aparece naturalmente na dedução de momentum, e como visto
anteriormente, o termo da equação da conservação de massa é nulo. Assim, a equação de
conservação da quantidade de movimento pode ser reduzida a:
∑F= ρDu
Dtdx dy dz (2.46)
As forças resultantes podem ser de dois tipos: de campo e de superfície. As forças de
campo são decorrentes de campos externos que agem sobre toda massa dentro do volume,
como a gravidade.
dFg r avi d ade = ρg dx dy dz (2.47)
Já as forças de superfícies ocorrem das tensões sobre os lados da superfície analisada.
Fig. 2.12: Componentes das tensões na superfície de controle
Fonte:Retirado de White [27]
27
Desta forma, o tensor de tensões pode ser definido como:
σi j =
−p +τxx τy x τzx
τx y −p +τy y τz y
τxz τy z −p +τzz
Fazendo então o mesmo balanço de entrada e saída realizado na obtenção da equação da
conservação da massa, tem-se que por exemplo a força σxx dy dz, que está orientada para
esquerda e atua sobre a face esquerda, está equilibrada pela força σxx dy dz para direita sobre
a face direita, deixando apenas o componente diferencial resultante da série de taylor truncada(∂σxx
∂x
)sobre a face direita.
Seguindo essa mesma linha de raciocínio para as outras quatros faces, tem-se que a força
de superfície líquida na direção x é dada por:
dFx,sup =[∂
∂x(σxx)+ ∂
∂y(σy x)+ ∂
∂z(σzx)
]dx dy dz (2.48)
Fig. 2.13: Componentes das tensões na direção x
Fonte: Retirado de Malalasekera [26]
28
A Equação 2.48 pode ser escrita utilizando termos do tensor de tensões como:
dFx
dV ol=−∂p
∂x+ ∂
∂x(τxx)+ ∂
∂y(τy x)+ ∂
∂z(τzx)
dFy
dV ol=−∂p
∂y+ ∂
∂x(τx y )+ ∂
∂y(τy y )+ ∂
∂z(τz y )
dFz
dV ol=−∂p
∂z+ ∂
∂x(τxz)+ ∂
∂y(τy z)+ ∂
∂z(τzz)
(2.49)
Para a forma mais compacta tem-se que:
(dF
dV ol
)sup
=−∇p +(
dFdV ol
)vi scosa
(2.50)
Onde as forças viscosas são dadas por:
(dF
dV ol
)vi scosa
=∇·(τ′
i j
)= di v.
τxx τy x τzx
τx y τy y τz y
τxz τy z −τzz
(2.51)
Na qual τ′i j representa o tensor de tensões viscosas agindo sobre o volume de controle
elementar. Substituindo as Equações 2.50 e 2.51 na Equação 2.46, obtem-se:
∂ρu
∂t+∇· (ρuu
)=−∇p +∇· (τ′)+ρg (2.52)
Onde ρ é a massa específica do fluido, u é o campo de velocidade do escoamento, p é a
pressão característica do fluido e τ′ representa as tensões viscosas. Escrevendo essa equação
da quantidade de movimento em sua forma completa, chega-se a:
−ρgx − ∂p
∂x+ ∂
∂x(τxx)+ ∂
∂y(τy x)+ ∂
∂z(τzx) = ρ
(∂u
∂t+u
∂u
∂x+ v
∂u
∂y+w
∂u
∂z
)(2.53)
−ρg y − ∂p
∂y+ ∂
∂x(τx y )+ ∂
∂y(τy y )+ ∂
∂z(τz y ) = ρ
(∂v
∂t+u
∂v
∂x+ v
∂v
∂y+w
∂v
∂z
)(2.54)
−ρgz − ∂p
∂z+ ∂
∂x(τxz)+ ∂
∂y(τy z)+ ∂
∂z(τzz) = ρ
(∂w
∂t+u
∂w
∂x+ v
∂w
∂y+w
∂w
∂z
)(2.55)
Onde os termos convectivos e não lineares das equações se encontram no lado direito.
29
Como a maioria dos escoamentos reais acontecem com fluidos newtonianos, as tensões para
um escoamento incompressível tridimensional são expressas por:
τxx = 2µdu
dx(2.56)
τy y = 2µdv
dy(2.57)
τzz = 2µdw
dz(2.58)
τy x +τx y =µ(
du
dy+ dv
dx
)(2.59)
τzx +τxz =µ(
du
dy+ dw
dx
)(2.60)
τy z +τz y =µ(
dv
dz+ dw
dy
)(2.61)
Substituindo na equação 2.55 obtem-se as equações de Navier-stokes:
−ρgx − ∂p
∂x+µ
(∂2u
∂x2+ ∂2u
∂y2+ ∂2u
∂z2
)= ρ
(∂u
∂t+u
∂u
∂x+ v
∂u
∂y+w
∂u
∂z
)(2.62)
−ρg y − ∂p
∂y+µ
(∂2v
∂x2+ ∂2v
∂y2+ ∂2v
∂z2
)= ρ
(∂v
∂t+u
∂v
∂x+ v
∂v
∂y+w
∂v
∂z
)(2.63)
−ρgz − ∂p
∂z+µ
(∂2w
∂x2+ ∂2w
∂y2+ ∂2w
∂z2
)= ρ
(∂w
∂t+u
∂w
∂x+ v
∂w
∂y+w
∂w
∂z
)(2.64)
2.6.3 Conservação da Energia
As equações de governo para a energia são obtidas de forma semelhante as equações de
momento linear. Deve-se também lembrar que a equação de conservação de momento angular
não é mencionada, pois possui uma simetria entre os valores de tensões. Assim, as equações
de energia são apenas apresentadas, e podem ser expressas tanto para energia interna quanto
para temperatura.
A equação da energia pode ser escrita para Energia interna (e) da seguinte forma indicial:
De
Dt= pE
ρ
∂ui
∂xi+ λ
ρ
(∂ui
∂xi
)2
+ 2µ
ρSi j : Si j + 1
ρ
∂
∂xi
[k(θ)
∂θ
∂xi
]+ qr ad (2.65)
30
Enquanto que para temperatura (θ) pode ser escrita da seguinte forma:
ρcvDθ
Dt=−
∣∣∣∣θ∂pE
∂θ
∣∣∣∣p
∂ui
∂xi+λ
(∂ui
∂xi
)2
+2µSi j : Si j + ∂
∂xi
[k(θ)
∂θ
∂xi
]+ρqr ad (2.66)
2.7 Turbulência
O fenômeno de turbulência é encontrado na maioria dos escoamentos e por este motivo
vem sendo estudado por vários séculos [23]. Teorias e conceitos têm sidos formulados na
tentativa de se obter uma descrição mais geral para este fenômeno que seja adequada para
qualquer problema de interesse. Enquanto esta descrição geral não é formulada, modelos
simplificados têm sido propostos como forma de análise para problemas em cada área de
interesse [28, 29].
A complexidade dos escoamentos turbulentos na análise dos fenômenos de transporte de
massa, momentum e energia não permite uma abordagem estritamente analítica do problema,
o que seria extramente vantajoso do ponto de vista físico e matemático [30].
O modelo matemático para análise de escoamentos (equações de conservação de massa,
energia e quantidade de movimento da mecânica dos fluidos) faz parte de um sistema complexo
de equações diferenciais parciais. Este sistema inclui equações constitutivas, condições
de contorno e condições iniciais a serem resolvidas em domínios de geometria variável,
aumentando a complexidade do problema. Além disso,vale ressaltar, que os escoamentos
turbulentos têm seu comportamento coerente nas grandes escalas que por sua vez influenciam
na forma global do escoamento. Já em menor escala, há um comportamento caótico que
contribui para dissipação de energia ao longo do tempo e espaço, e a transição de energia de
grande escala (que possui a priori mais energia que a de menor escala) para pequena escala
ocorre de forma difusiva, tridimensional e transiente [31].
As escalas envolvidas nestes escoamentos são normalmente muito maiores do que as
escalas de comprimento do movimento molecular [32], desta forma a turbulência pode ser
descrita como um fenômeno contínuo [31].
Além destas características, a turbulência também possui importantes propriedades em
problemas físicos, dos quais são citados abaixo:
31
• Facilidade para misturar ou transportar propriedades a taxas muito mais elevadas do
que aquelas alcançadas pela difusão molecular [31].
• São sempre dissipativos e necessitam de uma alimentação contínua de energia para
suprir as perdas viscosas de modo que, se nenhuma energia for fornecida ao escoamento,
a turbulência decai rapidamente [28].
• A turbulência é um fenômeno do contínuo, isto é, até mesmo os menores vórtices são
tipicamente muito maiores que a escala de comprimento molecular, fazendo com que
a Hipótese do Contínuo seja válida e o escoamento turbulento seja governado pelas
equações de Navier-Stokes [28].
De acordo com Eiger [29], para se avaliar turbulência deve-se analisar a vorticidade
pois esta desempenha um papel fundamental considerando que escoamentos turbulentos
são sempre rotacionais. Há uma forte dependência dos grandes vórtices da geometria do
escoamento e da maneira de como são gerados, assim estes vórtices possuem características
anisotrópicas e um comportamento um tanto determinístico [25]. A dependência da geometria
vai desaparecendo ao considerar vórtices cada vez menores e o padrão do campo turbulento
torna-se cada vez mais caótico e consequentemente mais isotrópico [25].
Com objetivo de estudar cada vez mais profundamente as escalas que envolvem o fenô-
meno da turbulência, e pela falta de um método analítico que forneça tais soluções, o uso de
novas metodologias numéricas vem sendo adotado ao longo dos últimos anos para a obtenção
de soluções aproximadas das equações de governo dos problemas físicos, fazendo-se uso
da fluidodinâmica computacional, que tem se tornado realidade devido ao avanço na área
computacional [33, 34].
32
Fig. 2.14: Fluxograma de modelagem
Fonte: Produzido pelos autores
A Figura 2.14 representa o fluxograma na qual se evidencia a abordagem da modelagem
na fluidodinâmica computacional, onde começa-se sempre com um problema real numa
dada condição física medida ou avaliada, que é então descrita por um modelo matemático
usando simplificações adequadas para sua descrição. Como muito destes problemas são
descritos por equações diferenciais parciais cujas soluções não são dadas analiticamente,
é então discretizado este modelo matemático de forma que as ferramentas computacionais
disponíveis possam interpretar e resolver estas equações numericamente, obtendo assim uma
solução numérica aproximada que é validada tanto numericamente quanto fisicamente.
2.7.1 Tratamento da Turbulência
Como os métodos de fluidodinâmica computacional apresentam sempre uma solução
aproximada, deve-se então adotar a melhor abordagem para o problema físico de interesse.
Como o presente trabalho se baseia na análise dos modelos de turbulência, são então apresen-
tadas as principais abordagens para a modelagem de turbulência. As abordagens apresentadas
pela Figura 2.15 destacam as escalas de energia simuladas pelos modelos em relação ao
33
comprimento de onda.
Fig. 2.15: Abordagens das simulações de escoamentos turbulentos disponíveis
Fonte: Retirado de Poinsot [35]
As equações de conservação de massa, energia e quantidade de movimento constituem
um modelo matemático completo, capaz de descrever os escoamentos turbulentos de fluidos.
Uma das abordagens que existem, que também seria a ideal, é a abordagem de resolução
direta dessas equações, também chamada de Simulação Numérica Direta (Direct Numerical
Simulation, DNS), que necessita de uma discretização espacial e temporal para a simulação
direta de todas as escalas de turbulência envolvidas que permitem a resolução do escoamento
turbulento em grande detalhe e com uma excelente precisão, porém a malha em cada direção
necessária para tal análise deverá ter um grau de refino muito alto, a fim de capturar toda
gama de escalas turbulentas para garantir a resolução de todos os detalhes da turbulência,
gerando um custo computacional muito elevado com tempos de processamento inaceitáveis,
mesmo em computadores avançados [36, 37].
Assim, na abordagem DNS a malha deve ser suficientemente refinada afim de captar os
efeitos das menores escalas de movimento, resultando num campo completo do escoamento
turbulento, tridimensional e transiente. Entretanto, devido as grandes exigências de resolução
espacial e temporal, a prática desta abordagem limita-se aos escoamentos com moderado à
34
baixo número de Reynolds (ReL) pois de acordo com Silveira [38], o mínimo de pontos de
discretização necessários para uma resolução especial completa do escoamento é proporcional
a ReL9/4.
Além disso, uma vez que a turbulência é marcada por diversas escalas temporais e espaci-
ais, estas aumentam consideravelmente com o número de Reynolds, tornando-se proibitivo
do ponto de vista prático para a grande maioria dos problemas de interesse em engenharia.
Como consequência da inviabilidade da solução via simulação direta, é necessário o emprego
de metodologias alternativas.
Estas metodologias alternativas se baseiam na decomposição das equações de governo
em campos médios e suas flutuações, a partir das quais surgem momentos de segunda ordem
ou mais, os quais envolvem flutuações, gerando ainda mais incógnitas que equações. Este
problema é conhecido como problema do fechamento, e ambas metodologias alternativas
possuem este problema.
A Simulação de Grandes Escalas (Large Eddy Simulation, LES) é uma opção alternativa
para a análise de escoamentos turbulentos cujo nível de detalhamento não é tão elevado mas
captura os efeitos mais significativos do escoamento turbulento. Esta simulação consiste em
resolver as equações de conservação (equações diferenciais parciais) em grandes escalas dos
escoamentos e empregar modelos para representar os efeitos das escalas inferiores à resolução
obtida pela discretização espacial do problema (escalas sub-malha) [39–41]. A técnica LES
utiliza malhas mais grosseiras que as da simulação numérica direta, porém finas o suficiente
para prever características de escoamento instantâneas até em escala inercial e resolver a
estrutura do fluxo turbulento.
Ao utilizar-se destas técnicas, podem se encontrar algumas dificuldades com as condições
de contorno, pois trata-se de uma análise que funciona melhor na região onde o efeito da
turbulência de grande escala é mais presente [25].
Os Modelos de submalha tornam-se importantes quando se utiliza LES para números de
Reynolds elevados quando malhas relativamente grossas são usadas em regiões de grande
cisalhamento [25, 42].
A maioria dos casos de interesse industrial se referem a valores médios (como de taxas de
35
transferência de calor, arrasto, pressão, esforços em geral) que devem ser obtidos usando o
mínimo de custo computacional possível. Assim, para obter as informações visadas, com uma
razoável precisão, no menor prazo e custo possível, uma abordagem utilizando valores médios
é então utilizada, chamada abordagem média de Reynolds (Reynolds Averaged Navier-Stokes,
RANS) que atualmente é a mais comum nas aplicações industriais.
As equações da técnica RANS são obtidas através de um conjunto de médias das equações
de Navier-Stokes e da continuidade. O elemento crítico da modelagem RANS é a representa-
ção das tensões de Reynolds ou tensões turbulentas que descrevem os efeitos das flutuações
turbulentas de pressão e velocidade.
Tanto o escoamento laminar quanto o escoamento turbulento satisfazem a Equação 2.41.
Porém para o caso de escoamentos turbulentos, por causa do comportamento caótico, há
flutuações para cada termo de velocidade e pressão, gerando uma função aleatória de variação
rápida do tempo e do espaço para esses termos. Como a resolução de variáveis flutuantes
instantâneas requer uma matemática muito complexa, como mostrado na Figura 2.16, a
abordagem utilizada é voltado para o uso de valores médios para escoamentos turbulentos
com alto número de Reynolds. Essa abordagem foi adotada pela primeira vez por Osborne
Reynolds em 1985 [27].
Fig. 2.16: Medição de velocidade de um escoamento turbulento
Fonte: Retirado de Malalasekera [26]
Observa-se na Figura 2.16 que para se trabalhar com a velocidade de um escoamento
36
turbulento deve-se trabalhar com uma média que varia no tempo. Foi assim que Reynolds
iniciou sua análise de escoamentos turbulentos propondo uma média temporal u definida por
u(x, y, z, t
)mostrada a seguir:
u= 1
T
∫ T
0ud t (2.67)
Onde T é um período usado para o calculo da média.
Assim, Reynolds fundamentou a base para o método de equações de médias de Reynolds
(RANS), onde é realizada a decomposição que descreve os valores instantâneos das variáveis
do movimento turbulento como uma variação randômica em torno de valores médios [19].
Dessa forma, a variável flutuante de velocidade u′ é então definida como desvio de u em
relação à sua média temporal:
ui′ = ui −ui (2.68)
A partir da média do quadrado da flutuação define-se a intensidade da turbulência [27].
ui′2 = 1
T
∫ T
0ui
′2d t (2.69)
Para então chegar na definição dos novos parâmetros para a velocidade turbulenta, Rey-
nolds propôs que cada velocidade é composta em uma média mais uma variável de flutuação
[27].
u = u +u′ (2.70)
v = v + v ′ (2.71)
w = w +w ′ (2.72)
Substituindo as equações acima na equação da continuidade 2.41, gera-se uma nova
equação da continuidade que contabiliza as flutuações no tempo e espaço do escoamento
turbulento [27].
∂u
∂x+ ∂v
∂y+ ∂w
∂z= 0 (2.73)
37
O mesmo método é válido para pressão:
ui = ui +ui′ (2.74)
pi = pi +pi′ (2.75)
E sua generalização:
φ=φ+φ′ (2.76)
Assim define-se o operador média como:
φ= 1
∆T
∫∆Tφd t (2.77)
A equivalência para quantidade de movimento irá conter valores médios mais três produtos
médios das velocidades de flutuação. Assim, a equação de quantidade de movimento fica:
ρdu
d t=−∂p
∂x+ρgx + ∂
∂x
(µ∂u
∂x−ρu′2
)+ ∂
∂y
(µ∂u
∂y−ρu′v ′
)+ ∂
∂z
(µ∂u
∂z−ρu′w ′
)(2.78)
ρd v
d t=−∂p
∂y+ρg y + ∂
∂y
(µ∂v
∂y−ρv ′2
)+ ∂
∂x
(µ∂v
∂x−ρv ′u′
)+ ∂
∂z
(µ∂v
∂z−ρv ′w ′
)(2.79)
ρd w
d t=−∂p
∂z+ρgz + ∂
∂z
(µ∂w
∂z−ρw ′2
)+ ∂
∂y
(µ∂w
∂y−ρw ′v ′
)+ ∂
∂x
(µ∂w
∂x−ρw ′u′
)(2.80)
Nas quais os termos −ρu′i
2, −ρu′i u′
j e −ρu′i u′
k são os termos de aceleração convectiva
porém são chamados de termos de tensões turbulentas, pois tem dimensão de tensão e
aparecem ao lado dos termos de tensão newtoniana [27]. Fazendo a mudança de notação:
ui =Ui , u′i = 0 , p = P , p ′ = 0 (2.81)
As equações podem ser reduzidas para:
∂ui
∂xi= 0 (2.82)
ρ∂Ui
∂t+ρ ∂
∂x j
(UiU j
)=− ∂P
∂xi+ ∂
∂x j
(2µSi j −ρU ′
iU ′j
)(2.83)
38
Onde Si j é o tensor de deformações médio definido por :
Si j = 1
2
(∂Ui
∂x j+ ∂U j
∂xi
)(2.84)
O termo τi j =−U ′iU ′
j é conhecido como tensor de tensão de Reynolds e possui duas
componentes para cada direção.
Reescrevendo as Equações 2.82:
∂Ui
∂t+U j
∂Ui
∂x j=− ∂P
∂xi+ν ∂2Ui
∂xi∂x j− ∂U ′
iU ′j
∂x j(2.85)
Pode-se observar que ao utilizar o conceito de média de Reynolds, foram adicionados três
incógnitas às equações, gerando assim um problema de fechamento para a determinação das
variáveis do escoamento.
Para resolver essas novas equações, foram criados diversos modelos de turbulência que
buscam a solução das equações de média de Reynolds.
2.7.2 Lei da Parede
As tensões turbulentas mencionadas anteriormente em 2.7 são geralmente desconhecidas e
devem ser relacionadas experimentalmente à geometria e às condições de escoamento. Porém,
para escoamento em dutos e em camada limite desenvolvida, o termo da tensão −ρu′v ′
associada à direção y normal à parede é dominante, obtendo então uma forma aproximada
mais simples e com excelente precisão para o momento na direção principal [27].
ρdu
dt≈−∂p
∂x+ρgx + ∂τ
∂y(2.86)
Na qual a tensão cisalhante é definida como:
τ=µ∂uy
∂y−ρu′v ′ = τl am +τtur b (2.87)
39
Fig. 2.17: Distribuições típicas de velocidade e tensões cisalhantes próximas à parede
Fonte: Retirado de White [27]
A Figura 2.17 mostra as distribuições das tensões cisalhantes e turbulentas típicas de
medições. A tensão laminar é dominante na região próxima à parede, também chamada de
subcamada viscosa, e a tensão turbulenta domina na camada externa. Existe também uma
camada intermediária na qual tanto a tensão laminar quanto a turbulenta são importantes.
Assim, em 1930, Prandtl deduziu para a camada junto à parede (subcamada viscosa) que
u deve ser independente da espessura da camada sob cisalhamento e utilizando a análise
dimensional, chegou a:
u = f(µ,τp ,ρ, y
)(2.88)
u+ = u
u∗ = F
(yu∗
ν
)u∗ =
(τp
ρ
)1/2
(2.89)
A Equação 2.88 é conhecida como lei de parede e u∗ é denominado velocidade de
atrito. Em 1933 Kármán demostrou que u na camada externa é independente da viscosidade
molecular, porém seu desvio em relação à velocidade de corrente U deve depender da
espessura da camada limite δ e outras propriedades. A Equação 2.90 mostra a avaliação da
dependência de u na camada externa:
(U −u)ext = g(δ, τp ,ρ , y
)(2.90)
40
Reescrevendo a equação acima chega-se à lei da diferença de velocidade para a camada
externa:U −u
u∗ =G( y
δ
)(2.91)
Já para a camada intermediária, C. B. Milikan mostrou em 1937 que tanto a lei da diferença
de velocidade para a camada externa quanto a lei de parede deveriam superpor-se na camada
intermediária e isso só poderia ser verdade se a velocidade na camada intermediária variasse
logaritmicamente com y:u
u∗ = 1
κln
(yu∗
ν
)+B (2.92)
Onde κ e B são constantes e a equação é chamada de lei logarítmica da camada interme-
diária. A Figura 2.18 mostra os perfis de velocidade com as camadas.
Fig. 2.18: Verificação experimental das leis das camadas interna, intermediária e externa
Fonte: Retirado de White [27]
Como observado, a lei de parede é única e segue um padrão viscoso linear:
u+ = u
u∗ =(
yu∗
ν
)= y+ (2.93)
41
Sendo válido desde a parede até em torno de y+ = 5, desviando-se a a partir daí para se
juntar com a lei logarítmica em torno de y+ = 30.
2.7.3 Funções de Parede
No ANSYS CFD® todos os modelos de turbulência são independentes do valor de y+,
porém a eficácia do resultado requer que a função de parede mais adequada seja utilizada
dependendo do nível de refinamento da malha na parede e das escalas relativas do escoamento.
Como as condições próximas a parede são normalmente previsíveis pelas leis de parede
citadas anteriormente, estas funções podem ser utilizadas nas regiões próximas da parede
para determinar as variáveis evitando a necessidade de malhas muito refinadas, levando a um
custo computacional mais baixo. A Figura 2.19 mostra como o uso de funções de parede pode
simplificar bastante a solução nestas regiões.
Fig. 2.19: Uso de funções de parede
Fonte: Retirado do Manual do Usuáro Ansys [43] e editado pelos autores
As funções de parede são somente válidas entre valores específicos de y+. Valores muito
altos de y+ farão com que o primeiro nó fique fora da região da camada limite impossibilitando
o uso das funções de parede enquanto valores muito baixos farão com que o primeiro nó fique
na região linear viscosa onde estas funções também não são válidas.
O uso da função de parede padrão (standard), para modelos de turbulência baseados
na variável ε consideram que a malha na camada limite está totalmente contida na região
42
governada pela lei logarítmica (y+ > 30). Para a maioria das aplicações industriais, porém,
costumam haver variações nas escalas geométricas e de velocidade, o que leva a valores
diferentes de y+ ao longo do domínio. Neste caso recomenda-se o uso da função scalable
wall function. Esta função faz com que o valor mínimo de y+ para toda a malha suba para
30, valor a partir do qual as funções de parede podem ser aplicadas. Quaisquer valores de y+
abaixo desse mínimo são ignorados.
Para modelos κ−ε o valor de y+ para que o uso das funções de parede seja válido deve
ser menor do que 300.
Para modelos κ−ω, a solução da camada limite é feita diretamente, sendo utilizadas as
funções de parede apenas para valores de y+ maiores que 2. Isto ocorre para que se possa
aproveitar a eficácia do método para malhas refinadas próximas a parede. Nesta situação
utiliza-se a função de parede automática que detecta as regiões onde é necessário ou não o
emprego das leis de parede.
Modelos de turbulência SST seguem a mesma linha do modelo κ−ω utilizando a função
de parede automática para detectar as regiões onde a camada limite será resolvida diretamente
pelo modelo.
2.8 Modelos de Turbulência
2.8.1 Modelos de Equação Zero
Modelo de Prandtl (1925)
Este modelo foi desenvolvido encima de uma hipótese de que a velocidade é proporcional
à escala de comprimento e à taxa de cisalhamento.
lm =∣∣∣∣∂u1
∂x2
∣∣∣∣ (2.94)
Relacionado a Equação 2.94 com a viscosidade turbulenta, tem-se que:
νt = lm2∣∣∣∣∂u1
∂x2
∣∣∣∣ (2.95)
43
Vantagens do modelo:
• Fácil implementação
• Pequeno custo computacional
• Ótimo para uso em escoamentos internos e externos simples
Desvantagens:
• Baixa precisão para escoamentos turbulentos com comprimento variável que tenha
separação de camada limite ou circulação
• Só calcula as propriedades médias do escoamento e da tensão cisalhante turbulenta
• Não capta bem os efeitos dos modelos de combustão, empuxo, etc.
2.8.2 Modelos de Uma equação
Resolve uma única equação de conservação (uma equação diferencial parcial) para a
viscosidade turbulenta. A equação é composta por termos convectivos e difusivos, além de
utilizar uma expressão para a produção e dissipação da viscosidade turbulenta νt . Geralmente
é utilizado com malhas não estruturadas em aplicações da indústrias aeroespacial.
2.8.3 Modelos de Duas Equações
O presente trabalho utiliza modelos de turbulência de duas equações, que são inclusive os
mais aplicados atualmente. Assim, estes modelos resolvem duas equações de conservação
para viscosidade turbulenta.
Modelo κ−ε
É um dos modelos de turbulência mais utilizados. Consiste na inclusão de duas equações
de transporte para modelar as propriedades turbulentas do escoamento, permitindo contabilizar
os efeitos de difusão e convecção de energia turbulenta.
44
As equações adicionadas estão relacionadas às variáveis κ e ε, onde κ é um parâmetro
associado à energia cinética presente nos vórtices e ε é um parâmetro associado à dissipação
de energia turbulenta.
O modelo κ−ε é composto pelas seguintes equações:
∂(ρκ)
∂t+ ∂(ρκui )
∂xi= ∂
∂x j
[(µ+ µt
σκ
)∂κ
∂x j
]+Pκ+Pb −ρε−YM +Sκ (2.96)
e
∂(ρε)
∂t+ ∂(ρεui )
∂xi= ∂
∂x j
[(µ+ µt
σε
)∂ε
∂x j
]+C1ε
ε
κ(Pκ+C3εPb)−C2ερ
ε2
κ+Sε (2.97)
A variável Pκ é dada por:
Pκ =µt S2 (2.98)
Onde S é o módulo do tensor taxa de deformação principal, definido como:
S =√
2Si j Si j (2.99)
A variável Pb é dada por:
Pb =βgiµt
Prt
∂T
∂xi(2.100)
Onde Prt é o número de Prandtl turbulento, gi é a componente do vetor aceleração
gravitacional na direção i e β é o coeficiente de expansão térmica, definido como:
β=− 1
ρ
(∂ρ
∂T
)p
(2.101)
O modelo utiliza as seguintes escalas:
• Escala de Comprimento
κt3/2
ε(2.102)
45
• Escala de Velocidade
κt1/2 (2.103)
Assim, a viscosidade turbulenta é dada por:
νt =Cµκt
2
ε(2.104)
As demais constantes desse modelo são: Cµ = 0,09, C1ε = 1,44, C2ε = 1,92, σκ = 1 e
σε = 1,3
Vantagens do modelo:
• Fácil implementação
• Boa estabilidade e convergência
• Razoável precisão para maioria dos escoamentos práticos
Desvantagens:
• Pouca precisão para escoamentos rotacionais e vórtices, escoamentos com grande
separação, jatos assimétricos, escoamentos não confinados e escoamentos totalmente
desenvolvidos em dutos não circulares
• Válido somente para escoamentos totalmente turbulentos
• Possui uma equação simplista demais da energia cinética viscosa ε
• Difícil obtenção de dados turbulentos
• Necessidade de tratamento por função de parede
Modelo κ−ω
Semelhante ao modelo κ−ε, este é um modelo no qual são adicionadas duas equações
diferenciais parciais de transporte, associadas às variáveis κ e ω que estão relacionadas à
energia cinética turbulenta e a taxa específica de dissipação de energia, respectivamente.
46
O modelo κ−ω é composto pelas seguintes equações:
∂κ
∂t+U j
∂κ
∂x j= τi j
∂Ui
∂x j−β∗κω+ ∂
∂x j
[(ν+σ∗νT )
∂κ
∂x j
](2.105)
e∂ω
∂t+U j
∂ω
∂x j=αω
κτi j
∂Ui
∂x j−βω2 + ∂
∂x j
[(ν+σνT )
∂ω
∂x j
](2.106)
A viscosidade cinemática dos vórtices νT é dada por:
νT = κ
ω(2.107)
Os coeficientes de fechamento e demais relações auxiliares são: α= 5
9; β= 3
40; β∗ = 9
100;
σ= 1
2; σ∗ = 1
2e ε=β∗ωκ.
Vantagens do modelo:
• Bom desempenho na análise próxima à parede
• Recomendável para escoamentos com camada limite complexa, sob gradientes de
pressão adversos e separação
Desvantagens:
• Efeitos de separação são tipicamente previstos excessivamente e muito cedo
• Requer alta resolução de malha próxima à parede (y+ < 2)
• Possui custo computacional mais alto que o modelo κ−ε
Shear Stress Tensor (SST)
Este modelo nada mais é do que uma combinação dos dois modelos anteriores, κ− ε e
κ−ω, de modo a otimizar a solução aproveitando-se das vantagens de cada um.
É usada a formulação do modelo κ−ω na região interna da camada limite próxima às
paredes, aproveitando-se da boa eficiência desse método nesta região. Na região de fluxo
47
livre onde o modelo κ−ω é sensível às propriedades turbulentas o modelo SST muda para a
formulação κ−ε adequando-se às condições da região.
O modelo SST é composto pelas seguintes equações:
∂κ
∂t+U j
∂κ
∂x j= Pκ−β∗κω+ ∂
∂x j
[(ν+σκνT )
∂κ
∂x j
](2.108)
e
∂ω
∂t+U j
∂ω
∂x j=αS2 −βω2 + ∂
∂x j
[(ν+σωνT )
∂ω
∂x j
]+2(1−F1)σω2
1
ω
∂κ
∂xi
∂ω
∂xi(2.109)
A viscosidade cinemática dos vórtices νT é dada por:
νT = α1κ
max(α1ω,SF2)(2.110)
Os coeficientes de fechamento e demais relações auxiliares são:
F1 = tanh
[[min
(max
( pκ
β∗ωy,
500ν
y2ω
),
4σω2κ
C Dκωy2
)]4]
(2.111)
F2 = tanh
[[max
(2pκ
β∗ωy,
500ν
y2ω
)]2]
(2.112)
Pκ = min
(τi j
∂Ui
∂x j,10β∗κω
)(2.113)
C Dκω = max
(2ρσω2
1
ω
∂κ
∂xi
∂ω
∂xi,10−10
)(2.114)
φ=φ1F1 +φ2(1−F1) (2.115)
48
α1 = 5
9α2 = 0,44 β1 = 3
40β2 = 0,0828 β∗ = 9
100
σκ1 = 0,85 σκ2 = 1 σω1 = 0,5 σω2 = 0,856
(2.116)
Este método tem como vantagens as mesmas citadas para os métodos κ−ε e κ−ω, tendo
como principal desvantagem seu custo computacional mais elevado.
Modelo RNG κ−ε
O modelo RNG usa métodos de grupos de renormalização (Re-Normalisation Group -
RNG) para renormalizar as equações de Navier-Stokes de modo a contabilizar efeitos que
ocorrem em escalas de movimento menores.
No modelo κ−ε a viscosidade cinemática é calculada para uma única escala de turbulência
específica, fazendo com que a difusão de energia seja calculada apenas naquela escala. Porém,
em modelos reais, todas as escalas contribuem para a difusão de energia tubulenta, portanto
o modelo RNG utiliza de técnicas matemáticas que permitem variar o termo de geração de
modo a contabilizar as diferentes escalas.
Para esse modelo, as equações que o compõem são as seguintes:
∂(ρκ)
∂t+ ∂(ρκui )
∂xi= ∂
∂x j
[(µ+ µt
σκ
)∂κ
∂x j
]+Pκ−ρε (2.117)
e∂(ρε)
∂t+ ∂(ρεui )
∂xi= ∂
∂x j
[(µ+ µt
σε
)∂ε
∂x j
]+C1ε
ε
κPκ−C∗
2ερε2
κ(2.118)
onde
C∗2ε =C2ε+
Cµη3(1−η/η0)
1+βη3(2.119)
e
η= Sκ
ε(2.120)
onde S é definido da mesma forma que na Equação 2.99.
As demais constantes do modelo são: Cµ = 0,0845; σκ =σε = 0,7194; Cε1 = 1,42; Cε2 =
49
1,68; η0 = 4,38 e β= 0,012.
Vantagens do modelo:
• Bom desempenho na análise de escoamentos rotacionais
• Recomendável para análises de escoamento de ar em domínios fechados (túnel de
vento)
Desvantagens:
• Para a maioria dos casos não apresenta uma melhora tão grande em relação ao modelo
κ−ε padrão que justifique seu uso
• Não prevê muito bem a evolução dos vórtices
Modelo do tensor de Reynolds
Fisicamente o mais completo pois leva em conta toda a parte anisotrópica tensorial do
escoamento, rotação, vórtices e grande taxas de deformação.
A desvantagem deste modelo é o alto custo computacional devido ao grande acoplamento
do momento e das equações de turbulência.
Usado em escoamentos de ciclones, combustão, escoamentos envolvendo separação da
camada limite, etc. Se o escoamento possuir elevado swirling, natureza 3D ou escoamento
rotativo, usar RSM (Reynolds Stress Model).
2.9 Modelagem
2.9.1 Métodos de Discretização na Análise Numérica Computacional
A maioria das equações de governo dos problemas físicos reais são expressados matemati-
camente por equações diferenciais parciais (EDPs) onde, em algumas situações, uma solução
analítica pode não ser encontrada para o problema de valor de contorno considerado. Assim,
quando não se pode determinar numa solução exata, seja pela impossibilidade de ser determi-
nada analiticamente ou pela não familiaridade com as técnicas matemáticas adequadas, surge
a necessidade de determinar uma função que seja uma aproximação adequada da solução.
50
Atualmente, diversos métodos numéricos são utilizados para aproximar soluções dessas
equações diferenciais parciais. Entre estes métodos podem-se citar: Método de Elementos
Finitos (MEF), Método das Diferenças Finitas (MDF) e Método dos Volumes Finitos (MVF).
No método de diferenças finitas (MDF) são utilizadas expressões algébricas para apro-
ximar cada termo do modelo matemático em cada nó da malha. Esta formulação é dada
através da série de taylor da função da derivada, já o método de elemento finitos (usado no
software ANSYS Workbench®) consiste na divisão do domínio do problema em subdomínios
de dimensões finitas no qual o conjunto de todos subdomínios é igual ao domínio original,
como evidenciado em Costa [44], e utiliza diferentes métodos numéricos que aproximam a
solução de problemas de valor de fronteira tanto por equações diferenciais ordinárias quanto
por equações diferenciais parciais dentro do subdomínio, podendo assim obter a solução exata
através da interpolação das soluções aproximadas.
Historicamente, o método de diferenças finitas sempre foi empregado na área de mecânica
dos fluidos, enquanto o método de elementos finitos é mais utilizado na área estrutural
na parte de problemas de elasticidade. Esses problemas são complementante diferentes do
ponto de vista físico, pois os escoamentos são altamente não lineares enquanto os problemas
envolvendo elasticidade não possuem os termos advectivos e assemelham-se com problemas
puramente difusivos de transferência de calor de característica linear. Porém, em problemas de
advecção dominantes, tanto o método de diferenças finitas quanto o método de elemento finitos
produziram grandes instabilidades, assim motivaram-se as pesquisas para o aprimoramento
do método de volumes finitos, onde equações aproximadas são obtidas através do balanço de
conservação no volume elementar [45].
Quando em meados de 1970 passou-se a ceder espaço do que era antes comum o uso
de sistemas coordenados ortogonais convencionais para o uso de sistemas coordenados
generalizados coincidentes com a fronteira do domínio, o método de volumes finitos começou
a resolver problemas também para geometrias irregulares, difundindo seu uso na área numérica
de fluidos [45].
Para os pacotes comerciais como a ANSYS®, o método de volumes finitos ainda é o mais
empregado para uso em âmbitos industriais. Sua preferência é em função da robustez, devido
51
às características conservativas do método [45]. Como o método de diferenças finitas (MDF)
e método de elementos finitos(MEF) não são conservativos a nível discreto e sim em pontos
de malha, o MVF se adapta melhor aos princípios conservativos dos fluidos.
Uma das vantagens do MVF em relações aos outros métodos é que como o método de
volume finitos ao criar suas equações aproximadas satisfaz os balanços de conservação em
volume elementares (nível discreto), todos os princípios de conservação serão satisfeitos
mesmo conferidos em uma malha bastante grosseira tornando o tempo de simulação menor
[45].
O método de volumes finitos basicamente se baseia em dividir o domínio computacional
em pequenos volumes de controle (volumes elementares) e integrar no espaço e tempo cada
equação de governo de cada volume elementar, resultando em uma equação discreta que
possui os mesmos princípios de conservação que estes volumes. Assim, a aproximação de cada
incógnita é obtida em pontos específicos do domínio, podendo então ter uma representação
completa do comportamento do escoamento [43].
No presente trabalho, foi utilizado o software comercial ANSYS CFX® 15.0, que utiliza
o método de volume finitos baseado em elementos que foi desenvolvido originalmente para
resolver escoamentos descritos pelas equações de Navier-Stokes.
O conceito inicial do método surgiu com Baliga e Patankar [46] para resolver as equações
de advecção-difusão, porém logo mais tarde foi estendido para resolução de problemas mais
gerais de mecânica dos fluidos e transferência de calor.
O método de discretização utilizado neste trabalho é o método de volumes finitos base-
ado em elementos (Element-based Control Volume Methods - EbFVM) usado pelo pacote
comercial ANSYS CFX®. Este método engloba principalmente malhas não estruturadas e
cria volumes de controle para o balanço de discretização triangulares ou paralelogramos para
2D, e de tetraedros e hexaedros para 3D [45].
No método de volumes finitos baseado em elementos (EbFVm), a própria malha é utilizada
para definir os volumes elementares que serão criados, conectando por planos todos os centros
dos elementos vizinhos com os correspondentes centros das suas arestas, formando assim um
poliedro. Esta abordagem apresenta maior precisão numérica do que o método tradicional
52
de volumes finitos, pois possui mais pontos de integração por volume de controle (24 para
elementos hexaédricos e 60 para elementos tetraédricos) [45, 47].
Após definido o método de aproximação das equações diferenciais parciais que governam
os movimentos dos fluidos, deve-se selecionar o método adequado para e resolução do sistema
de equações algébricas lineares originado a partir da aplicação do método numérico. Esta
etapa de seleção do método numérico é de extrema importância pois o tempo computacional
necessário para a solução de um determinado problema representa em torno de 60%−70% do
tempo da solução do sistema linear.
O métodos diretos trabalham com matrizes completas e necessitam de processos equiva-
lentes à inversão da matriz, porém, ao contrário dos métodos indiretos, não necessitam de uma
estimativa inicial. Como para o caso das análises de escoamentos as matrizes obtidas com
aplicação de métodos numéricos são geralmente bastantes esparsas e grandes em tamanho, as
operações realizadas para o processo de inversão de matrizes trabalham com os elementos
zeros da matriz. Assim, o esforço computacional é elevado e por se tratar muitas vezes de
equações diferenciais não lineares, a matriz de coeficiente do sistema linear algébrico deve ser
atualizada constantemente, portanto não faz sentido utilizar os métodos diretos para resolução
de sistemas lineares, preferindo-se assim para aplicações de fluidodinâmica computacional os
métodos iterativos [45].
Discretização da equação de transporte
Assumindo-se como φ uma propriedade qualquer de interesse, a equação geral de trans-
porte como é utilizada no software comercial é dada por:
∂(ρφ)
∂t︸ ︷︷ ︸Termo Transiente
+ ∂(ρu jφ)
∂x j︸ ︷︷ ︸Termo Convectivo
= ∂
∂x j
(Γφ
∂φ
∂x j
)︸ ︷︷ ︸Termo Difusivo
+ Sφ︸︷︷︸Termo Fonte
(2.121)
Sendo Γφ o coeficiente de difusão relativo à φ e o Sφ é o termo de geração de φ por
unidade de volume.
Pode-se obter as equações de governo variando os valores de φ,Γφ,Sφ como mostrado na
Tabela 2.2.
53
Tab. 2.2: Valores das propriedades φ,Γφ,Sφ
Equação φ Γφ Sφ
Continuidade 1 0 0
Momentum em x u µ Bx + ∂
∂x
(µ∂u
∂x− 2
3µ∇· (~u))+ pd y
(µ∂v
∂x
)+ ∂
∂z
(µ∂w
∂z
)− ∂P
∂x
Momentum em y v µ By + ∂
∂y
(µ∂u
∂y− 2
3µ∇· (~u))+ pd x
(µ∂v
∂y
)+ ∂
∂z
(µ∂w
∂z
)− ∂P
∂y
Momentum em z w µ Bz + ∂
∂z
(µ∂u
∂z− 2
3µ∇· (~u))+ pd x
(µ∂v
∂z
)+ ∂
∂y
(µ∂v
∂z
)− ∂P
∂z
Energia T KCp
1Cp
DPDt + µ
Cpφ
A partir das equações de conservação, a discretização espacial ocorre por meio do cálculo
dos valores discretos da variável escalar φ nos nós dos volumes de controles representado na
Figura 2.20. No entanto, os termos advectivos ou convectivos necessitam ser conhecidos, e
assim são interpolados por meio de funções de interpolação, como será explicado na seção de
aspectos numéricos 2.9.6.
Fig. 2.20: Representação do cálculo usando volumes de controle
Fonte: Retirado de Hurtado [49]
Esquemas de interpolação dos termos advectivos
O ANSYS-CFX® dispõe de vários esquemas de interpolação de termos advectivos como
o esquema de diferenças centrais (Central Differencing Scheme, CDS), esquema upwind de
54
primeira ordem (Upwind Difference Scheme, UDS), fator de mistura específico (Specified
Blend Factor) e o de alta resolução (High Resolution).
Todos esquemas possuem a mesma equação básica:
φpi =φup + (β∇φ · ∆r
)(2.122)
Onde φi p é o valor da propriedade no ponto de integração (integration point) e φup é o
valor do termo à montante (upstream), β é uma constante chamada Blend Factor e ∆−→r é o
vetor do vértice vizinho que aponta para o ponto de integração. O termo(βr · ∇φ)
é também
chamado de Numerical Advection Correction (NAC) que pode ser interpretado como uma
correção anti-difusiva aplicada ao esquema Upwind [47].
1. Esquema de diferenças centrais (Central Differencing Scheme, CDS)
Este esquema possui precisão de segunda ordem e tem problemas de oscilações numéri-
cas em problemas advectivos dominantes, como representados na Figura 2.21. O valor
do gradiente ∇φ é o valor do gradiente analisado localmente no elemento da malha, já
o valor de β é para a aproximação de primeira ordem igual a duas unidades.
Fig. 2.21: Oscilações numéricas
Fonte: Retirado de Eidt [50]
2. Esquema Upwind (Upwind Difference Scheme, UDS)
É um método de discretização muito usado para a resolução de equações diferenciais
55
hiperbólicas, sendo um esquema robusto porém introduz erros difusivos devido a
suavização dos altos gradientes como mostrado na Figura 2.22. Quando β é zero, o
segundo termo da equação se torna nulo, resultando no esquema Upwind de primeira
ordem e quando β é 1, obtêm-se o esquema Upwind de segunda ordem.
Fig. 2.22: Gráfico sobre a suavização dos altos gradientes
Fonte: Retirado de Eidt [50]
3. Fator de Mistura Específico (Specified Blend Factor)
Neste método faz-se uma variação de 0 ≤ β ≤ 1 e o termo ∇φ recebe a média dos
gradientes nos nós adjacentes. Assim, os erros causados pela suavização dos altos
gradientes são reduzidos, entretanto este método introduz oscilações numéricas ainda
menores que as oscilações causadas pelo método de diferenças centrais.
4. Esquema de alta resolução (High Resolution)
Utiliza-se de um esquema não linear para designar um valor de β em cada nó da malha
baseado em limitadores de fluxo, com o objetivo de evitar oscilações das regiões com
fortes gradientes. Nas demais regiões, faz-se o valor de β aproximar-se de 1 localmente
para assegurar precisão de segunda ordem e utiliza-se o valor do gradiente do volume
de controle à montante para o termo ∇φ. A metodologia usada para este método de
discretização numérica começa pela avaliação do fluxo advectivo usando os valores
de β e ∇φ pelo método Upwind e então calcula-se β a partir de um φmi n e φmax
56
em todos os nós adjacentes ao ponto nodal em questão. Depois, para cada ponto de
integração, calcula-se o φi p e atualiza-se o valor de β de modo a garantir seu valor
entre os extremos.
Após estes cálculos, o valor de β no ponto nodal inicial é utilizado como valor mínimo
permitido para o cálculo de todos os pontos de integração nas vizinhas deste nó. Assim, o
esquema de alta resolução gera soluções de segunda ordem e evita a oscilação numérica,
como evidenciado por [47], sendo o método mais utilizado pelo ANSYS CFX® e o
método escolhido para a utilização no presente trabalho por ter mais vantagens que os
demais métodos.
Métodos Numéricos Iterativos
Os métodos numéricos iterativos são um conjunto de técnicas usadas para resolver os
sistemas de equações lineares da forma A x = b. Assim, como os computadores trabalham
apenas de forma algébrica, ou seja, somando e subtraindo, deve-se então adequar o modelo
matemático proveniente do modelo físico que representa o problema real, em equações
algébricas para o uso da ferramenta computacional para obtenção das soluções numéricas
aproximadas.
Estes métodos iterativos são aqueles que requerem um chute inicial para começar o
processo de solução e são classificados como ponto a ponto, linha a linha ou plano a plano.
Um método ponto a ponto é um método direto se a malha tiver apenas um volume elementar.
Consequentemente, se o problema for unidimensional, o método linha a linha será um método
direto, e o plano a plano será direto para um problema bidimensional [45].
1. Método de Jacobi
É um método ponto a ponto que resolve o sistema linear, atualizando o valor das
variáveis anteriores a cada iteração. Possui lenta convergência e requer uma matriz com
dominância diagonal para que possa convergir [45].
2. Método de Gauss-Seidel
57
Método similar ao anterior, porém faz uso dos valores das variáveis calculadas anterior-
mente na mesma iteração, acelerando a convergência. Apesar de ter uma convergência
mais rápida que o método de jacobi, este método ainda possui os pontos fracos de um
método iterativo ponto a ponto [45].
3. Método das Sobrerelaxações Sucessivas - S.O.R.
Este método procura acelerar ainda mais o processo de convergência, aplicando-se um
coeficiente de sobrerrelaxação nos valores obtidos a partir do método de Gauss-Seidel,
avançando mais rapidamente para a solução, o que por outro lado em demasia pode
causar divergência [45].
Os métodos de convergência ponto a ponto são lentos na transmissão de informação
advinda das condições de contorno, por isso deve-se partir de condições de contornos “fortes”
como a condição tipo Dirichlet, sendo a dominância diagonal da matriz dos coeficientes
o principal requisito para convergência deste método iterativo [45]. Contudo, o programa
utilizado neste trabalho utiliza um método iterativo algébrico chamado Algebric Multigrid
(AMG) a partir de uma extrapolação usando malhas mais grosseiras que a adotada no problema,
para a suavização dos erros numéricos e otimização na obtenção de resultados.
2.9.2 Malha
Como está demonstrado anteriormente, as equações de governo para fluidos geralmente
não possuem soluções analíticas e são discretizadas em volumes de controle pequenos. Para
tal, em problemas de fluidodinâmica computacional, o domínio a ser analisado é subdividido
em domínios menores nos quais serão aplicadas e resolvidas as equações de governo já
discretizadas. Estes subdomínios gerados são chamados de elementos ou células (cells) e o
conjunto de todos eles é chamado de malha. A Figura 2.23 mostra um exemplo de malha 2D
obtida a partir de um domínio de geometria circular.
58
Fig. 2.23: Subdivisão de um domínio em elementos de malha
Fonte: Retirado do site [51] e editado pelos autores
Os elementos são as regiões que geralmente são utilizadas como volume de controle onde
é feito o balanço com as equações de conservação. Os vértices dos elementos, os nós, são os
pontos onde são calculadas as variáveis do problema [45].
Os métodos tipicamente utilizados para aproximar e resolver as equações de governo
dentro de cada elemento são os já citados Método de Elementos Finitos (MEF), Método das
Diferenças Finitas (MDF) e Método dos Volumes Finitos (MVF). Deve-se tomar cuidado
para que se escolha um método que possa garantir uma continuidade das soluções através das
interfaces comuns entre cada elemento de forma que as soluções aproximadas dentro de cada
subdivisão possam ser agrupadas para fornecer uma visão completa do comportamento do
fluido ao longo do domínio inteiro.
Para a escolha da malha apropriada, devem-se considerar os seguintes parâmetros:
• Método de solução
• Resultados que se busca obter
• Tempo de configuração
• Esforço Computacional
• Geometria do domínio
O método de solução deve ser considerado porque, como citado acima, deve garantir a
continuidade das soluções através das células da malha em quastão.
59
Dependendo dos resultados que se deseja obter em uma simulação, diferentes configura-
ções de malha podem ser utilizadas. Por exemplo, se em uma simulação é importante a análise
dos efeitos de parede, a malha naquela região deve ser mais fina para que se possa captar com
mais precisão estes efeitos. A Figura 2.24 mostra um caso deste tipo que é referente ao caso
que é analisado neste trabalho.
Fig. 2.24: Malha refinada nas paredes para captação de efeitos específicos
Fonte: Retirado de Migliorini [8]
O tempo de configuração diz respeito ao tempo hábil que se leva para se estruturar a
malha. Malhas não automáticas exigem do usuário processos muita vezes demorados para se
construir, principalmente quando se necessita alterar a malha em regiões específicas como no
caso citado anteriormente onde a malha é refinada em regiões próximas a parede.
O esforço computacional para se criar uma malha também deve ser levado em considera-
ção. Mesmo para problemas que exigem uma malha bem fina, muitas vezes o desempenho
computacional para computar tal malha é tão caro que a torna inviável. Deve-se sempre fazer
um balanço entre a precisão necessária de resultados e o custo computacional requerido para
obtê-los.
A geometria do domínio está totalmente atrelada aos tipos de malha a serem utilizados.
Estes são definidos na subsecção a seguir.
60
2.9.3 Tipos de Malha
Malha Estruturada
Uma malha estruturada é caracterizada pela regularidade da sua estrutura em duas ou
três dimensões. Esta coesão geométrica restringe a escolha dos elementos em quadriláteros
em duas dimensões (2D) ou hexaedros em três dimensões (3D). Esta malha geralmente é
criada de forma com que se alinha na direção do escoamento fazendo com que exista menor
difusão numérica e consequentemente maior acurácia de resultados. Possui também um maior
controle sobre a qualidade da malha e da contagem de células ajudando a captar todos efeitos
das condições de contorno. Em contrapartida, pela limitação geométrica de seus elementos,
esse tipo de malha não se adequa bem à geometrias mais complexas.
Malha Não Estruturada
Uma malha não estruturada é caracterizada pela irregularidade em sua estrutura podendo
também ser em duas ou três dimensões. Esta configuração irregular permite que os ele-
mentos de malha tenham uma maior flexibilidade de formas que podem assumir, sendo
então geralmente aplicada para geometria complexas, economizando assim, esforço do usuá-
rio em configurar a malha e possui como característica uma grande quantidade de células,
possibilitando um grande poder computacional para resolução numérica.
Malha Híbrida
Uma malha híbrida contém tanto configuração estruturada quanto configuração não
estruturada. Esta configuração híbrida é usada para se otimizar os resultados e esforço
computacional, aproveitando-se das vantagens de ambos os tipos de malha. Geralmente utiliza-
se malha estruturada em regiões de geometria regular e em que se deseja obter resultados
mais precisos enquanto regiões de geometria complexa ou não muito relevantes se utiliza a
malha não estruturada. A Figura 2.25 mostra um exemplo de malha na região da cauda de um
aerofólio. Na região próxima à parede onde ocorrem efeitos provenientes da camada limite
e ondas de choque é ultilizada uma malha estruturada que obtém resultados mais precisos
61
enquanto na região de fluxo livre é utilizada malha tetraédrica automática já que nessa região
não ocorre nenhum fenômeno relevante que se deseja analisar.
Fig. 2.25: Malha híbrida na cauda de um aerofólio
Fonte: Retirado do site [52]
2.9.4 Classificação de Elementos
As malhas também podem ser classificadas com base na dimensão e nos tipos de elementos
presentes. A malha pode ser gerada em 3D ou 2D dependendo do tipo de análise e resolução
necessária para o problema de interesse. Para malha em duas dimensões, os elementos são
geralmente triangulares (tris) ou quadriláteros (quad), já para três dimensões, os elementos
comuns são tetraédricos (tets), hexaédricos (hex elements), pirâmides (pyramids) e prismáticos
(prisms).
Como todos os elementos em uma malha 3D são visualizados como elementos em 2D,
há uma semelhança entre as malhas em três dimensões e duas dimensões considerando os
contornos. Muitos algoritmos começam gerando malha na superfície antes de preencher o
interior do domínio com a malha (tais algoritmos também são conhecidos como boundary to
interior algorithms).
A qualidade da malha desempenha um papel importante na precisão e estabilidade numé-
rica do problema simulado. Os parâmetros mais importantes para se avaliar a qualidade da
malha são:
• Distribuição de nós
62
• Smoothness
• Skewness
• Razão de Aspecto
Estes parâmetros são ilustrados na Figura 2.26 e descritos a seguir.
Fig. 2.26: Exemplos dos principais parâmetros de malha
Fonte: Retirado do site [53] e editado pelos autores
A distribuição de nós diz respeito a discretização do domínio em regiões de interesse,
normalmente deve-se utilizar o refinamento (uso de células menores) onde houver detalhes
geométricos ou altos gradientes de propriedade. Assim, a distribuição de nós deve ser adequada
para que a solução do problema de interesse tenha uma relativa independência do número
de nós. Assim, normalmente é feito um estudo de convergência de malha obtido a partir do
uso de malha de diversos graus de refinamento que a partir de um determinado valor de nós e
elementos a precisão se torna quase constante, obtendo assim sua “convergência”.
O parâmetro de smoothness indica o quão brusca é a mudança de tamanho entre elementos
de malha subsequentes. Quanto menos brusca for essa variação, ou seja, quanto maior o
smoothness, melhor será a qualidade da malha por representar transições mais suaves entre os
elementos, o que evita erros de interpolação.
O parâmetro skewness indica o quão oblíquos estão os elementos de malha. O valor do
skewness varia de zero a um, sendo zero equivalente a um elemento de malha reto, que é o
ideal, e o valor 1 equivalente a uma célula totalmente torcida, o que nunca deve estar presente
em uma malha de boa qualidade.
63
A razão de aspecto é a razão entre o maior e o menor lado de um elemento de malha.
Quanto maior esse valor, mais "‘achatada"’ é a célula, o que não é desejável por poder
acarretar erros de interpolação durante a solução.
2.9.5 Estratégias para Uso de Malha Estruturada
Embora simples dutos possam ser modelados utilizando-se apenas um simples bloco, a
maioria das geometrias encontradas em casos reais devem ser modeladas utilizando uma
estratégia de multi blocos sempre que possível.
Existem diversos tipos de metodologia que se pode adotar na estratégia multi blocos
implementado usando o ANSYS ICEM CFD® como:
• H-Grid
• O-grid
• C-Grid
Estes tipos são exemplificados na Figura 2.28 e descritos abaixo:
Fig. 2.27: Tipos usuais de malha estruturada
Fonte: Retirado do site [54] e editado pelos autores
O H-Grid representa qualquer malha estruturada que apresente elementos unicamente
quadrangulares. É o tipo mais simples de malha estruturada porém não se adapta bem à
geometrias curvas.
O O-Grid representa malhas estruturadas nas quais as linhas de pontos se arranjam de
tal forma que o primeiro ponto encontre o último, formando linhas de malha em formato de
64
O (circulares). É utilizado em torno de geometrias circulares como furos ou extremidades
abauladas.
No C-Grid as linhas de pontos se curvam, em forma de C, e depois seguem direções
paralelas. Entre as principais aplicações estão as malhas utilizadas para análise de aerofólios,
como foi mostrado na Figura 2.28.
2.9.6 Aspecto Numéricos
Como mencionado na Seção 2.9.1 o núcleo numérico do ANSYS CFX® é baseado no
Método de Volumes Finitos Baseado em Elemento (EbFVm) que gera um volume de controle
em torno de cada nó, possuindo então maiores pontos de integração por volume de controle e
consequentemente mais precisão numérica mantendo o caráter conservativo do método de
volumes finitos aplicado para tratamento de malhas não-estruturadas.
O CFX® resolve por padrão as equações de governo Navier-Stokes. É realizado o acopla-
mento pressão-velocidade, construindo um sistema não-linear para os campos de velocidades
e pressão e resolvendo-os de uma vez só. Este método de solução pode ser tanto baseado em
pressão (Pressure-based Solver) quanto em massa específica (Density-based Solver).
O método de resolução baseado em pressão obtém o campo de velocidades a partir das
equações de conservação de momentum e o campo de pressão é obtido através da equação
acoplada de pressão e velocidade. Já o método baseado em massa específica obtém o campo
de massa específica a partir da equação da continuidade e assim como o método baseado em
pressão, o campo de velocidades é obtido a partir das equações de conservação de momentum
e a pressão e provêm da solução da equação de estado.
O CFX® utiliza como padrão uma malha co-localizada, isto é, a localização das variáveis
na malha computacional, ou arranjo de variáveis, acontece no centro de cada volume de
controle. Esta escolha resulta na simplicidade da implementação do algoritmo computacional,
utilizando assim um único volume de controle para realizar todas as integrações necessárias
para a discretização do domínio.
Para a resolução dos sistemas de equações lineares proveniente da discretização do
domínio o programa utiliza o método algébrico Multigrid (AMG) otimizado com a técnica
65
Incomplete LowerUpper (ILU) [55], cuja abordagem se baseia na suavização dos erros de alta
frequência utilizando de algumas iterações de outros métodos iterativos e de uma extrapolação
entre malhas mais grossas e mais finas, resolvendo as equações para uma malha mais grosseira
e corrigindo com malhas cada vez mais refinadas como mostrado na Figura 2.28, evitando
assim os erros de alta frequência, que são os erros da ordem do espaçamento da malha adotada
(λ≈∆x).
Fig. 2.28: Extrapolação usada no método Multigrid
Fonte: Retirado de Fonseca [55]
No ANSYS CFX®, ao resolver as equações de conservação, é utilizado um procedimento
de marcha no tempo, chamado falso transiente. No caso específico de um regime transiente por
exemplo, pode-se configurar o tamanho do passo do tempo e no caso de regime permanente,
como a formulação do tempo é implícita, um passo de tempo grande é especificado de forma
a obter rapidamente a solução em regime permanente [47].
2.9.7 Análise de Erros
Quando é inserida a linguagem matemática para dentro do ambiente computacional,
visando a resolução de equações diferenciais e integrais, introduzem-se erros ao utilizar
aproximações numéricas.
Estas aproximações numéricas podem ser minimizadas para cada tipo de problema e seus
erros devem ser analisados para sua validação física, caso contrário o resultado pode não ser
adequado para o problema de interesse. Erros sempre vão existir, porém deve-se utilizar de
um critério para garantir uma precisão suficientemente boa.
66
Erro de Discretização
A resolução das equações governantes do problema de interesse são resolvidas para pontos
discretos do domínio, estes pontos, para o caso de volume finitos, são volumes de controle
elementares que formam todo este domínio. Assim, devido à incapacidade computacional de
discretizar o domínio em um número infinito de volumes, considera-se que haverá sempre
um erro pois a simulação resolve o sistema de equações para estes volumes discretos e para
o restante utiliza-se um método de interpolação como diferenças centrais, Upwind, Quick,
Power-Law, entre outros.
Estes métodos introduzem um erro associado, e estes erro diminui com o refinamento de
malha ou com o aumento da ordem do método de interpolação.
Erro de Linearização
As equações de Navier-Stokes possuem termos não-lineares representados pelo produto
de duas variáveis não conhecidas da equação de momentum 2.123:
(u∂u
∂x+ v
∂u
∂y+w
∂u
∂z
)︸ ︷︷ ︸
ter mos não l i near es
(2.123)
Para resolver o sistema algébrico de termos não lineares, supõem-se valores para as
velocidades e atualizam-se estes valores através de um método iterativo.
Considerando um termo em particular, u1, e considerando que o termo u em qualquer
região do espaço pode ser descrito por um valor “advinhado” (ug ) mais um valor de correção
(∆u), chega-se a seguinte equação:
u = ug +∆u (2.124)
E se houver uma função relativa a u, podem-se desenvolver os termos a partir de séries de
Taylor:
f (u) = f(ug +∆u
)= f
(ug
)+∆u f ′ (ug)+ ∆u2
2f ′′ (ug
)+ ter mos de al t a or dem
(2.125)
67
Para os valores das funções:
f (u) = u2 f(ug
)= ug2 f ′ (ug
)= 2ug f ′ (u) = 2u (2.126)
Substituindo na Equação 2.125, obtem-se a seguinte solução:
u2 = ug2 +∆u
(ug
)+Er r o (2.127)
O erro mostrado na Equação 2.127 é obtido pois o termo f ′′ (ug)
é desprezado da equação.
Este erro é minimizado à medida que o valor de ug se torna a solução real, mostrando a
importância do método iterativo para esta atualização do valor de u. Este método é também
chamado de Método Newton-Raphson.
Panorama Geral de Erros
Observou-se que os erros são gerados ao implementar o modelo matemático ao mundo
computacional. Assim, é responsabilidade do encarregado pela simulação avaliar tais erros
para garantir o máximo de precisão nos resultados.
Para isto, deve-se ter um cuidado especial na minimização dos erros, pois ao refinar a
malha buscando a diminuição dos erros por discretização, acaba-se dificultando a convergência
dos métodos iterativos usados na linearização, o que faz aumentar este tipo de erro. Assim,
busca-se um balanço entre estes erros para se obter a melhor solução possível.
3 – Metodologia
Este capítulo apresenta o procedimento realizado para a modelagem do problema de
interesse, sendo analisados os resultados das modelagens de turbulência utilizando os modelos
κ−ε, κ−ω, RNG κ−ε e SST , em um escoamento de ar que escoa numa folga de selo mecânico
tipo hole-pattern cujas características são analisadas e assumindo-se regime estacionário e
compressibilidade do gás.
Para a simulação dos modelos de turbulência, é escolhido o uso do software comercial
ANSYS CFX® 15.0. O método numérico usado neste trabalho para resolver as equações de
governo do problema é o método de volumes finitos baseado em elemento (EbFVm) [45],
que facilita a utilização de malhas não-estruturadas em coordenadas cartesianas, mantendo
também o caráter conservativo do método de volumes finitos.
É utilizado como método iterativo, para a resolução numérica das equações diferenciais
parciais, o método High resolution. A solução então será comparada com os resultados de [8],
utilizando o modelo de turbulência κ−ε nas mesmas condições de contorno.
3.1 Geometria
A geometria do objeto de estudo provêm da região entre o selo e o rotor na qual haverá
o escoamento representado pela Figura 3.1 com base no texto do artigo [8]. Como esta
geometria do selo hole-pattern é simetricamente periódica ao longo da direção circunferencial,
pode-se utilizar apenas um trecho de toda geometria para a análise. Assim como Migliorini
[8], é utilizado um trecho de 6,2 graus do selo.
68
69
Fig. 3.1: Geometria do domínio de estudo e trecho de 6.2º graus da geometria
Fonte: Retirado de Migliorini [8]
O fluxo do escoamento através do domínio é demonstrado na Figura 3.2, demonstrando as
zonas de recirculação que ocorrem dentro dos furos cilíndricos.
Fig. 3.2: Demonstração do escoamento no domínio analisado
Fonte: Produzido pelos autores
A modelagem deste trecho do selo mecânico é feita utilizando o Design Modeler®,
seguindo as dimensões utilizadas em [8], mostradas a seguir na tabela 3.1.
70
Tab. 3.1: Geometria do Selo
Parâmetros Valor
Raio do Selo 57.37 mm
Comprimento do Selo 86 mm
Folga Radial 0.20 mm
Razão da Área do Furo 0.685
Diâmetro do Furo 3.175 mm
Profundidade do Furo 3.302 mm
3.2 Condições de Contorno
Para a configuração das condições de contorno que descrevem as condições físicas do
problema real, foi utilizado o programa CFX Pre®.
Para as condições de contorno, foram usadas as mesmas do artigo do Migliorini [8],
apresentadas na tabela 3.2. A direção do escoamento segue o sentido da direção z:
Tab. 3.2: Configurações das Condições de Contorno
Condições de Contorno
Entrada Pressão Total 7 MPa
Saída Pressão Estática 3.15 MPa
Rotor Rotação 20200 rpm
Estator Velocidade Nula V=0
Lados do Selo Periódico Rotacional
3.3 Pré-Análise
É muito comum ao se realizar uma simulação, tentar sempre buscar o maior número
de informações possível para o melhor conhecimento do problema e da física associada,
assim podem-se fazer simplificações e escolher o modelo matemático mais adequado para o
problema de interesse. Desta forma, antes de começar uma simulação, é necessário definir
71
todas as características básicas do escoamento de forma a simplificá-lo de todas as maneiras
possíveis e garantido que a capacidade computacional disponível atenda a precisão preterida
e mantendo o caráter físico do problema em questão.
Nesta etapa, é utilizado uma hipótese do absurdo de modo a dar mais confiabilidade para
a especificação das propriedades acerca do escoamento analisado no presente trabalho.
O problema de interesse envolve um escoamento dado por um gradiente de pressão elevado
(são conhecidas a pressão de entrada e a pressão na saída), em cavidades com geometria
complexa com uma das paredes em movimento periódico.
Neste problema, em que é utilizado o ar como fluido de trabalho, é complicado avaliar
a natureza do escoamento. Por não se ter as variáveis como velocidade, comprimento de
entrada nem como o fluido irá se comportar devido as mudanças bruscas de geometria. Como
mencionado anteriormente, foram usadas as mesmas condições de contorno baseadas em [8].
No entanto, o uso de simplificações é de grande importância para a melhor compreensão e
para a verificação do modelo matemático. A priori, como a maioria dos escoamentos reais
são turbulentos e por se tratar de um escoamento em cavidade (geometria complexa) com alto
gradiente de pressão, a hipótese de escoamento turbulento parece razoável.
Porém, para melhor entendimento físico do problema, pode-se realizar uma aproxima-
ção grosseira do problema para um duto simples com paredes estacionárias, isotérmico e
incompressível, o que fornece o escoamento de Poiseuille, cuja solução pode ser encontrada
analiticamente.
Foi utilizada a hipótese absurda do escoamento como incompressível em um tubo circular
reto de raio R. Assume-se também um escoamento totalmente desenvolvido, ou seja, a região
de interesse está suficientemente distante da entrada, de modo que o escoamento é puramente
axial.
Supondo também que há simetria axial,∂
∂θ= 0, a equação da continuidade para coordena-
das cilíndricas fica na forma:
∂vz
∂z= 0 ou vz = vz (r ) (3.1)
A equação de movimento linear (Navier-Stokes) fica simplificada a partir dos termos
72
∂p
∂r= 0, ou p = p (r ). Assim, a equação da quantidade de movimento em z reduz-se a:
ρvz∂vz
∂z=−dp
dz+µ52vz =−dp
dz+ µ
r
d
dr
(r
dvz
dr
)(3.2)
Como visto anteriormente em 3.1, tem-se que∂
∂z (vz) = 0, logo:
µ
r
d
dr
(r
dvz
dr
)= dp
dz= const ante (3.3)
Integrando duas vezes a equação 3.3, chega-se a:
vz = dp
dz
r 2
4µ+C1l n(r )+C2 (3.4)
Onde C1 e C2 são constantes. Empregando agora as condições de contorno devido à
geometria e à viscosidade do fluido:
1. Condição de não escorregamento em r = R: vz(R) = 0 = dpdz
R2
4µ +C1ln(R)+C2
2. Velocidade máxima finita (Para evitar a singularidade logarítmica e a respeito das
condições físicas do problema) em r = 0: vz(0) = 0 =C1ln(0)+C2
Assim, chega-se à conhecida solução analítica para o escoamento de Hagen-Poiseuille
totalmente desenvolvido:
vz =(−dp
dz
)1
4µ
(R2 − r 2) (3.5)
Para este escoamento, a velocidade média é dada por:
Vméd = 1
A
∫vz dA =
(−dp
dz
)R2
8µ(3.6)
Sabendo-se a aproximação de que:(−dp
dz
)= ∆p
L , o valor médio da geometria complexa
deve ser aproximado para um Rmédi o e uma viscosidade µ. Assim, para µ do ar, utiliza-se a
Lei de Sutherland, para o escoamento isotérmico com a temperatura de contorno 17,4°C dada
73
por:
µ
µ0≈
(T
T0
) 23(
T0 +S
T +S
)onde Sar ≈ 110,4K T0 = 273K ; µ0 = 1,71x10−5
[kg
m s
](3.7)
Resultando em:
µ≈ 1,70455x10−5[
kg
m s
](3.8)
Utilizando uma aproximação grosseira, visando a maior área possível, a largura do modelo
(na direção x) é transformada em um diâmetro, assim área média será dada por:
R = ∆x
2= 0,00318579[m] Améd =πR2 = 3,18848x10−5 [
m2] (3.9)
Tem-se então que a massa específica pode ser aproximada para a massa específica de
entrada (pois está simplificado para um escoamento incompressível):
ρ = P
RTonde R ≈ 287
[J
kg . K
](3.10)
Resultando em ρ≈ 84[
kgm3
]. Assim, a velocidade média do escoamento onde L = 0,082530871 m
é descrita por:
Vméd ≈(−∆p
∆L
)R2
8µ≈ 3471992,572
[m
s
](3.11)
Aparentemente, o valor da velocidade média se encontra elevado demais, o que já alerta
para algo fora do comum. A partir da referência [27], obtem-se a validação deste modelo para
que o escoamento seja laminar, caso o número adimensional de Reynolds do escoamento seja
ReD < 2100. Verificando o número de Reynolds:
Re = ρVméd (2R)
µ≈ 1,09002x1011 (3.12)
Este resultado é um caso clássico da hipótese do absurdo, na qual dentro de duas hipóteses,
escolhe-se uma a priori, e desenvolve-se um raciocínio acerca desta, tentando obter uma
validação para o caso. Caso contrário, assume-se a outra hipótese como a correta. Assim,
como o resultado da Equação 3.12 invalida a primeira hipótese de escoamento laminar,
74
assume-se a hipótese do escoamento ser turbulento (o que é razoável pois a maioria dos
escoamentos reais são de natureza turbulenta).
3.4 Discretização do Domínio
Para o estudo do escoamento em torno da cavidade do selo mecânico, deve-se primeiro
criar uma geometria que descreva as principais características dimensionais do problema.
Após a definição da geometria, como mostra a seção anterior, faz-se necessário a criação
de uma malha que descreva o mais precisamente possível a geometria do problema e as
regiões que necessitam de uma melhor precisão.
Nesta etapa do pré-processamento é feita a discretização do domínio computacional a
partir da geração de malha.
As malhas usadas neste trabalho foram criadas utilizando o programa Pointwise® 17R4,
onde foram obtidas malhas predominantemente hexaédricas (M1 e M2), o que de acordo com
[56] garantem maior precisão quanto maior a quantidade de elementos hexaédricos na malha.
Estas malhas não foram desenvolvidas pelos autores tendo sido fornecidas por colaboradores
do projeto da PETROBRAS "‘Desenvolvimento de metodologias e ferramentas numéricas
para a obtenção de coeficientes dinâmicos de selos internos de compressores centrífugos"’
(processo 2012/00251-0). Assim, as características das malhas M1 e M2 como números de
nós e elementos são mostradas na tabela 3.3:
Tab. 3.3: Discretização do domínio computacional
Malha Nós Elementos
M1 122152 111192
M2 521038 502208
Não apenas o tipo de elemento poliédrico faz diferença na precisão, como também o
refinamento nas regiões de interesse. Como é possível ver nas Figuras 3.3 e 3.4, as regiões
próximas à parede e no centro dos cilindros estão mais refinadas, pois representam regiões de
interesse do problema.
75
Fig. 3.3: Visualização da malha próxima à parede - detalhe 1
Fonte: Produzido pelos autores
Fig. 3.4: Visualização da malha nos centros dos cilindros - detalhe 2
Fonte: Produzido pelos autores
Após fazer uma rápida análise física do escoamento estudado, é uma boa prática sempre
tentar adequar ao máximo a modelagem ao problema físico, visando reduzir custo computaci-
onal e garantir maior precisão dos resultados.
76
3.5 Pré-Processamento
Nesta seção são discutidas todas as decisões e considerações realizadas na abordagem
do problema de interesse. As considerações realizadas para adequar o modelo matemático
ao problema físico do problema em um ambiente computacional proporcionado pelo uso do
software em questão são mostradas, assim como a natureza do escoamento, as condições de
contorno, os métodos iterativos para resolução do sistema de equações, o fluido a ser simulado
e os critérios de convergência que serão usados pelo Solver do ANSYS CFX® a partir da
configuração através da interface CFX-Pre®.
Com o intuito de estudar os efeitos dos modelos de turbulência na precisão da captura dos
resultados dos efeitos rotodinâmicos em um escoamento em um trecho de selo mecânico tipo
hole-pattern, desprezam-se os efeitos de qualquer excentricidade do rotor em relação ao eixo,
de maneira a garantir a simetria da geometria possibilitando a análise simplificada em um
trecho de selo de 6,2°, garantindo menor custo computacional e ainda mantendo o objetivo do
trabalho.
Nesta primeira simulação é resolvido um escoamento em regime estacionário, supondo
que a turbomáquina que contém o selo mecânico em questão se encontra funcionando a um
período de tempo suficientemente grande para que se possa considerar o escoamento como
estacionário. Baseado também no texto de [8], pode-se prosseguir considerando o escoamento
baseado em regime permanente.
A partir da pré-analise, onde é estimado que a hipótese de assumir o escoamento como
turbulento se parece razoável e utilizando como base o artigo do experimento realizado na
universidade do Texas [57], chega-se a uma primeira configuração física do problema.
As configurações foram ajustadas utilizando o pre-solver do CFX-Pre® onde primeira-
mente o regime estacionário foi ajustado de modo a se assemelhar com os trabalhos anteriores
como Migliorini [8] e Yan [12].
A Figura 3.5 apresenta a janela de seleção mostrando a configuração básica escolhida.
77
Fig. 3.5: Configuração do Regime Estacionário
Fonte: Produzido pelos autores
Durante a criação da malha, foi criado um corpo que seria exatamente o domínio de solução
do escoamento. Este mesmo corpo, é selecionado como domínio de trabalho na interface
CFX-Pre® definido como “domínio de fluido” garantindo que a região de escoamento do
fluido de interesse seja representada pela própria geometria.
O fluido de interesse foi configurado como ar ideal mostrado na Figura 3.6. Foi configurada
a pressão de referência como uma atmosfera e selecionada a morfologia do escoamento como
continuous fluid para usar uma abordagem euleriana para o tratamento do fluido.
Fig. 3.6: Configuração do Fluido de trabalho
Fonte: Produzido pelos autores
A seguir, para a primeira simulação, é assumido o uso para a primeira configuração
baseada em um modelo de turbulência κ− ε, onde foi contabilizada a equação geral da
energia incluindo o efeito de dissipação viscosa para tentar captar todos os efeitos físicos do
escoamento como mostrado na Figura 3.7.
78
Fig. 3.7: Configuração do modelo de turbulência e energia
Fonte: Produzido pelos autores
3.5.1 Condições de Contorno
Após definir a base principal para o solver, definem-se as configurações que representam
fisicamente as condições de contorno mostradas anteriormente na Tabela 3.2.
Para o uso das determinações de contorno, existem algumas configurações recomendadas
pelo manual de modelagem do ANSYS (Modelling Guide) na qual são estabelececidas as
condições de contorno que são normalmente usadas no ANSYS CFX e listadas em ordem de
mais robusto para menos robusto, como é mostrado a seguir.
• Mais robusto: Uso da condição de contorno de velocidade ou vazão na entrada e
pressão estática na saída. A pressão total na entrada é um resultado implícito das
condições de contorno anteriores.
• Robusto: Uso da pressão total na entrada e velocidade ou vazão na saída. A pressão
estática na saída e a velocidade na entrada são parte da solução.
• Sensível ao chute inicial: Pressão total na entrada e pressão estática na saída. A vazão
mássica do sistema é parte da solução.
• Não-confiável: Pressão estática na entrada e pressão estática na saída. Esta combinação
não é recomendada, pois como a pressão total e a vazão mássica são resultados implíci-
tos da predição dos resultados, resulta em condições de contorno que dão margem para
diversas possíveis soluções.
• Não recomendado: A pressão total não pode ser especificada para à saída pois é uma
condição de contorno instável para o escoamento que sai do domínio de solução
79
Assim, devido às condições de processamento do solver, as condições de contorno utiliza-
das baseadas em [8], são sensíveis ao chute inicial, devendo então haver muito cuidado ao
configurá-las para garantir sua estabilidade.
A configuração da pressão total descrita na região de entrada do domínio é feita selecio-
nando a região que representa o inlet e assumindo a condição de contorno essencial como
7 MPa, como mostrado na Figura 3.8. A direção do escoamento é normal à região de entrada
e a intensidade de turbulência foi utilizada a do padrão do ANSYS®.
(a) Condição de contorno usada na entrada (b) Região da entrada
Fig. 3.8: Configuração da Condições de Contorno na Entrada
Fonte: Produzido pelos autores
Já a configuração na saída do trecho do selo foi especificada para ter uma pressão estática
de 3.15 MPa e para ser tipo opening, possibilitando fluxos tanto entrando quanto saindo,
assim como mostrado em na Figura 3.9.
(a) Condição de contorno usada na saída (b) Região da saída
Fig. 3.9: Configuração da condições de contorno na saída
Fonte: Produzido pelos autores
As condições de contorno para o rotor e o estator se baseiam na mesma usada pela
Universidade do Texas [57], apresentadas também na Tabela 3.2 e mostradas nas Figuras 3.10
80
e 3.11 para um rotor concêntrico ao eixo, girando a 20.200 rotações por minuto.
(a) Condição de contorno usada no estator (b) Região do estator
Fig. 3.10: Configuração da região do estator
Fonte: Produzido pelos autores
(a) Condição de contorno usada no rotor (b) Região do rotor
Fig. 3.11: Configuração das condições de contorno no rotor
Fonte: Produzido pelos autores
Como foi utilizado um trecho do selo total, as condições de contorno para os lados desta
geometria deveriam refletir as condições físicas, como a possibilidade de entrada e saída e
também o movimento rotacional do rotor. A partir de outros trabalhos semelhantes, utilizou-se
a condição de periodicidade rotacional para cada interface, como mostrado na Figura 3.12:
81
(a) Condição de contorno usada nas laterais (b) Região lateral da geometria
Fig. 3.12: Configuração das condições de contorno da interface lateral da geometria
Fonte: Produzido pelos autores
3.5.2 Controle do Processamento
Na etapa de processamento, o ANSYS CFX® permite a visualização e configuração
de diversos parâmetros, entre os quais os resíduos das variáveis durante o transcorrer das
iterações e seus respectivos critérios de parada. Para o caso da primeira simulação, onde
há maior interesse em verificar a validade do caráter físico do problema, foi utilizado um
esquemático mostrado na Figura 3.13.
82
Fig. 3.13: Configuração do controle de processamento - Solver Control
Fonte: Produzido pelos autores
O ANSYS CFX® recomenda o uso do High resolution como esquemático para os termos
advectivos para a maioria dos casos, pois possui precisão de segunda ordem e é bastante
adaptável nas situações de fortes gradientes.
A opção do uso de residuais máximo é justificada pela posição conservadora da busca por
precisão, pois os residuais máximos são os valores máximos para cada iteração a partir das
maiores variações para todos os pontos da malha. Assim, o critério de parada foi estabelecido
como resíduo máximo de 0,001 ou 100 000 iterações.
O programa executa as simulações em precisão simples (ponto flutuante com uso de 4
bytes) ou por dupla precisão (uso de 8 bytes). Como o presente estudo irá analisar apenas
um trecho da geometria e sendo um escoamento monofásico com limitações de hardware
com o objetivo de analisar a priori apenas as diferenças dos modelos de turbulência, opta-se
por utilizar precisão simples. As simulações foram executadas em um sistema operacional
Windows 10 com processador de 64-bit i7-4500U utilizando 8GB de RAM.
É importante avaliar a quantificação das incertezas, pois é um aspecto fundamental
para a utilização segura das ferramentas de CFD. Surge assim a necessidade da verificação e
validação cuidadosa dos resultados para verificação de inconsistências físicas que invalidariam
a solução numérica.
Para a validação física, têm-se a necessidade de observar se as condições de contorno se
83
mantiveram na resolução do problema, por exemplo, se as velocidades obedeceram o princípio
de não escorregamento nas superfícies dos corpos, e para os casos de regime permanente, em
que há saídas e entradas de fluido, uma das verificações importantes a serem realizadas é o
cálculo das vazões em massa das entradas que devem ter valores rigorosamente próximos aos
das saídas. Assim, os procedimentos usados para a validação numérica de cada simulação
incluem a confirmação de que a massa é conservada, sendo esta análise inclusa nos monitores
do CFX-Solver Manager® 15.0 a partir da criação de um monitor de vazão de entrada e vazão
de saída mostrado na Figura 3.14:
Fig. 3.14: Monitoração dos parâmetros físicos do solver
Fonte: Produzido pelos autores
O pacote Ansys Workbench® 15.0 oferece uma interface de pós-processamento chamada
CFD-Post® 15.0, assim pode-se fazer a visualização de parâmetros do escoamento através de
imagens, gráficos, tabelas e animações.
Os parâmetros analisados são as pressões médias que serão obtidas por meio de planos
transversais que distanciam-se 0,003614 m um do outro, apresentados no Apêndice A 6, o
vazamento do selo mecânico que é representado pela vazão mássica na saída do selo, obtido
pelo uso da formula =ave(MassFlow)@Outlet dentro do próprio CFX-Post®, e pelos valores
de tensão cisalhantes, obtidas pelo opção de inserir imagem de contorno na região do rotor.
84
Fig. 3.15: Planos transversais na geometria
Fonte: Produzido pelos autores
Utiliza-se também a ferramenta de pós processamento para avaliar os resultados e averiguar
as inconsistências físicas que invalidariam a solução numérica.
4 – Resultados e Discussão
Este capítulo apresenta os resultados da validação do modelo isotérmico e das propriedades
obtidas utilizando a mesma configuração mostrada no capítulo anterior de metodologia, para
os diferentes modelos de turbulência nas malhas M1 e M2 que são então comparados com os
valores obtidos por trabalhos anteriores para sua validação e análise.
4.0.1 Verificação e Validação Preliminar
Nesta subseção, são analisados e verificados os resultados da simulação utilizando o
modelo κ−ε com a equação da energia completa utilizando a malha M1 com o objetivo de
validar o modelo numérico a partir da comparação dos resultados mostrados pelo CFD-Post®
15.0 com os dados obtidos experimentalmente [12, 57, 58].
Na pré-análise da validação do modelo isotérmico, não foi obtida a convergência a partir
dos critérios estabelecidos anteriormente, assim foi realizada uma análise do caráter físico para
a simulação envolvendo o modelo κ−ε utilizando a equação completa de energia considerando
o termo relativo à dissipação viscosa.
A figura 4.1 obtida no CFX-Solver Manager® 15.0 mostra a variação dos residuais com
as iterações, podendo-se observar que a maioria dos residuais atenderam o critério de residual
máximo, porém o maior resíduo final para este caso foi a velocidade na direção z, W −Mom,
estabilizada na ordem de 10−2.
85
86
Fig. 4.1: Residuais do modelo κ−ε com equação da energia completa
Fonte: Produzido pelos autores
O grande valor obtido pelo residual W −Mom pode ser explicado em parte pela grande va-
riação de pressão e pela mudança abrupta de geometria proporcionado pelas formas cilíndricas,
gerando uma perturbação instável na direção z do escoamento.
A solução analisada foi obtida com 74 905 iterações após 129,34 horas de simulação
utilizando dois processadores em paralelo.
Um dos parâmetros mais importantes na comparação com o problema experimental é a
vedação obtida por meio do uso do selo mecânico tipo Hole Pattern. Esta análise da vedação
é realizada através do valor da vazão na saída da geometria do problema.
Foi obtido o valor de vazão mássica 0.00657444[ kgs ] e ao multiplicar por 58 seções do
selo obtem-se a vazão mássica total de 0,381318[ kgs ]. Ao analisar o valor obtido pela vazão
mássica na saída e comparado com o valor obtido pelo experimento da Universidade do Texas
[1] é visível que, mesmo não atendendo os critérios de convergência, foi obtido um resultado
razoável com 6,77% de erro comparado ao valor experimental de 0,409[ kgs ].
A figura 4.2 mostra a variação total de temperatura ao longo da geometria do problema,
que demonstra a obtenção de uma variação de 5,73% de temperatura, o que é bem próximo
dos valores obtidos experimentalmente em [12, 57, 58].
87
A partir da similaridade com os dados experimentais, valida-se o modelo utilizado e a
simplificação do uso do modelo isotérmico que já é utilizada por [8] para a redução do custo
computacional.
Fig. 4.2: Variação total de temperatura do modelo κ−ε com equação da energia completa
Fonte: Produzido pelos autores
A Figura 4.2 mostra o descréscimo da temperatura ao longo da direção axial da seção do
domínio, que ocorre devido às sucessivas expansões que o fluido sofre e também é influenciado
pelas forças de atrito.
Desta etapa adiante, é utilizada a simplificação do modelo de energia isotérmico, ou seja,
não é usada a equação de energia para a resolução do escoamento.
4.0.2 Resultados dos Modelos de Turbulência
Para a validação numérica, foi configurada no pós-processamento a visualização da
pressão máxima de 7 MPa, assim os valores de pressão acima são considerados como erros
gerados numericamente.
Foram utilizadas as mesmas configurações básicas, em ambas as malhas M1 e M2, na
análise de cada modelo de turbulência. Pode-se então analisar se a configuração imposta
representa corretamente o caráter físico do problema. Assim, pode ser visto a partir do
resultado mostrado na Figura 4.3 da primeira simulação para o modelo κ− ε na malha
M1 que as condições de contorno físicas do problema, como a velocidade de rotação do
88
rotor, introduzem um movimento rotacional nas linhas de corrente deste problema o que era
esperado.
Fig. 4.3: Visualização das linhas de correntes no modelo κ−εFonte: Produzido pelos autores
Também é feita uma análise da variação do tamanho dos elementos ao longo de cada
malha. As figuras 4.4a e 4.4b mostram essa variação a partir da determinação dos valores de
y+ ao longo de cada malha. Estes dados são muito importantes para a análise do desempenho
de certos modelos de turbulência.
(a) Visualização Y-Plus da malha M1 (b) Visualização Y-Plus da malha M2
Fig. 4.4: Visualização Y-Plus em ambas malhas M1 e M2
Fonte: Produzido pelos autores
89
Pode-se ver também que a simulação possui erros numéricos, que são causadas pelas
condições de contorno impostas em elementos interpolados da malha, assim é visto na Figura
4.5 um pequeno salto da pressão na região próxima da entrada do escoamento, que foi imposto
para ter uma pressão total especificada e bem ao lado foi configurada uma condição periódica
de rotação o que possivelmente gerou a instabilidade vista na figura.
Fig. 4.5: Erro numérico apresentado pelo modelo κ−ε na região de entrada na malha M1
Fonte: Produzido pelos autores
A Figura 4.6 mostra que os vetores de velocidade são calculadas em cada nó da malha,
o que demonstra uma característica intrínseca do programa comercial ANSYS CFX® que
resolve os balanços do método de volume finitos nos elementos, ou seja, são baseados nos
nós da malha como visto na figura.
90
Fig. 4.6: Visualização dos vetores de velocidade na parte lateral do selo - resultado obtido apartir do modelo SST na malha M2
Fonte: Produzido pelos autores
A Figura 4.7 evidencia uma característica esperada para o escoamento, pode-se ver a
formação de recirculação proveniente da geometria do selo. Esta zona de circulação é marcada
fisicamente pela presenta de gradientes adversos e favoráveis como mostrado pela figura.
Fig. 4.7: Visualização dos vetores de velocidade na parte lateral do selo - resultado obtido apartir do modelo SST na malha M2
Fonte: Produzido pelos autores
Após uma breve analise dos resultados obtidos, as configurações numéricas usadas e os
resultados obtidos demostram ser razoavelmente adequados para o problema de interesse.
91
Assim, pode-se prosseguir para a simulação dos demais modelos de turbulência visando sua
comparação com o modelo experimental.
É mostrado nas Tabelas 4.1 e 4.2 o tempo computacional usado em cada modelo de
turbulência nas malhas M1 e M2, assim como o número de iterações realizadas e se a
convergência foi obtida, usando como critério a obtenção de residuais máximos da ordem de
10−3.
Tab. 4.1: Tabela dos parâmetros da resolução dos modelos de turbulência para malha M1
Modelo de Turbulência Tempo Computacional Numero de Iterações Convergiu?
κ−ε 11,57 horas 6590 Sim
κ−ω 0,4 horas 116 Sim
RNG κ−ε 54 horas 44682 Não
SST 12,78 horas 4400 Sim
Tab. 4.2: Tabela dos parâmetros da resolução dos modelos de turbulência para malha M2
Modelo de Turbulência Tempo Computacional Numero de Iterações Convergiu?
κ−ε 139 horas 21848 Sim
κ−ω 55 horas 5597 Sim
RNG κ−ε 69 horas 6683 Sim
SST 54 horas 8504 Sim
Estes resultados demostram que as simulações processadas na malha M2 chegam a levar
até 136 vezes mais tempo para a simulação utilizando o modelo κ−ω que convergiu excepci-
onalmente rápido na malha M1, contudo mesmo não considerando essa excepcionalidade, o
modelo κ−ε da malha M2 levou 12 vezes o tempo exigido para o mesmo modelo processado
na malha M1.
Pode-se também perceber que o número de iterações realizados na malha M2 também
é razoavelmente maior, porém a redução se torna muito mais branda do que a diferença em
termos de tempo computacional, tendo uma diferença de apenas 48 vezes mais iterações para
o modelo κ−ω da malha M2 e apenas três vezes mais para o modelo κ−ε da malha M2.
92
Assim os resultados obtidos confirmam a generalização de que uma malha mais refinada,
em geral, leva mais tempo de processamento (tempo computacional) pois a matriz do sistema
de equações lineares obtidos pela matriz mais refinada é sempre maior que a da obtida na
malha grosseira, contudo, devido ao menor uso de interpolação entre o domínio discreto
computacional, por ter mais nós e elementos, a convergência é obtida mais rapidamente.
A seguir são mostrados os residuais obtidos nas simulações realizadas no presente tra-
balho, demonstrando a estabilidade numérica dos modelos sendo processados na malha
correspondente e sua convergência para o passo de tempo acumulado.
4.0.3 Análise dos Residuais
Os residuais são mostrados nas figuras a seguir evidenciando uma estabilidade ou linha de
tendência atingida após um certo número de iterações.
Residuais do modelo κ−ε
A Figura 4.8a mostra que o modelo κ−ε para a malha M1 convergiu após cerca de 6700
iterações, já a Figura 4.8b mostra que no modelo κ− ε para a malha M2, as velocidades
atingiram o critério de parada após cerca de 12000 iterações, porém o modelo só convergiu
após derca de 22000 iterações como é mostrado anteriormente na Tabela 4.2.
(a) Residuais do κ−ε da malha M1 (b) Residuais do κ−ε da malha M2
Fig. 4.8: Residuais obtidos utilizando o modelo κ−ε em ambas malhas M1 e M2
Fonte: Produzido pelos autores
93
Residuais do modelo κ−ω
A figura 4.9a mostra que o modelo κ−ω para a malha M1 convergiu após cerca de 175
iterações, já para malha M2 a figura 4.9b mostra que o modelo κ−ω convergiu após cerca de
5600 iterações:
(a) Residuais do κ−ω da malha M1 (b) Residuais do κ−ω da malha M2
Fig. 4.9: Residuais obtidos utilizando o modelo κ−ω em ambas malhas M1 e M2
Fonte: Produzido pelos autores
Residuais do modelo RNG κ−ε
A simulação usando o modelo RNG κ−ε só obteve a convergência na malha M2. A figura
4.10a mostra que os residuais não atingiram o critério de convergência estabelecido de 10−3,
e a figura 4.10b mostra a convergência obtida para a malha M2 após cerca de 6600 iterações:
94
(a) Residuais do RNG κ−ε da malha M1 (b) Residuais do RNG κ−ε da malha M2
Fig. 4.10: Residuais obtidos utilizando o modelo RNG κ−ε em ambas malhas M1 e M2
Fonte: Produzido pelos autores
Residuais do modelo SST
A figura 4.11a mostra que para a malha M1 o modelo SST teve sua convergência após
cerca de 4400 iterações.
No caso da malha mais fina, M2, os residuais atingiram o critério de convergência com
cerca de 6000 iterações conforme mostrado na figura 4.11b, entretanto, o modelo em si só
convergiu após cerca de 8500 iterações como é indicado anteriormente pela tabela 4.2.
95
(a) Residuais do SST da malha M1 (b) Residuais do SST da malha M2
Fig. 4.11: Residuais obtidos utilizando o modelo SST em ambas malhas M1 e M2
Fonte: Produzido pelos autores
4.0.4 Tensões Cisalhantes
Um dos parâmetros importantes para a verificação das diferenças básicas na análise dos
diferentes modelos de turbulência é a tensão cisalhante no rotor, pois a maioria dos modelos
de turbulência foram formulados para avaliar o escoamento em uma determinada distância da
parede. As imagens de contorno para as malhas M1 e M2 são apresentadas para o modelo
κ−ε nas Figuras 4.12a e 4.12b, para o modelo κ−ω nas Figuras 4.13a e 4.13b, para o modelo
RNG κ−ε nas Figuras 4.14a e 4.14b e para o modelo SST nas Figuras 4.15a e 4.15b.
96
Imagem do contorno do modelo κ−ε
(a) Tensões cisalhantes modelo κ−ε da M1 (b) Tensões cisalhantes modelo κ−ε da M2
Fig. 4.12: Visualização das tensões cisalhantes na parede do rotor do modelo κ−ε em ambasmalhas M1 e M2
Fonte: Produzido pelos autores
Imagem do contorno do modelo κ−ω
(a) Tensões cisalhantes modelo κ−ω da M1 (b) Tensões cisalhantes modelo κ−ω da M2
Fig. 4.13: Visualização das tensões cisalhantes na parede do rotor do modelo κ−ω em ambasmalhas M1 e M2
Fonte: Produzido pelos autores
97
Imagem do contorno do modelo RNG κ−ε
(a) Tensões cisalhantes modelo RNG κ−ε da M1 (b) Tensões cisalhantes modelo RNG κ−ε da M2
Fig. 4.14: Visualização das tensões cisalhantes na parede do rotor do modelo RNG κ−ε emambas malhas M1 e M2
Fonte: Produzido pelos autores
Imagem do contorno do modelo SST
(a) Tensões cisalhantes modelo SST da M1 (b) Tensões cisalhantes modelo SST da M2
Fig. 4.15: Visualização das tensões cisalhantes na parede do rotor do modelo SST em ambasmalhas M1 e M2
Fonte: Produzido pelos autores
O modelo Shear Stress Transport utiliza-se de um fator de combinação, chamado também
de blend factor, para o uso de otimização do modelo de turbulência na adaptação da malha
usada para a geometria do problema. Assim, nas regiões mais próximas da parede este modelo
se comporta como o κ−ω e nas regiões mais distantes como o κ−ε.Ao utilizar este modelo, surge a possibilidade de visualização da imagem de contorno
deste fator de combinação, assim podem-se usar estas imagens para a verificação da malha e
do uso deste modelo de turbulência.
98
A seguir são mostradas as imagens de contorno dos blend factors (Figuras 4.17a e 4.17b)
e da wall scale (Figuras 4.16a e 4.16b).
Imagem do contorno das escalas de parede do modelo SST
(a) Escalas de parede modelo SST da M1 (b) Escalas de parede modelo SST da M2
Fig. 4.16: Visualização das escalas de parede do modelo SST em ambas malhas M1 e M2
Fonte: Produzido pelos autores
Imagem do contorno do first blend factor no modelo SST
(a) First blend factor modelo SST da M1 (b) First blend factor modelo SST da M2
Fig. 4.17: Visualização do uso do primeiro blend SST em ambas malhas M1 e M2
Fonte: Produzido pelos autores
99
Imagem do contorno do Second blend factor no modelo SST
(a) Second blend factor modelo SST da M1 (b) Second blend factor modelo SST da M2
Fig. 4.18: Visualização do uso do segundo blend SST em ambas malhas M1 e M2
Fonte: Produzido pelos autores
4.0.5 Comparações entre os Modelos de Turbulência
A seguir, são apresentadas a partir das figuras 4.19 e 4.20 e das tabelas 4.3, 4.4, 4.5 e 4.6 as
comparações diretas entre as variáveis de pressão total média obtida nos planos transversais,
entre as vazões mássicas obtidas comparadas ao valor obtido pelo modelo experimental e
entre as tensões cisalhantes nas direções x, y e z, respectivamente.
Fig. 4.19: Gráfico das pressões totais média obtidas nos planos transversais para a malha M1
Fonte: Produzido pelos autores
100
Fig. 4.20: Gráfico das pressões totais média obtidas nos planos transversais para a malha M2
Fonte: Produzido pelos autores
Tab. 4.3: Tabela das vazões obtidas em cada simulação para malha M1
Modelo de Turbulência Vazão Total Obtida (kgs ) Erro (%)
κ−ε 0,38311 6,329
κ−ω 0,37696 7,831
RNG κ−ε 0,53079 29,78
SST 0,39969 2,275
Experimental 0,409 0
Tab. 4.4: Tabela das vazões obtidas em cada simulação para malha M2
Modelo de Turbulência Vazão Total Obtida (kgs ) Erro (%)
κ−ε 0,37872 7,402
κ−ω 0,35757 12,572
RNG κ−ε 0,42157 3,075
SST 0,40153 1,824
Experimental 0,409 0
101
Tab. 4.5: Tensões Cisalhantes no rotor para a malha M1
Tensão Cisalhante no Rotor (Wall Shear)
Modelo de Turbulência Direção X Direção Y Direção Z Resultante
κ−ε 1645.12 Pa -0.269025 Pa 1678.23 Pa 2372.47 Pa
κ−ω 1589.99 Pa -0.193298 Pa 1457.82 Pa 2174.99 Pa
RNG κ−ε 1121.4 Pa 0.0234617 Pa 1889.52 Pa 2228.6 Pa
SST 1558.78 Pa -0.206694 Pa 1580.97 Pa 2239.69 Pa
Tab. 4.6: Tensões Cisalhantes no rotor para a malha M2
Tensão Cisalhante no Rotor (Wall Shear)
Modelo de Turbulência Direção X Direção Y Direção Z Resultante
κ−ε 1633.62 Pa -0.52437 Pa 1669.92 Pa 2358.59 Pa
κ−ω 1619.3 Pa -0.371526 Pa 1389.26 Pa 2150.66 Pa
RNG κ−ε 1409.56 Pa -0.375696 Pa 1741.08 Pa 2264.98 Pa
SST 1565.11 Pa -0.318568 Pa 1599.92 Pa 2257.17 Pa
5 – Conclusões
A partir dos resultados obtidos, chega-se a algumas conclusões a respeito do uso de dife-
rentes modelos de turbulência e malhas de diferentes níveis de refinamento para o problema
em questão.
Ao analisar os gráficos de distribuição de pressão obtidos para os modelos analisados
em ambas as malhas, pode-se ver que os resultados ficaram muito próximos, com exceção
do modelo RNG κ−ε para a malha M1, já que este não atingiu o critério de convergência.
Neste caso pode-se concluir que para fins de obtenção da distribuição de pressão, que é
um dado importante caso se queira calcular os coeficientes rotodinâmicos, o modelo de
turbulência utilizado não fará muita diferença, permitindo que se opte por aquele que apresenta
convergência mais rápida.
Em relação aos dados obtidos de tensão cisalhante no rotor, pode-se observar que os
valores na direção y são quase nulos, o que era esperado uma vez que a face do rotor está
praticamente toda contida no plano xz, apesar da leve curvatura. As intensidades da tensão
cisalhante nas outras direções para os modelos avaliados tiveram algumas diferenças, mas
as tensões resultantes para todos os casos nas duas malhas foram relativamente próximas,
o que demonstra que as funções de parede utilizadas pelos métodos κ− ε e RNG κ− ε se
aproximam bem dos resultados calculados diretamente pelos métodos κ−ω e SST . Desta
forma pode-se concluir que o refinamento próximo às paredes não é tão primordial para
a obtenção das tensões na situação analisada, já que os modelos que se beneficiam deste
refinamento apresentaram resultados próximos dos modelos que aproximam a região por
funções de parede.
A respeito dos valores de vazão mássica obtidos, os modelos SST apresentaram resultados
finais mais próximos do valor experimental, com erros relativos de 2,275% e 1,824% nas
malhas M1 e M2 respectivamente, o que era esperado uma vez que este modelo se comporta
bem em malhas com variações grandes de y+, que é o caso das malhas utilizadas neste
trabalho, como pode ser visto nas figuras 4.4a e 4.4b. O modelo RNG κ− ε para a malha
M1 apresentou o pior resultado, 29,78%, o que é explicado, em parte, pelo fato deste não
ter atingido o critério de convergência. O modelo κ−ω apresentou erros relativos maiores,
102
103
7,831% e 12,572% nas malhas M1 e M2 respectivamente, pois este método só é vantajoso
em regiões de y+ < 2, e ao analisar as figuras 4.4a e 4.4b pode-se ver que poucas regiões das
duas malhas obedecem este critério.
Como esperado, os tempos de simulação na malha M2, que possui mais elementos, foi
maior. Apesar disto, os resultados obtidos para vazão mássica não foram muito melhores de
forma a compensar o grande aumento no tempo de simulação, tendo sido inclusive piores
para os modelos κ−ε e κ−ω.
Portanto, pode-se concluir que para a análise em questão o modelo SST é o mais reco-
mendável, por ter apresentado resultados melhores para a vazão mássica. Seu resultado na
malha M2 foi apenas 0,4% melhor do que na malha M1 porém o tempo de simulação foi mais
de quatro vezes maior, portanto a malha M1 é mais recomendável para o uso em conjunto
com o modelo SST por apresentarem um bom balanço entre qualidade do resultado e tempo
de processamento.
5.1 Sugestões para Trabalhos Futuros
Como foi dito na seção Rotodinâmica do capítulo de Revisão Bibliográfica, os coeficientes
rotodinâmicos não puderam ser calculados pelo método analítico descrito porém poderiam
ter sido determinados, como foi feito em [8], com a utilização do modelo bulk-flow. Portanto
fica como sugestão para trabalhos futuros o uso dos dados de distribuição de pressão e de
tensões cisalhantes calculados neste trabalho para o desenvolvimento do modelo bulk-flow
com o objetivo de se determinarem os coeficientes rotodinâmicos deste selo.
Também fica como sugestão um trabalho que analise como alterações na geometria do
selo (diâmetro dos furos, comprimento do selo, etc.) influencia nas características calculadas
neste trabalho.
6 – Apêndice A
A tabela 6.1 contém os valores da coordenada z para cada plano transversal usado para
avaliar as pressões médias.
104
105
Tab. 6.1: Configuração dos planos transversais para análise de pressão média
Nº Planos Transversais Coordenada em Z (m)
Plano 1 0.0029990
Plano 2 0.0045865
Plano 3 0.0082005
Plano 4 0.0118145
Plano 5 0.0154285
Plano 6 0.0190425
Plano 7 0.0226565
Plano 8 0.0262705
Plano 9 0.0298845
Plano 10 0.0334985
Plano 11 0.0371125
Plano 12 0.0407265
Plano 13 0.0443405
Plano 14 0.0479545
Plano 15 0.0515685
Plano 16 0.0551825
Plano 17 0.0587965
Plano 18 0.0624105
Plano 19 0.0660245
Plano 20 0.0696385
Plano 21 0.0732525
Plano 22 0.0768665
Plano 23 0.0804805
Plano 24 0.0840945
7 – Referências Bibliográficas
[1] Dara W. Childs. Annular gas seals and rotordynamics of compressors and turbines.
Texas A&M University - Proceedings of the 26TH turbomachinery symposium, 26TH,
1997.
[2] Lomakin A. Calculation of critical number of revolutions and the contidions necessary
for dynamic stability of rotors in high-pressure hydraulic machines when taking into
account forces originating in sealing. Power and Mechanical Engineering, 1(1):1, 1958.
[3] H. BLACK. Effects of hydraulic forces in annular pressure seals on the vibrations of
centrifugal pump rotors. Journal of Mechanical Engineering Science, 11(2):206–213,
1969.
[4] H BLACK e D JENSSEN. Dynamic hybrid bearing characteristics of annular controlled
leakage seals. Paper 9, 1969.
[5] D. W.. CHILDS. Dynamic analysis of turbulent annular seals based on hirs’ lubrication
equation. Journal of Tribology, 105(3):429–436, 1983.
[6] Daniela Gonçalves. Análise da dinâmica de selos de fluxo em compressores. Dissertação
de graduação, Universidade Federal do Rio de Janeiro, 2015.
[7] Eagleburgmann International. Shf seal. Disponível em www.eagleburgmann.com, 15-
07-2016.
[8] Patrick J. Migliorini. A computational fluid dynamics/bulk-flow hybrid method for
determining rotordynamic coefficients of annular gas seals. Journal of Tribology, 1(9):
1–9, 2012.
[9] Dr. Rajiv Tiwari. A Brief History and State of the Art of Rotor Dynamics. Indian Institute
of Technology Guwahati, 2008.
[10] Maurice L. Adams. Rotating Machinery Vibration from Analysis to Troubleshooting.
Marcel DekkerNew York, 2001.
106
107
[11] Mancini M. Increasing pump reliability and life. Pumps and Systems, 1(1):1, 2009.
[12] Xin Yan e Zhenping Feng Jun Li. Investigations on the rotordynamic characteristics
of a hole-pattern seal using transient cfd and periodic circular orbit model. Journal of
Vibration and Acoustics, 133(041007-1):1–9, 2011.
[13] Jim Fitch Robert Scott e Lloyd Leugner. The Practical Handbook of Machinery Lubri-
cation - Fourth Edition. Noria Corporation, 2012.
[14] Et al Daniele Massinia. Analysis of flat plate honeycomb seals aerodynamic losses: ef-
fects of clearance. Conference of the Italian Thermal Machines Engineering Association,
68(1):502–511, 2014.
[15] Giuseppe Vannini Diego Saba, Paola Forte. Review and upgrade of a bulk flow model
for the analysis of honeycomb gas seals based on new high pressure experimental data.
Journal of Mechanic Engineering, pages 1–2, 2014.
[16] Hale e Elrod Childs. Annular honeycomb seals: Test results for leakage and rotordynamic
coefficients; comparisons to labyrinth and smooth configurations. ASME Journal of
Tribology, 1(111):293–300, 1989.
[17] D.W. Childs e Kim C.H. Analysis and testing for rotordynamic coefficients of turbulent
annular seals with different,directionally homogeneous surface roughness treatment for
rotor and stator elements. ASME Journal of Tribology, 1(107):296–306, 1985.
[18] Gocha Chochua. Computations of gas annular damper seal flows. Dissertação de pos
graduação, University of Florida, 2002.
[19] Andre Luiz Tenorio Rezende. Análise numérica da bola de separação do escoamento tur-
bulento sobre a placa fina inclinada. Dissertação de doutorado, Pontifícia Uuniversidade
Católica do Rio de Janeiro - PUC-RIO, 2009.
[20] C. Nelson. Rotordynamic coefficients for compressible flow in tapered annular seals.
ASME Journal of Tribology, 1(107):318–325, 1985.
108
[21] Mihai Arghir A.Z. Rotordynamic analysis of textured annular seals with multiphase
bubbly flow. Workshop dynamic sealing under Severe working conditions, 1:3–13, 2009.
[22] Luis Fernando Azevedo. Mecânica dos Fluidos 1 Capitulo 2. Laboratório de Engenharia
de Fluidos DEM/PUC-Rio, 2016.
[23] Stephen B. Pope. Turbulent Flows - 1 ed. Cambridge University Press, 2000.
[24] Suhas V. Patankar. Numerical Heat Transfer and Fluid Flowl. CRC Press, 1980.
[25] Francieli Aparecida Vaz. Modelagem e simulação de chamas difusivas turbulentas de
etanol. Dissertação de doutorado, Universidade Federal do Rio Grande do Sul, 2013.
[26] Versteeg H K e Malalasekera W. Introduction to Computational Fluid Dynamics the
Finite Volume Method. Longman, 1995.
[27] Frank M. White. Mecânica dos Fluidos 6a Edição. Mc Graw Hill, 2011.
[28] José Francisco Almeida SOUZA. Uma revisão sobre a turbulência e sua modelagem.
Revista Brasileira de Geofísica, 29(1):21–41, 2011.
[29] EIGER S. Métodos numéricos em recursos hídricos. Coleção da ABRH-Associação
Brasileira de Recursos Hídricos, 1(2):84–155, 2011.
[30] Kundu e Cohen. Fluid Mechanics. Academic Press, 2002.
[31] Henk Tennekes e John L. Lumley. A first course in turbulence. M.I.T Press, 1:816– 819,
1973.
[32] Hinze J. O. Turbulence: Second Edition. McGraw-Hill Book Company, 1975.
[34] Reddy e D.K. Gartling. The Finite Element Method in Heat Transfer and Fluid Dynamics.
CRC Press, 1994.
[35] et al. T. Poinsot. Simulation tools for 3d reacting flows. page 5, 2013.
[36] G. Grotzbach. Direct numerical and large eddy simulation of turbulent channel flows.
Encyclopedia of Fluid Mechanics, 1987.
109
[37] Kim W. e Menon S. An unsteady incompressible navier-stokes solver for large eddy
simulation of turbulent flows. International Journal Numerical Meth. Fluids, 31(1):
983–1017, 1999.
[38] A. Silveira Neto. Simulação de grandes escalas de escoamentos turbulentos. Escola de
Primavera,Transiçã e turbulência, 1(1):159–190, 1998.
[39] Ferziger J.H. Simulation of complex turbulent flow: recent advances and prospects in
wind engineering. Journal of Wind Engineering and Industrial Aerodynamics, 46 and
47:195–212, 1993.
[40] Rogallo R.S. e Moin P. Numerical simulation of turbulent flows. Annual Review Fluid
Mechanics, 16:99–137, 1984.
[41] Lesieur M. Comte e Métais. Numerical simulations of coherent vortices in turbulence.
ASME - Application Mechanical Review, 48(4), 1995.
[42] Haxworth D. C. Large-eddy simulation of in-cylinder flow. Oil and Gas Science and
Technology- Rev. IFP, 54(2):175–185, 1999.
[43] Ansys. Ansys user manual. ANSYS, 1(1), 2010.
[44] BILESKY Marcelo Sidney Mendes, COSTA e Luciano Rossi. Método dos elementos
finitos aplicado a engenharia civil. Faculdade de Ciências Sociais e Agrárias de Itapeva,
1(1):1–11, 2014.
[45] Clovis R. Maliska. Transfência de Calor e Mecânica dos Fluidos Computacional. LTC
- Livros Técnicos e Cientícos, 2014.
[46] BALIGA e PATANKAR. A control-volume finite-element method for two-dimensional
fluid flow and heat transfer. Numerical Heat Transfer, 6(1):245–261, 1983.
[47] Ricardo Barbosa Damian. Acoplamento de balanço populacional a simulção compu-
tacional de escoamentos multifásicos polidispersos. Dissertação de pos graduação,
Universidade Federal do Rio de Janeiro, 2007.
110
[49] Fernando Sandro Velasco Hurtado. Uma formulação de volume finitos baseada em
elementos para a simulação do deslocamento bifásico imiscível em meios porosos.
Dissertação de mestrado, Universidade Federal de Santa Catarina, 2005.
[50] Henrique Krainer Eidt. Análise da influência do bocal no desempenho da camara de
expansão do separador bifásico líquido-gás tipos vasps. Dissertação de graduação,
Universidade tecnológica da federal do Paraná, 2014.
[51] Inc. XYZ Scientific Applications. Truegrid quality mesh. Disponível em
www.truegrid.com, 17-05-2016.
[52] John Chawner. Aerodynamic optimization with pointwise and friendship-framework.
Disponível em blog.pointwise.com, 18-05-2016.
[53] Inc. SAS IP. Index of software ansys. Disponível em www.sharcnet.ca, 18-05-2016.
[54] Chris Nelson. Making sense of cfd grid types. Disponível em www.innovative-cfd.com,
18-05-2016.
[55] Bruno da Costa Fonseca. Projeto de ensino de graduação (peg) "o ensino de multigrid
em álgebra linear numérica". Universidade Federal de Minas Gerais, 2009.
[56] Nigel P. Weatherill Joe F. Thompson, Bharat K. Soni. Handbook of Grid Generation.
CRC Press, 1999.
[57] J. L Wade. Test versus predictions for rotordynamic coefficients and leakage rates of
hole-pattern gas seals at two clearances in choked and unchoked condition. Master
thesis, Texas A&M University - Department of Mechanical Engineering, 2004.
[58] J Childs, D e Wade. Rotordynamic-coefficient and leakage characteristics for hole-
pattern stator annular gas seals measurements versus predictions. Journal of Tribology,
1(126):326–333, 2004.
