UNIVERSIDADE FEDERAL DE SANTA CATARINADEPARTAMENTO DE ENGENHARIA MECÂNICA

Diogo Nardelli Siebert

MODELOS CINÉTICOS DISCRETOS PARA FLUIDOSNÃO-IDEAIS COM TRANSIÇÃO DE FASES

Florianópolis(SC)

2013

Diogo Nardelli Siebert

MODELOS CINÉTICOS DISCRETOS PARA FLUIDOSNÃO-IDEAIS COM TRANSIÇÃO DE FASES

Tese submetida ao Programa de Pós-Graduação em Engenharia Mecânicapara a obtenção do Grau de Doutorem Engenharia Mecânica.Orientador: Paulo César Philippi, Dr.

Florianópolis(SC)

2013

Ficha de identi�cação da obra elaborada pelo autor,atravésdo Programa de Geração Automática da Biblioteca Universitária da UFSC

Nardelli Siebert, DiogoModelos cinéticos discretos para �uidos não-ideais com tran-

sição de fases / Diogo Nardelli Siebert ; orientador, Paulo CésarPhilippi - Florianópolis, SC, 2013.

151 p.

Tese (doutorado) - Universidade Federal de Santa Catarina,Centro Tecnológico. Programa de Pós-Graduação em EngenhariaMecânica.

Inclui referências

1. Engenharia Mecânica. 2. transição de fases. 3. �uidosnão-ideais. 4. interfaces líquido-vapor. I. Philippi, Paulo César.II. Universidade Federal de Santa Catarina. Programa de Pós-Graduação em Engenharia Mecânica. III. Título.

Diogo Nardelli Siebert

MODELOS CINÉTICOS DISCRETOS PARA FLUIDOSNÃO-IDEAIS COM TRANSIÇÃO DE FASES

Esta Tese foi julgada aprovada para a obtenção do Título de�Doutor em Engenharia Mecânica�, e aprovada em sua forma �nal peloPrograma de Pós-Graduação em Engenharia Mecânica.

Florianópolis(SC), 05 de Abril 2013.

Júlio César Passos, Dr.Coordenador

Paulo César Philippi, Dr.Orientador

Banca Examinadora:

Prof. Paulo César Philippi, Dr.Presidente

Prof. Felix Sharipov, Dr. - RelatorUniversidade Federal do Paraná

Prof. Aristeu da Silveira Neto, Dr.Universidade Federal de Uberlândia

Prof. Amir Antônio Martins de Oliveira Jr., Ph.D.Universidade Federal de Santa Catarina

Prof. Celso Peres Fernandes, Dr.Eng.Universidade Federal de Santa Catarina

Prof. Júlio César Passos, Dr.Universidade Federal de Santa Catarina

Prof. Luis Orlando Emerich dos Santos, Dr.Eng.Universidade Federal de Santa Catarina

RESUMO

De forma distinta aos métodos de CFD convencionais, os métodos LB(Lattice-Boltzmann) baseados na discretização de equações cinéticas,historicamente denominadas de equações de Boltzmann, descrevem ocomportamento dinâmico dos �uidos na escala mesoscópica. Na escalacitada, fenômenos macroscópicos complexos, como o surgimento e co-lapso de interfaces, podem ser naturalmente descritos como associadosa termos fontes incorporados às equações cinéticas. Utilizando estacaracterística do método LB, o presente trabalho propõe um novo mé-todo numérico para simulação de �uidos monocomponentes não-ideais,onde as fases líquida e vapor podem coexistir. O modelo cinético con-tínuo para este �uido tem como base a primeira equação da hierarquiaBBGKY, na qual a repulsão a curta distância é representada pelo mo-delo de esferas rígidas de Enskog e a atração a longa distância é des-crita utilizando um potencial médio. Por intermédio da análise multi-escala de Chapman-Enskog, demonstra-se que as equações de balançode massa, quantidade de movimento e energia são compatíveis comaquelas do modelo da interface difusa. Objetivando a obtenção de ummodelo mínimo, são realizadas simpli�cações nas equações cinéticas quenão comprometem o comportamento macroscópico. O estudo de umainterface líquido-vapor em equilíbrio permite a derivação de expressõespara o per�l de densidade na interface e para a tensão super�cial, assimcomo comprova que o modelo contínuo recupera corretamente as den-sidades de saturação previstas pelas leis da termodinâmica. A formadiscreta para a equação cinética é então obtida aplicando-se o métododas abscissas prescritas em conjunto com uma nova forma discreta deterceira ordem para o termo advectivo da equação de Boltzmann. A va-lidação numérica do método discreto é realizada por meio da medida davelocidade do som, da obtenção da curva de coexistência, da conferênciada lei de Laplace assim como da relação teórica para a tensão super�-cial. Tais simulações são realizadas utilizando as equações de estado deVan der Waals, Redlich-Kwong, Redlich-Kwong-Soave, Peng-Robinsone Carnahan-Starling. Os resultados demonstram que a inconsistênciatermodinâmica apresentada por métodos anteriores é eliminada com autilização do termo de advecção de terceira ordem. Por �m, o métodoé aplicado no estudo da coalescência entre duas gotas e no crescimentode bolhas em um meio líquido. Observa-se uma boa concordância comas correlações teóricas e dados experimentais disponíveis.

Palavras-chave: transição de fases. �uidos não-ideais. interfaceslíquido-vapor.

ABSTRACT

Unlike conventional CFD methods, the Lattice Boltzmann Method(LBM) describes the dynamic behavior of �uids in a mesoscopic scaleusing discrete forms of the kinetic equations. In this scale, complexmacroscopic phenomena as the formation and collapse of interfaces canbe naturally described as related to source terms incorporated intothe kinetic equations. Based on this characteristic, a novel numericalmethod for the simulation of a single component multiphase �ow is pro-posed. A continuum kinetic model for this �uid is obtained from the�rst equation of the BBGKY hierarchy using a short-range repulsionrepresented by the Enskog hard sphere model and long-range attrac-tion described by a mean-�eld potential. Using the Chapman-Enskogmultiscale analysis, it is shown that the mass, momentum and energybalance equations are compatible with those of the di�use interfacemodel. A minimal version of the model is then proposed by meansof simpli�cations which do no a�ect the macroscopic behavior. Thestudy of a �at liquid-vapor interface demonstrates that the continuummodel recovers the saturation densities predicted by the laws of ther-modynamics and expressions for the density pro�le and surface tensionare derived. A discrete form of the kinectic equation is then obtai-ned by applying the quadrature method based on prescribed abscissastogether with a new third order scheme for the discretization of theadvection term in the Boltzmann equation. The numerical validationof the method is performed by measuring the speed of sound, retrie-ving the coexistence curve and verifying both the Laplace's law and thetheoretical expression for the surface tension. The simulations are per-formed with the equations of state of Van der Waals, Redlich-Kwong,Redlich-Kwong-Soave, Peng-Robinson and Carnahan-Starling. The re-sults show that the thermodynamic consistency problem presented inearlier multiphase models are greatly reduced by using the third orderdiscretization scheme for the advection term. Simulations of the coa-lescence process of drops and of the growth of a bubble in a liquid poolare, �nally, presented with good agreement with available theoreticalpredictions and experimental data.Keywords: phase change. non-ideal �uids. liquid-vapor interface.

NOMENCLATURA

Letras Gregras

β Parâmetro Livre do modelo Lee-Fischer

γ Tensão Super�cial

Γi Massa por unidade de área armazenada na interface

ρc Massa especí�ca no ponto crítico

δ Passo de tempo

ϵ Ângulo que relaciona a posição relativa das duas partículas an-tes da colisão.

η Viscosidade Dinâmica

Θ Temperatura adimensional utilizada no LBM

κ Parâmetro de interação / Momento de segunda ordem do po-tencial

λ Segundo Coe�ciente de Viscosidade

λI Espessura da Interface

Λs Termo de geração de entropia

µ Potencial Químico

ξo Velocidade molecular adimensional

ξo,i Velocidade molecular adimensional discreta

ξ̄ Velocidade molecular pós-colisão

ρ Massa especí�ca

ρr Massa especí�ca na forma reduzida

ϱ Massa especí�ca adimensional utilizada nos métodos LB e LGA

σ Raio da esfera de ação

¯̄τ Tensor tensão viscosa

τ Tempo de Relaxação

χi Força resultante sobre a partícila i

χ Função correlação ou Função distribuição radial

Ψ Função massa efetiva do modelo Shan-Chen

Ψ Função massa efetiva do modelo Shan-Chen

ψ Energia livre de Helmholtz por unidade de volume

ψ0 Energia livre de Helmholtz por unidade de volume para umsistema homogêneo

Ω Operador de Colisão

Ωb Operador de colisão de Boltzmann

Ωi Operador de colisão discreto

ωi Peso da quadratura relacionada à velocidade ci

Ωl Termo de interação de longa distância do operador de colisão

Ωs Termo de interação de curta distância do operador de colisão

Letras Romanas

a Fator de escala

b Volume especí�co de uma partícula do gás

ci Velocidade discreta adimensional utilizada no LGA e LBM

cs Velocidade do Som

D Número de Dimensões

Dt Derivada Advectiva

¯̄E(n) Tensor isotrópico de ordem n

e Energia interna por unidade de massa

e0 Energia Interna de um sistema homogêneo

f̂ Função distribuição modi�cada de segunda ordem

f Função distribuição mássica

fi Função distribuição discreta

feqi Função distribuição de equilíbrio discreta

FN Função densidade de probabilidade de N partículas

FR Função densidade de Probabilidade reduzida

g Aceleração devido a força externa

G Parâmetro de interação do modelo Shan-Chen

H(n)rn Polinômio de Hermite de ordem n e conjunto de índices rn

h Distância entre sítios adjacentes

h Entalpia especí�ca ou entalpia por unidade de massa

¯̄I Tensor Identidade

k̂ Versor que indica a direção relativa entre duas partículas antesda colisão

k Condutividade Térmica

kx Número de Onda

ṁ Vazão mássica

m Massa da Partícula

n̂ Vetor Normal

ni Variável booleana que representa a ocupação de um sítio noLGA

p Pressão Termodinâmica

pr Pressão na forma reduzida

peq Pressão de Saturação

q Vetor �uxo de Calor

qs Vetor �uxo difusivo de entropia

R Constante do gás

r Parâmetro de impacto

S Área da região junto a interface

s Entropia por unidade de massa

¯̄T Tensor tensão

t̂ Versor tangente à superfície entre as fases

T Temperatura

Tc Temperatura no ponto crítico

To Temperatura de referência

Tr Temperatura na forma reduzida

u Velocidade adimensional utilizada nos métodos LB e LGA, co-mumente denominada de velocidade em unidades de rede

v Velocidade macroscópica do �uido

vi Velocidade da interface

V Volume da região junto a interface

V Potencial Intermolecular por unidade de massa

xi Posição da partícula i

SUMÁRIO

1 INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 DESCRIÇÃO MACROSCÓPICA DE ESCOAMEN-TOS MULTIFÁSICOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.1 EQUAÇÕES DA TERMO-HIDRODINÂMICA . . . . . . . . . . 212.2 MODELO DA INTERFACE SINGULAR . . . . . . . . . . . . . . . 232.2.1 Casos especí�cos e Simpli�cações das Condi-

ções de Contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.2.2 Limitações do modelo . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.3 MODELO DA INTERFACE DIFUSA . . . . . . . . . . . . . . . . . . 293 FUNDAMENTOS DA TEORIA CINÉTICA. . . . . . . . . 353.1 FUNÇÃO DISTRIBUIÇÃO E EQUAÇÃO DE LIOUVILLE 363.2 HIERARQUIA BBGKY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.3 EQUAÇÃO DE BOLTZMANN . . . . . . . . . . . . . . . . . . . . . . . . 383.3.1 Propriedades da Colisão e Distribuição de Equi-

líbrio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.3.2 Limite Macroscópico . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.3.3 Modelos de colisão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.4 EQUAÇÃO DE ENSKOG . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434 EQUAÇÃO LB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 474.1 MODELOS DE GÁS EM REDE . . . . . . . . . . . . . . . . . . . . . . . 474.2 ORIGEM DA EQUAÇÃO LB . . . . . . . . . . . . . . . . . . . . . . . . . 504.3 CONEXÃO COM A TEORIA CINÉTICA . . . . . . . . . . . . . . 534.3.1 Termo Advectivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.3.2 Espaço de velocidades e funcão distribuição . . . . 56

4.3.2.1 Redes e distribuições de equilíbrio . . . . . . . . . . . . . . . . . . . . 605 REVISÃO DOS MODELOS MULTIFÁSICOS . . . . . . . 635.1 MODELO SHAN-CHEN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 635.1.1 Modi�cações do Modelo Shan-Chen . . . . . . . . . . . . 67

5.2 MODELO SWIFT-ORLANDINI . . . . . . . . . . . . . . . . . . . . . . . 695.3 MODELO DE LEE & FISCHER . . . . . . . . . . . . . . . . . . . . . . 71

6 MODELO PROPOSTO . . . . . . . . . . . . . . . . . . . . . . . . . . . . 756.1 MODELO CINÉTICO COMPLETO . . . . . . . . . . . . . . . . . . . 756.2 ANÁLISE MULTI-ESCALA . . . . . . . . . . . . . . . . . . . . . . . . . . 796.2.1 Análise do Modelo Completo . . . . . . . . . . . . . . . . . . . 83

6.2.1.1 Balanço de Massa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 856.2.1.2 Balanço da Quantidade de Movimento . . . . . . . . . . . . . . . . 856.2.1.3 Balanço da Energia Interna . . . . . . . . . . . . . . . . . . . . . . . . . . 86

6.3 MODELO ISOTÉRMICO SIMPLIFICADO . . . . . . . . . . . . . 886.4 MODELO TÉRMICO SIMPLIFICADO . . . . . . . . . . . . . . . . 896.5 COMPORTAMENTO DA INTERFACE . . . . . . . . . . . . . . . . 916.6 FORMA DISCRETA DOS MODELOS SIMPLIFICADOS 966.6.1 Forma discreta do termo advectivo . . . . . . . . . . . . . 966.6.2 Forma discreta no espaço de velocidades . . . . . . . 986.6.3 Forma discreta dos operadores diferenciais . . . . . 101

7 RESULTADOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1037.1 VELOCIDADE DO SOM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1047.1.1 Solução Analítica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1047.1.2 Resultados Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . 106

7.2 CURVA DE COEXISTÊNCIA E PERFIL DE MASSAESPECÍFICA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

7.3 LEI DE LAPLACE E TENSÃO SUPERFICAL. . . . . . . . . . 1137.4 CORRENTES ESPÚRIAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1167.5 COALESCÊNCIA DE GOTAS . . . . . . . . . . . . . . . . . . . . . . . . 1217.6 CRESCIMENTO DE BOLHA EM MEIO LÍQUIDO . . . . . 123

8 CONCLUSÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125Referências Bibliográ�cas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127APÊNDICE A -- Polinômios de Hermite . . . . . . . . . . . . . . . 137APÊNDICE B -- Redes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141APÊNDICE C -- Equações de Estado . . . . . . . . . . . . . . . . . . 145APÊNDICE D -- Lista de Publicações . . . . . . . . . . . . . . . . . 149

17

1 INTRODUÇÃO

O estudo de �uidos multifásicos assim como o de escoamentosonde ocorrem transições de fase são de grande interesse tanto do pontode vista cientí�co quanto tecnológico. Entre os exemplos mais comunsde aplicações tecnológicas que envolvem esta área da mecânica dos �ui-dos estão o desenvolvimento de trocadores de calor e de sistemas derefrigeração, a recuperação de petróleo em rochas reservatório[1], o pro-jeto de usinas nucleares[2] e a criação de dispositivos micro�uidos[3]

(micro�uidics). Na área cientí�ca, escoamentos multifásicos são en-contrados em problemas que vão desde a obtenção da distribuição dotamanho de gotas da chuva na meteorologia[4] até o estudo do meca-nismo de hidratação de animais na biologia[5].

Porém, a descrição de escoamentos multifásicos encontra-se en-tre as áreas de maior complexidade da mecânica dos �uidos, tanto sobo aspecto físico quanto computacional. Esta di�culdade provém prin-cipalmente da dinâmica da interface, uma vez que esta sofre mudançasem sua topologia durante a evolução do escoamento, podendo ser criadaou destruída no decorrer de tal evolução.

O rápido avanço na capacidade de processamento e armazena-mento computacional transformou o uso da computação em uma cons-tante na área da mecânica dos �uidos. As mais conhecidas e utilizadasdentre as ferramentas de mecânica dos �uidos computacional (CDF)são aquelas que utilizam métodos numéricos para solução das equaçõesmacroscópicas de Navier-Stokes, como os métodos de volumes �nitos,elementos �nitos e diferenças �nitas. Porém, a aplicação desta abor-dagem a problemas multifásicos resulta em diversas di�culdades rela-cionadas à complexa dinâmica da interface, assim como à transição defase que ocorre nesta superfície.

Criado no início da década de 90, o método de Boltzmann pararedes ou método Lattice-Boltzmann (LB) tem se mostrado bastante pro-missor na simulação de escoamentos complexos, possuindo uma sériede vantagens em relação aos métodos CFD convencionais na descriçãode escoamentos imiscíveis[6]. Tais resultados se devem à representaçãomesoscópica utilizada em LB, sendo esta uma escala intermediária en-tre a descrição microscópica, onde a posição e velocidade de todas asmoléculas do �uido são conhecidas, e a macroscópica, a qual utiliza ahipótese do contínuo e descreve o �uido por meio dos campos de massaespecí�ca, velocidade e temperatura.

Embora com relativo sucesso na área, grande parte dos modelos

18

multifásicos existentes em LB não possuem conexão direta com a teo-ria cinética, sendo a descrição numérica a nível mesoscópico livrementeajustada com o objetivo de recuperar o comportamento macroscópicodesejado, o que resulta em equações de estado arti�ciais, na inconsistên-cia com as leis da termodinâmica e na quebra da invariância galileana.Tal abordagem foi comumente utilizada em LB, uma vez que a cone-xão detalhada entre o método e a teoria cinética foi estabelecida apenasrecentemente[7, 8]. Adicionalmente, somente na última década foi ela-borada uma descrição mesoscópica simples e consistente para �uidosnão-ideais[9] que possibilita a elaboração de métodos computacionais.

Dentro deste contexto, o presente trabalho alia o procedimentode obtenção de formas discretas das equações da teoria cinética coma descrição mesoscópica citada, propondo assim um novo modelo LBpara a descrição de �uidos não-ideais que apresentam transição de fase.A proposta origina-se de uma descrição baseada na equação de primeiraordem da hierarquia BBGKY, na qual o termo de interação é separadoem curta e longa distância, sendo estes aproximados, respectivamente,pelo operador de colisão da equação de Enskog e por um potencialmédio. Como trabalhos recentes[10] demonstram que a forma atual dométodo LB é insu�ciente para descrever tal equação cinética, uma novaforma discreta para o termo advectivo é proposta, o que em conjuntocom o método das abcissas prescritas proposta por Philippi et al.[7]

resulta em um novo método computacional.De maneira a apresentar o desenvolvimento do modelo e seus

resultados, o presente documento está organizado da seguinte forma.Inicialmente, no Capítulo 2 realiza-se uma revisão sobre escoamentosmultifásicos do ponto de vista macroscópico, abordando-se tanto o mo-delo convencional, no qual a interface é descrita como uma condiçãode contorno, quanto o modelo da interface difusa, que representa a in-terface como uma região de variação contínua da massa especí�ca. OCapítulo 3 foca a representação mesoscópica do �uido, apresentandoos fundamentos da teoria cinética , no qual as equações de Boltzmann,Enskog e a hierarquia BBGKY são expostas e discutidas. As origens,estrutura e fundamentação teórica do método LB são abordadas noCapítulo 4, dando especial ênfase a conexão deste método com a teo-ria cinética. Em sequência, o Capítulo 5 discute os principais modelosmultifásicos que fazem uso da estrutura LB, focando principalmente nospontos fortes e limitações de cada equação. A derivação do modelo me-soscópico para �uidos não-ideais é exposta no Capítulo 6. Este capítuloapresenta também a análise multiescala da equação cinética e discutesimpli�cações que não afetam o comportamento macroscópico do esco-

19

amento nas primeiras ordens do número de Knudsen. Por meio destesmodelos simpli�cados, a forma discreta é então obtida. No capítulo 7,a equação LB proposta é então validada por intermédio de simulaçõesque visam demonstrar a incorporação da equação de estado, a consis-tência termodinâmica e a validade da equação de Young-Laplace. Nestemesmo capítulo, o método é utilizado na simulação dos problemas decoalescência entre gostas e do crescimento de uma bolha em um meiolíquido a baixa pressão. Por �m, o Capítulo 8 resume as principaisconclusões do trabalho, apresentando também as possíveis extensõesdo modelo assim como indicando procedimentos que contornem as suaslimitações.

20

21

2 DESCRIÇÃO MACROSCÓPICA DE ESCOAMENTOSMULTIFÁSICOS

Como mencionado anteriormente, a descrição do escoamento comduas ou mais fases distintas é bastante complexo. Isto ocorre principal-mente em virtude da dinâmica da interface localizada entre as diversasfases, assim como da transferência de quantidade de movimento e ener-gia que ocorre por meio desta superfície. Neste capítulo, é abordadoo tratamento matemático em escala macroscópica para esta classe deescoamentos.

Primeiramente, são apresentadas, de forma sucinta, as equaçõesque descrevem a dinâmica de um �uido newtoniano composto de umaúnica fase. Com base nestas equações, são discutidas as extensõesque permitem descrever o comportamento de escoamentos multifásicosatravés dos modelos da interface singular e da interface difusa.

2.1 EQUAÇÕES DA TERMO-HIDRODINÂMICA

Sob uma ótica macroscópica, o estado de um �uido é caracteri-zado através dos campos contínuos de massa especí�ca (ρ) , velocidade(v) e temperatura (T ) , sendo desnecessário inferir diretamente sobrea estrutura molecular do mesmo. A in�uência dessa estrutura micros-cópica é observada somente de forma indireta, através dos valores doscoe�cientes de transporte como a viscosidade e a difusividade térmica.

As equações que determinam tais campos são obtidas por inter-médio da aplicação dos balanços de massa, quantidade de movimento eenergia[11, 12] em um volume material de �uido. Como resultado destebalanço, uma variação destas quantidades no interior do volume sópode ocorrer por meio de um �uxo difusivo pelas fronteiras do volumeou em decorrência de uma geração interna.

Do balanço de massa aplicado a um volume material de �uido,conforme descrito no parágrafo anterior, obtém-se a equação da conti-nuidade ou equação da conservação da massa,

∂ρ

∂t+∇ · (ρv) = 0. (2.1)

A aplicação do procedimento citado à quantidade de movimentoleva à equação

22

∂(ρv)

∂t+∇ · (ρvv) = ∇ · ¯̄T + ρg (2.2)

onde ¯̄T é o tensor tensão e g é o campo de aceleração resultante dapresença de uma força externa. O tensor tensão representa a forçapor unidade de área realizada pelos elementos de �uidos adjacentes.Do ponto de vista microscópico, este tensor está relacionado com o�uxo difusivo de quantidade de movimento proveniente do movimentodas partículas que compõem o �uido. Logo, uma expressão para estetensor em termos das variáveis ρ, v e T demanda o conhecimento daestrutura microscópica, ou, alternativamente, a utilização de relaçõesempíricas. Quando o �uido é newtoniano, o tensor �uxo de quantidadede movimento, ou tensor tensão, é a soma do tensor pressão com otensor tensão viscosa, sendo que este último é diretamente proporcionalao tensor taxa de deformação. Ao utilizar esta premissa, chega-se àexpressão

¯̄T = −p¯̄I +[η(∇v +∇Tv

)− λ(∇ · v)¯̄I

], (2.3)

onde a constante de proporcionalidade η é a viscosidade dinâmica, ¯̄Ié o tensor identidade, p é a pressão termodinâmica e λ é o segundocoe�ciente de viscosidade, o qual em geral é aproximado por λ = 2Dη.

Se é razoável assumir que a temperatura permanece constanteno problema a ser estudado, então as Eqs (2.1) e (2.2) em conjuntocom a equação de estado (EOS) do �uido formam um sistema fechadopara a determinação das variáveis ρ e v. Entretanto, se o campo detemperatura também é uma incógnita do problema, torna-se necessáriaa adição da equação de balanço de energia ao conjunto de equações,

∂(ρe)

∂t+∇ · (ρev) = −∇ · q + ¯̄T : ∇v, (2.4)

sendo e ≡ e(ρ, T ) a energia interna por unidade de massa e q o vetor�uxo de calor. De modo análogo ao tensor tensão, uma relação parao vetor q requer o conhecimento do �uxo microscópico de energia oua aplicação de uma relação constitutiva. Em geral, na ausência deradiação térmica utiliza-se para este �m a lei de Fourier, que propõeque o vetor �uxo de calor é diretamente proporcional ao gradiente detemperatura, i.e.,

q = −k∇T, (2.5)

23

na qual k é a condutividade térmica.Deste modo, as Eqs. (2.1) , (2.2) e (2.4), em conjunto com a equa-

ção de estado e com as relações constitutivas encontradas nas Eqs. (2.3)e (2.5), formam um sistema fechado para determinação dos campos demassa especí�ca, velocidade e temperatura. Como a solução analíticaexata deste conjunto de equações não é possível para a grande maio-ria dos problemas, torna-se então necessária a utilização de métodosnuméricos. Dentre os métodos numéricos utilizados na solução destasequações destacam-se os métodos de elementos �nitos, diferenças �ni-tas e volumes �nitos[11]. Enfatiza-se aqui que o método LB não consisteem uma forma discreta das equações da termo-hidrodinâmica, mas simdas equações da teoria cinética.

2.2 MODELO DA INTERFACE SINGULAR

Em um caso que envolva o escoamento de duas ou mais fases,tem-se um problema adicional relacionado à dinâmica da interface e àtransferência de massa, energia e quantidade de movimento que acon-tecem nesta região. A abordagem comumente utilizada é assumir quea interface entre as duas fases pode ser modelada por meio de umasuperfície de descontinuidade, ou seja, a região de transição não possuiespessura. Este modelo é conhecido como modelo da interface singularou modelo de Gibbs para a Interface.

Do ponto de vista microscópico, experimento e teoria demons-tram que na realidade a interface entre dois �uidos constitui uma região�na onde a massa especí�ca varia de forma contínua[2]. Porém, estazona de transição possui um comprimento típico da ordem de dezenasde Angströms[13]. Consequentemente, o comprimento da interface será106 vezes menor que o diâmetro característico de uma bolha ou umagota, o que torna razoável considerar a interface como uma superfíciede descontinuidade.

Matematicamente, o modelo da interface singular atribui domí-nios distintos a cada fase, nos quais as equações da termo-hidrodinâmicasão aplicadas utilizando os respectivos coe�cientes de transporte. A in-terface terá o papel de relacionar as variáveis nos diferentes domíniosatravés das condições nos contornos de cada domínio. Para obtençãodestas condições, o modelo emprega a hipótese que a interface está su-jeita a uma tensão denominada de tensão super�cial e denotada porγ. Esta grandeza é usualmente de�nida como a força por unidade decomprimento realizada pela interface na direção tangencial a superfície.

24

Alternativamente, o valor de γ pode ser interpretado como a energiapor unidade de área armazenada na superfície entre as fases.

Figura 1: Volume de controle junto à interface

Assim como as equações da termo-hidrodinâmica, as condiçõesde contorno na interface são obtidas através do balanço das grandezasmacroscópicas em um volume junto a região interfacial, representadona Figura 1, conjuntamente com a de�nição da tensão super�cial.

Denotando a velocidade da interface e o vetor normal à superfíciedo volume respectivamente por vi e n̂, o balanço de massa sobre ovolume V resulta na relação

d

dt

∫VρdV = −

∮Sρ(v − vi) · n̂dS, (2.6)

na qual S é a superfície do volume. A expressão anterior foi obtidaconsiderando-se que não há massa acumulada na interface1. No limiteem que a altura ϵ do volume se aproxima de zero, o lado esquerdo daexpressão anterior torna-se nulo e a Eq. (2.6) se restringe as integraisde superfície, que fornecem então a condição de contorno

ρ1(v1 − vi) · n̂ = ρ2(v2 − vi) · n̂ ≡ ṁ. (2.7)

Conforme ilustrado na Figura. 1, os índices 1 e 2 encontrados na velo-cidade e massa especí�ca indicam o lado da interface em que a variávelé avaliada. Logo, a relação anterior expressa que a vazão mássica ṁdeve ser igual em ambos os lados da interface.

Na aplicação do balanço da quantidade de movimento ao volumeda Figura. 1, nota-se que além de considerar o �uxo advectivo destaquantidade nas fronteiras do volume, é necessário adicionar a variaçãode quantidade de movimento resultante da atuação de forças sobre ovolume. Sobre esse volume de controle, identi�ca-se forças oriundas do:

1A massa por unidade de área acumulada na interface Γi depende da escolhada posição da interface o que pode ser constatado a partir do modelo da interface

difusa. Em geral de�ni-se esta posição como aquela em que Γi = 0.

25

Campo de tensão do �uido, ¯̄T · n̂;

Tensão super�cial, γt̂, onde t̂ é um vetor unitário tangente asuperfície entre as duas fases e perpendicular a linha de contatoentre o volume e a interface;

Campo conservativo ρg.

Ao utilizar as considerações expostas acima, o balanço da quantidadede movimento para o volume de interesse assume a forma

d

dt

∫VρvdV = −

∮Sρv(v − vi) · n̂dS

+

∮S

¯̄T · n̂dS +∮Cγt̂dl +

∫VρgdV (2.8)

A integral de linha contida no terceiro termo do lado direito destaequação computa a força total realizada pela interface sobre o volumede �uido, onde C indica a região de contato entre o volume e a inter-face. Este termo pode ser reescrito como uma integral de superfície,para tanto, emprega-se o Teorema de Stokes juntamente com algumasmanipulações algébricas[14], o que resulta em∮

Cγt̂dl =

∮S[∇γ − (γ∇ · n̂)n̂]dS.

Quando a expressão anterior é substituída na Eq. (2.8) e o limite deϵ → 0 é aplicado, as integrais volumétricas se anulam e os termosremanescentes fornecem a igualdade,

¯̄T2 · n̂− ¯̄T1 · n̂ = ṁ(v2 − v1)− (γ∇ · n̂)n̂+∇γ, (2.9)

que pode ser interpretada como um balanço de forças na interface.Analogamente, o balanço de energia na proximidade da interface

leva a relação,

d

dt

∫Vρ

(v2

2+ e

)dV =−

∮Sρ

(v2

2+ e

)(v − vi) · n̂dS −

∮q · n̂dS

+

∮S(v · ¯̄T ) · n̂dS +

∫Vρg · vdV, (2.10)

26

que gera uma nova condição de contorno no limite em que a altura ϵdo volume tende a zero,

ṁ

[(v222

+ e2

)−(v212

+ e1

)]= (q1 − q2) · n̂+

(v2 · ¯̄T2 − v1 · ¯̄T1

)· n̂.

(2.11)

2.2.1 Casos especí�cos e Simpli�cações das Condições de Con-torno

O uso das condições de contorno descritas nas Eqs. (2.7), (2.9)e (2.11) não é su�ciente para determinar os campos de massa especí-�ca, velocidade e temperatura, assim como a dinâmica da interface,caracterizada pela velocidade da interface vi

[15], em razão do elevadonúmero de incógnitas. Portanto, torna-se necessário realizar algumassuposições adicionais acerca do comportamento da interface.

Na situação em que o escoamento não apresenta transição de faseou quando este é composto por duas fases completamente imiscíveis, o�uxo de massa através da interface é nulo, ou seja, ṁ = 0. Ao aplicaresta premissa, a condição de contorno em (2.6), obtida do balanço demassa na interface, assume a forma

v1 · n̂ = v2 · n̂ = vi · n̂. (2.12)

Esta relação expressa a continuidade da componente normal da velo-cidade através da interface e também que a interface se moverá nadireção normal com a mesma velocidade que o �uido. Adicionalmente,aplica-se a condição de não escorregamento na interface, o que impõea continuidade da componente tangencial à interface, ou seja,

v1 − (v1 · n̂)n̂ = v2 − (v2 · n̂)n̂.

Enfatiza-se que a condição anterior não pode ser obtida somente depressupostos macroscópicos, mas é resultado da interação microscópicana interface e apoia-se em observações experimentais[15]. Com tal hi-pótese, tem-se que tanto a componente normal quanto a tangencial sãocontínuas, logo, a própria velocidade é contínua na interface.

Se a imposição de �uxo mássico nulo através da interface é em-pregada, a condição obtida do balanço da quantidade de movimentopode ser simpli�cada como

27

¯̄T2 · n̂− ¯̄T1 · n̂ = −(γ∇ · n̂)n̂+∇γ. (2.13)

A relação anterior expressa que a tensão não é contínua na interface,mas sofre um salto cujo valor é descrito pelo lado direito da equação.

O comportamento da tensão na interface pode ser interpretadode maneira mais clara quando a relação de descontinuidade é decom-posta nas direções normal e tangencial à interface. Como a tensãosuper�cial γ é de�nida somente na interface, um gradiente desta gran-deza deve ser tangente a esta superfície. Consequentemente, o produtoescalar da equação anterior com o versor normal n̂ resulta em

n̂ · ¯̄T2 · n̂− n̂ · ¯̄T1 · n̂ = −γ∇ · n̂.

O termo ∇ · n̂ é conhecido como a curvatura média da superfície quetambém pode ser calculado através da média dos inversos dos raiosprincipais de curvatura. Deste modo, conclui-se que o salto na compo-nente normal da tensão será conjuntamente proporcional a curvaturada superfície e a tensão interfacial. Especi�camente na situação emque o sistema se encontra em equilíbrio, a tensão resume-se à pres-são termodinâmica, conforme Eq. (2.3), e a relação anterior assume aforma:

p2 − p1 = γ∇ · n̂, (2.14)

também conhecida como equação de Young-Laplace.De forma semelhante, o produto escalar da Eq. (2.13) com o

versor t̂ resulta em uma relação para o salto da componente tangencial,

t̂ · ¯̄T2 · n̂− t̂ · ¯̄T1 · n̂ = ∇γ · t̂. (2.15)

Nota-se pelo lado direito da expressão anterior que uma descontinui-dade na componente tangencial da tensão só ocorrerá quando houverum gradiente da tensão super�cial. Esta força gerada por uma varia-ção da tensão super�cial ao longo da interface é denominada de efeitoMarangoni, sendo que tal variação é ocasionada pela existência de gra-dientes de temperatura ou de concentração na interface.

Por �m, a condição de contorno para a transferência de energiaatravés da interface, Eq. (2.11), pode ser simpli�cada desprezando-seos termos relativos ao trabalho devido às forças viscosas, ou seja,

ṁ(h2 − h1) = (q1 − q2) · n̂, (2.16)

onde a energia interna especí�ca foi escrita em termos da entalpia es-

28

pecí�ca h ≡ e + pρ . Se o valor das entalpias não difere sensivelmentedos seus valores na saturação, a diferença na entalpia entre as duasfases pode ser aproximada pela calor latente, L = hsat2 − hsat1 . Comtal aproximação, a Eq. (2.16) expressa que há uma diferença entre o�uxo de calor que entra e que sai da interface, causada pela geração ouabsorção de energia proveniente da transformação de fase que ocorrenesta superfície. No caso onde não há transição de fase, a condição naEq. (2.16) impõe a continuidade do �uxo de calor.

2.2.2 Limitações do modelo

Seja na versão completa ou simpli�cada, as condições de con-torno juntamente com as equações da termo-hidrodinâmica descrevemo comportamento dinâmico de um escoamento com múltiplas fases. Po-rém, esta descrição possui diversas limitações tanto do ponto de vistafísico como numérico.

Ao descrever a interface como uma descontinuidade entre doismeios, opta-se por simpli�car os processos que ali ocorrem. Esta hipó-tese é adequada para uma grande gama de situações onde o compri-mento da região interfacial é várias ordens de grandeza menor do que ocomprimento característico do problema. Entretanto, na proximidadedo ponto crítico a espessura da interface pode assumir valores da ordemde magnitude deste comprimento característico, o que torna o modeloda interface singular inválido neste caso[13].

Outra limitação nesta abordagem ocorre na descrição do pro-cesso de coalescência entre duas bolhas ou gotas. Isto ocorre pois naetapa inicial, ou seja, quando as gostas se tocam, uma singularidadesurge no ponto de encontro entre as interfaces, o que caracteriza umainde�nição da curvatura. Tem-se também o desaparecimento de umadas interfaces, já que passa a existir uma única interface após a união,o que torna o tratamento matemático bastante difícil.

Similar à situação descrita no parágrafo anterior é a do surgi-mento de uma gota ou bolha, fenômeno chamado também de nuclea-ção. Neste caso tem-se a criação de uma fase no interior de outra, oque acarreta o aparecimento de uma nova interface que por sua vezdelimita um novo domínio.

Além dos aspectos físicos do modelo, há de se considerar tambémas restrições geradas pela aplicação de métodos numéricos à solução dasequações, visto que a obtenção de resultados analíticos é possível apenasem um número restrito de casos. Quando a interface é descrita como

29

uma descontinuidade, divide-se o espaço em domínios distintos, umpara cada fase, onde a solução é acoplada através da interface. Poréma interface é livre para se mover e deformar o que gera uma condição decontorno móvel. Logo o algoritmo numérico deve ser capaz de rastreara posição da interface a cada passo de tempo, e impor as condições decontorno através de uma nova geração de malha ou de uma condição decontorno imersa. Tais abordagens geram ummétodo bastante complexoe, em muitos casos, sem base matemática consistente[13].

2.3 MODELO DA INTERFACE DIFUSA

Quando observada em escala apropriada, nota-se que a interfaceconsiste em uma região onde a massa especí�ca varia continuamenteentre as fases líquida e vapor. O modelo da interface difusa consideraesta natureza contínua, representando o escoamento líquido-vapor atra-vés de um único domínio, onde a in�uência da interface no comporta-mento dinâmico do escoamento ocorre através do surgimento de novostermos nas equações da termo-hidrodinâmica. Estas modi�cações nasrelações macroscópicas são obtidas com base na análise do sistema emequilíbrio e na imposição da consistência termodinâmica das equaçõesresultantes.

Em decorrência da característica não homogênea da interface, adescrição do estado de equilíbrio deste sistema físico requer uma exten-são dos conceitos da termodinâmica, haja vista que o sistema não podeser mais representado por um único valor da massa especí�ca. Adi-cionalmente, é possível demonstrar que uma descrição termodinâmicalimitada as variáveis locais resulta em uma tensão super�cial nula[16].Torna-se então necessário acrescentar uma dependência dos potenciaistermodinâmicos com o gradiente da massa especí�ca, que será denomi-nado de termo não local.

Dentro deste contexto, Van de Waals[16] propôs uma abordagemonde a energia livre de Helmholtz por unidade de volume, ψ, é escritacomo

ψ = ψ0(ρ, T ) +κ

2|∇ρ|2. (2.17)

O termo ψ0(ρ, T ) representa a energia livre por unidade de volume deum sistema homogêneo de massa especí�ca ρ e temperatura T . Já o se-gundo termo incorpora uma dependência da energia de Helmholtz comtermos não locais, sendo κ um parâmetro relacionado com o potencial

30

intermolecular.Com a adição desta nova variável na função ψ, o diferencial desta

função passa a ser calculado através da relação

dψ =

(∂ψ

∂T

)ρ

dT +

(∂ψ

∂ρ

)T

dρ+

[∂ψ

∂(∇ρ)

]· d(∇ρ). (2.18)

Reconhecendo a derivada no primeiro termo como a entropia por uni-dade de volume −ρs, a derivada parcial com relação à massa especí�cacomo o potencial químico µ e derivando explicitamente ψ em relaçãoao gradiente da massa especí�ca, obtém-se a expressão

dψ = −ρsdT + µdρ+ κ∇ρ · d(∇ρ). (2.19)

Com o objetivo de obter uma expressão para a variação da ener-gia interna, aplica-se a transformada de Legendre ψ = ρe − ρsT àrelação anterior, o que resulta na diferencial,

ρde = ρTds+ (ρµ− ψ) dρρ

+ κ∇ρ · d(∇ρ). (2.20)

O segundo termo do lado direito da equação anterior pode ser reformu-lado em termos da pressão termodinâmica, o que pode ser constatadoda de�nição dessa gradeza em função da energia livre de Helmholtz1

P = −∂(ψ/ρ)∂(1/ρ)

= ρ∂ψ

∂ρ− ψ = ρµ− ψ (2.21)

onde ψ/ρ é a energia livre de Helmholtz por unidade de massa e 1/ρ éo volume especí�co.

Utilizando a diferencial na Eq. (2.20), é possível obter uma rela-ção para a derivada total no tempo, também conhecida como derivadamaterial, da energia interna

ρde

dt= ρT

ds

dt+P

ρ

dρ

dt+ κ∇ρ · d(∇ρ)

dt. (2.22)

É importante observar que as Eq. (2.1), (2.2) e (2.4) continuamválidas neste caso, pois estas resultam apenas do balanço das grandezasmacroscópicas sem imposições adicionais. Consequentemente, possíveis

1Neste trabalho denotar-se-á a pressão termodinâmica de um sistema não homo-

gêneo por P , enquanto a pressão para um sistema homogêneo, que corresponde aequação de estado, será representada por p.

31

alterações no comportamento dinâmico �cam restritas à forma do ten-sor tensão e no �uxo de calor. A relação anterior pode ser utilizada parainferir a forma para q e ¯̄T que resultem em equações termodinamica-mente consistentes. Para tal �m, utiliza-se a relação para a variaçãoda entropia dada por

ρds

dt= −∇ · qs + Λs, (2.23)

que a�rma que a variação da entropia em um volume material seráresultado do balanço entre o �uxo de entropia difusiva pela fronteirado volume (qs) e a geração de entropia no interior do volume (Λs). Paraque as equações obtidas estejam em consonância com a segunda lei datermodinâmica, ou seja, que a evolução no tempo leve a um sistemacom entropia crescente, espera-se que o termo de geração de entropiaΛs seja positivo.

A substituição tanto da Eq. (2.23), que descreve a evolução notempo da entropia, quanto da equação para o balanço da energia internaem (2.22) resulta na seguinte expressão para geração de entropia,

−∇ · q − ¯̄T : ∇v = −T∇ · qs + TΛs − P∇ · v + κ∇ρ ·d(∇ρ)dt

(2.24)

O último termo da equação anterior pode ser reescrito utilizando ade�nição de derivada material,

κ∇ρ · d(∇ρ)dt

= ∇ ·[κ∇ρdρ

dt

]− dρdt

∇ · (κ∇ρ)− κ∇ρ∇ρ : ∇v

(2.25)

Isolando o termo de geração na expressão (2.24) e substituindoa Eq. (2.25) para o termo não-homogêneo obtém-se a expressão �nalpara Λs

Λs =1

T

[¯̄T +

(P − κρ∇2ρ

) ¯̄I + κ∇ρ∇ρ] : ∇v+∇ ·

[qs −

1

T

(q + κ∇ρdρ

dt

)]− 1T 2

(q + κ∇ρdρ

dt

)·∇T.

(2.26)

Uma das maneiras de restringir que o termo de geração seja sempre

32

positivo é escolhendo um �uxo de calor e um tensor tensão[13, 17, 18] que,ao ser substituído na equação anterior, resultem na mesma expressãode Λs encontrada para um �uido Newtoniano que obedece a Lei deFourier. O uso desta abordagem resulta nas leis constitutivas

¯̄T = (−P + κρ∇2ρ)¯̄I − κ∇ρ∇ρ+ ¯̄τ , (2.27)

q = −κ∇ρdρdt

− k∇T. (2.28)

onde ¯̄τ é o tensor de tensão viscosa, cuja forma é encontrada no segundotermo da Eq. (2.3).

Vale notar que, apesar de existirem outras escolhas de ¯̄T e qque estejam de acordo com a segunda lei da termodinâmica, a rela-ção acima constitui a alternativa mais simples que resulta nas relaçõesconstitutivas tradicionais longe da interface, ou seja, quando κ∇ρ podeser desprezado. A dedução aqui realizada faz uso de uma abordagempuramente macroscópica, porém a forma apresentada para o �uxo decalor e tensor tensão estão em consonância com resultados obtidos apartir de pressupostos microscópicos. Utilizando como ponto de par-tida a equação de Lioville, é possível inferir que em �uidos constituídosde partículas que interagem a longa distância deverão existir termosadicionais no �uxo de calor e no tensor tensão[19].

Da análise da Eq. (2.27), nota-se a presença de novos termosem relação à expressão para o tensor tensão encontrada na equação deNavier-Stokes. Em especial, tem-se a adição do termo não isotrópico

¯̄S = κ∇ρ∇ρ.

Este tensor foi originalmente postulado por Korteweg [Korteweg[20],1901 apud Anderson, McFadden e Wheeler[18], 1998, p. 142] e por estemotivo recebeu seu nome, sendo também conhecido como tensor de ca-pilaridade. Ressalta-se também que a pressão P não se resume a pressãode um sistema homogêneo com uma massa especí�ca ρ. Isto pode serconstado substituindo a expressão para a energia livre de Helmhotz porunidade de volume, Eq. (2.19), na de�nição da pressão encontrada naEq. (2.21),

P = ρ∂ψ0

∂ρ− ψ0 − κ

2(∇ρ)2 = p− κ

2(∇ρ)2,

na qual p a é pressão para um sistema homogêneo. Com esta consta-tação a tensão assume a forma �nal

33

¯̄T =[−p+ κ

2(∇ρ)2 + κρ∇2ρ

]¯̄I − κ∇ρ∇ρ+ ¯̄τ . (2.29)

Conforme constata-se na Eq. (2.28), o �uxo de calor contém umtermo adicional quando comparado a relação clássica da Lei de Fou-rier. Do ponto de vista microscópico, a interpretação física deste termonão corresponde a de um �uxo de energia têrmica mas sim a de umatransferência de energia devido ao trabalho realizado pela força a longadistância[Dunn[21], 1986 apud Anderson, McFadden eWheeler[18], 1998,pg. 146].

Uma compreensão mais clara acerca da in�uência dos termosnão clássicos na dinâmica da energia interna e, consequentemente, datemperatura pode ser obtida quando a energia interna é reescrita como

ρe = ψ + ρsT = ψ0 + ρs0T +κ

2(∇ρ)2 = ρe0 +

κ

2(∇ρ)2

onde e0 representa a energia interna por unidade de massa de um sis-tema homogêneo, ou seja, onde não é considerada a dependência comgradientes. Através da relação acima é possível obter que a derivadatotal no tempo da energia interna é dada por

ρde

dt= ρ

de0dt

− κ2ρ

(∇ρ)2 dρdt

+ κ(∇ρ) · d(∇ρ)dt

= ρde0dt

− κ2(∇ρ)2(∇ · v) +∇ ·

[κ∇ρdρ

dt

]− dρdt

∇ · (κ∇ρ)− κ∇ρ∇ρ : ∇v (2.30)

sendo aplicada a relação obtida na Eq. (2.25). Utilizando esta relaçãona Eq. (2.10) que expressa o balanço de energia em um volume material,juntamente com as expressões para o tensor tensão e �uxo de calor,contidas respectivamente nas Eq. (2.29) e (2.28), resulta na igualdade

ρde0dt

= ∇ · (k∇T ) + (−p¯̄I + ¯̄τ ) : ∇v, (2.31)

Observa-se que a evolução no tempo da variável e0, que é a energiainterna sem considerar termos não-homogêneos, independe dos termosnão clássicos. Na realidade a relação anterior é exatamente a forma

34

da equação clássica para a energia interna. Deste modo é possívelconcluir que os termos não clássicos não in�uenciam na evolução docampo de temperatura, o que por sua vez leva a interpretação de queo termo adicional no �uxo de calor, Eq (2.28), não representa um �uxodifusivo de energia térmica, mas consiste em um trabalho cuja atuaçãoin�uencia somente os termos não homogêneos da energia interna.

A utilização do tensor tensão da Eq. (2.29) juntamente com arelação anterior para a energia possibilita descrever um sistema líquido-vapor[22] sem a necessidade da imposição de condições de contorno emfronteiras móveis ou da separação das fases em domínios distintos, oque torna seu uso bastante vantajoso em métodos numéricos.

35

3 FUNDAMENTOS DA TEORIA CINÉTICA

Na descrição de um sistema físico, seja em problemas puramentecientí�cos ou de engenharia, almeja-se obter uma representação ma-temática para o sistema em questão, onde as grandezas do problemapossam ser representadas quantitativamente. Em geral, tal represen-tação pode ser realizada de diferentes formas, dependendo do nível dedetalhe necessário. Para o problema da dinâmica de um �uido, estadescrição pode ser realizada em três níveis ou escalas distintas: a mi-croscópica, a mesoscópica e a macroscópica.

Do ponto de vista microscópico, um �uido é composto por umgrande número de partículas que interagem entre si. A descrição ma-temática nesta escala é realizada através da posição e velocidade decada uma das partículas, sendo estas variáveis governadas pelas leisda mecânica. Quando analisada nesta escala, tanto a dinâmica dainterface quanto o processo de transição de fase são fenômenos queocorrem naturalmente como resultado das interações intermoleculares,o que permite observar tais fenômenos com grande riqueza de detalhe.Em contrapartida, a aplicação prática desta abordagem é bastante li-mitada, em virtude do elevado número de graus de liberdade, o quetorna proibitiva a simulação computacional mesmo de problemas comcomprimentos característicos da ordem do micrometro.

Conforme abordado no Capítulo 2, a representação macroscópicautiliza os campos de massa especí�ca, velocidade e temperatura paracaracterizar quantitativamente o �uido. Sob a ótica macroscópica, o�uido é tratado como um meio contínuo, sendo desnecessário inferir so-bre sua estrutura molecular. Por este motivo, o número de variáveis émenor quando comparado com a abordagem microscópica. No entanto,parte da dinâmica do �uido é consequência direta do comportamentomicroscópico e não pode ser obtida somente de pressupostos macroscó-picos, o que torna necessária a utilização de equações constitutivas parao cálculo dos �uxos difusivos, assim como o uso de hipóteses quanto aocomportamento da interface.

Entre a escala macroscópica, que representa o �uido como sendocontínuo, e a escala microscópica, que segue individualmente cada par-tícula, existe uma escala que pode ser considerada como intermediária.Essa escala, denominada mesoscópica, ainda considera o �uido comoconstituído por partículas, porém, ao invés de acompanhar a trajetó-ria de cada partícula, leva em conta o comportamento microscópicoapenas de forma probabilística. Consequentemente, esta abordagem se

36

bene�cia de uma representação matemática reduzida, enquanto retém acaracterística de representar fenômenos complexos como consequêncianatural das interações intermoleculares.

Tendo em vista tais características da escala mesoscópica, estetrabalho utilizou esta como base para a estruturação do modelo paraescoamentos multifásicos. Assim, neste capítulo são apresentados osconceitos e equações básicas utilizadas na teoria cinética, já que estacontém o arcabouço teórico para descrição de �uidos nesta escala.

3.1 FUNÇÃO DISTRIBUIÇÃO E EQUAÇÃO DE LIOUVILLE

Na análise microscópica, um �uido pode ser considerado comosendo um sistema de N partículas que interagem entre si. A caracte-rização deste sistema é realizada através das posições x1,x2, . . . ,xN evelocidades ξ1, ξ2, ξ3, . . . , ξN das N partículas, que devem ser conheci-das para todos os instantes de tempo.

Na teoria cinética, utiliza-se ainda a premissa de que o �uido écomposto por partículas, porém opta-se por uma abordagem probabi-lística em detrimento da descrição das posições e velocidades em funçãodo tempo. Em um contexto geral, a função que caracteriza o estadodo sistema é a função densidade de probabilidade de N partículas FN ,onde a quantidade

FN (x1, . . . ,xN , ξ1, . . . , ξN , t)dx1 . . . dxNdξ1 . . . dξN

representa a probabilidade do sistema se encontrar na região do es-paço de fase com volume dx1 . . . dxNdξ1 . . . dξN localizada no pontox1, . . . ,xN , ξ1, . . . , ξN no tempo t.

Pode-se mostrar a partir das leis da mecânica que a evoluçãodessa função ocorre segundo a equação de Liouville[23],

∂FN∂t

+N∑i=1

[∂FN∂xi

· ξim

+∂FN∂ξi

· χi]= 0. (3.1)

onde χi é a força resultante que atua sobre a i-ésima partícula e m éa massa da partícula[24]. Interpreta-se a equação de Liouville como aa�rmação de que a probabilidade é constante ao longo de uma trajetóriado sistema no espaço de fase.

A Eq. (3.1) fornece a base para a descrição probabilística docomportamento das partículas que compõem o �uido. Porém, assim

37

como ocorre na abordagem puramente microscópica, não é possível ob-ter soluções desta equação para problemas realistas, uma vez que afunção FN possui 6N +1 variáveis independentes. Contudo, é possíveldiminuir o número de variáveis por meio de uma série de simpli�ca-ções e modi�cações que reduzem o grau de detalhe com que o sistemaé descrito. Tais simpli�cações da equação de Liouville resultam nasdiferentes equações da teoria cinética, dentre as quais destacam-se aequação de Boltzamnn, a hierarquia BBGKY e a equação de Enskog.

3.2 HIERARQUIA BBGKY

Uma das formas de reduzir o grau de detalhe na descrição dosistema é diminuindo o número de variáveis na função distribuição.Com esta �nalidade, de�ne-se a função distribuição reduzida1

FR(x1, . . . ,xR, ξ1, . . . , ξR, t) =

∫FNdxR+1 . . . dxNdξR+1 . . . dξN .

(3.2)O papel da integração é contrair a descrição de modo que FR forneçaa densidade de probabilidade de encontrar as R primeiras partículasem certas posições e velocidades, independente da posição e velocidadeque se encontram as partículas de R + 1 até N . No caso especí�co deR = 1,

F1(x1, ξ1, t) =

∫FNdx2 . . . dxNdξ2 . . . dξN .

Esta função é denominada de função distribuição para uma partículae representa a densidade de probabilidade de encontrar a partícula 1na posição x1 com velocidade ξ1 no tempo t independente da posiçãoe velocidade das outras partículas.

De modo a obter uma equação que determine a evolução dafunção distribuição reduzida para R partículas, integra-se a equaçãode Liouville, Eq. (3.1), sobre as posições e velocidades das (N − R)últimas partículas. Aplicando o procedimento descrito no caso R = 1e assumindo que as partículas são idênticas obtém-se a equação paraevolução da função distribuição de uma partícula,

1Quando não especi�cado os limites de integração, convenciona-se que a integral

é realizada sobre todo o domínio das variáveis de integração

38

∂F1∂t

+∂F1∂x1

· ξ1 +∂F1∂ξ1

· χ1em

= − (N − 1)m

∫χ12 ·

∂F2∂ξ1

dx2dξ2. (3.3)

As forças χ1e e χ12 indicam respectivamente, a força externa realizadasobre a partícula 1 e a força que a partícula 2 faz sobre a partícula 1

Nota-se que a equação anterior depende da função distribuiçãode probabilidade para duas partículas F2(x1,x2, ξ1, ξ2, t). Assim comoa equação para F1 depende de F2, de maneira geral a relação obtidaatravés deste procedimento para FR será dependente de FR+1. Estadependência gera um conjunto de equações, conhecido como hierarquiade equações BBGKY[23].

Do ponto vista prático, a hierarquia BBGKY não representa umganho real. Isto ocorre pois troca-se a descrição através de uma funçãode 6N+1 variáveis da equação de Liouville por um sistema de N equa-ções diferenciais. Entretanto, empregando simpli�cações que permitamescrever a função de probabilidade conjunta F2 em termos da funçãoF1, é possível obter uma equação fechada para a função distribuiçãode uma partícula. Esta abordagem será utilizada posteriormente naderivação do modelo de transição de fase.

3.3 EQUAÇÃO DE BOLTZMANN

A primeira equação de evolução fechada para a função distri-buição F1(x, ξ, t) foi proposta por Ludwig Boltzmann no século XIX.Ainda que esta equação possa ser derivada por outro meio, a deduçãoa partir da equação de Liouville permite analisar em maior detalhe oslimites de aplicabilidade da equação de Boltzmann. Por se tratar deum procedimento onde as hipóteses e simpli�cações são evidentes, estaabordagem permite também observar os pontos que podem ser alte-rados na dedução de modo a expandir o limite de aplicabilidade daequação resultante.

Na dedução da equação de Boltzmann, as seguintes hipóteses sãoutilizadas:

a massa especí�ca do gás é su�cientemente baixa, de modo quecolisões envolvendo três partículas ou mais são altamente impro-váveis. Logo, somente colisões binárias são consideradas. Poroutro lado, é necessário que o número de partículas N contidasno gás seja su�cientemente grande para que a descrição probabi-lística seja válida;

39

o potencial interação é de curto alcance. Pois mesmo em gasesrarefeitos o conceito de colisão binária perde o sentido se o poten-cial tiver um alcance longo o su�ciente, já que uma dada partículainteragiria com muitas outras em um dado instante;

a hipótese do caos molecular é válida. Esta hipótese consiste ema�rmar que antes de duas partículas colidirem, suas velocidadese posições não estão correlacionadas.

Tais hipóteses permitem considerar que a interação entre duas partí-culas está restrita a uma região de raio σ em torno da partícula, que édenominada de esfera de ação. Fora da esfera não há interação e, conse-quentemente, a força entre partículas χ12 é nula quando |x1−x2| > σ.

Como a esfera de ação possui um volume desprezível quandocomparada ao volume total do gás, é razoável a�rmar que a probabi-lidade de encontrar uma partícula dentro da esfera de ação de outra énegligenciável. Assim, ao realizar a contração da função distribuiçãode N partículas pode-se desprezar a região do espaço de fase onde aspartículas estão separadas de uma distância inferior a σ, ou seja,

F1(x1, ξ1, t) =

∫|x2−x1|≥σ

F2(x1,x2, ξ1, ξ2, t)dx2dξ2.

É possível demonstrar que a utilização da aproximação anterior noprocesso de contração da Equação de Liouville faz com que o valor dafunção distribuição conjunta F2 seja necessário somente na superfícieda esfera de ação.

A dependência da função F2(x1,x2, ξ1, ξ2, t) pode ser completa-mente eliminada da equação por meio da hipótese do caos molecular.Na situação em que as velocidades das duas partículas ξ1 e ξ2 indicamuma situação pré-colisional, a hipótese a�rma que as probabilidadesnão são correlacionadas, logo

F2(x1,x2, ξ1, ξ2, t) = F1(x1, ξ1, t)F1(x2, ξ2, t).

Já na situação onde ξ1, ξ2, x1 e x2 especi�cam um estado pós-colisional,a não correlação das velocidades não pode ser utilizada. Para contornaresta situação, usa-se o princípio de que a probabilidade é constante emuma trajetória no espaço de fase. Através das leis da mecânica é possívelobter as velocidades que as partículas devem possuir antes da colisão,denotadas de ξ̄1 e ξ̄2, para que as velocidades �nais sejam ξ1 e ξ2. Estaconstatação permite expressar a distribuição conjunta na forma

40

F2(x1,x2, ξ1, ξ2, t) = F2(x1,x2, ξ̄1, ξ̄2, t) = F1(x1, ξ̄1, t)F1(x2, ξ̄2, t).

Onde a função F2(x1,x2, ξ̄1, ξ̄2, t) foi escrita em termos da distribuiçãode probabilidade para uma partícula F1, já que esta é calculada em umestado pré-colisional e, por hipótese, este não é correlacionado.

Com as hipóteses e simpli�cações expostas, a redução da equaçãode Liouville a uma única partícula resulta na celebrada equação deBoltzmann,

∂f

∂t+ ξ1 ·

∂f

∂x1=

1

m

∫|ξ1 − ξ|

[f(x1, ξ̄1, t)f(x1, ξ̄2, t)

− f(x1, ξ1, t)f(x1, ξ2, t)] rdrdϵdξ2 (3.4)

onde a função distribuição f é de�nida como,

f(x1, ξ1, t) = (Nm)F1(x1, ξ1, t).

Se a hipótese ergódica é considerada verdadeira, esta nova função dis-tribuição pode ser interpretada como a massa de partículas por unidadede volume do espaço de fase que se encontram na posição x1 com veloci-dade ξ1. Valendo-se desta interpretação, é possível obter uma conexãosimples com as variáveis macroscópicas por meio dos momentos destadistribuição,

ρ =

∫f(x1, ξ1, t)dξ1, (3.5)

ρv =

∫f(x1, ξ1, t)ξ1dξ1, (3.6)

D

2ρRT =

1

2

∫f(x1, ξ1, t)|ξ1 − v|2dξ1. (3.7)

na qual D é o número de dimensões e R é a constante do gás.O lado direito da Eq. (3.4) é comumente denominado de operador

de colisão, já que é este termo o responsável por descrever a ação dainteração molecular. As variáveis de integração r e ϵ correspondemrespectivamente ao parâmetro de impacto e ao ângulo de entrada dacolisão, conforme ilustrado na Figura 2. De modo simpli�cado, pode-seargumentar que esta integral percorre todas as possíveis con�guraçõesde uma colisão binária, calculando dessa maneira a taxa de variação

41

Figura 2: Ilustração da esfera de ação com as variáveis de integraçãoutilizadas no operador de colisão.

da probabilidade provocada pelas colisões com as outras moléculas queconstituem o gás.

A equação integro-diferencial obtida por Boltzmann é de difícilsolução, em razão principalmente da forma complexa do operador decolisão. Consequentemente, a utilização desta equação �ca restrita aocálculo teórico dos coe�cientes de transportes, como a viscosidade e adifusividade térmica, à derivação de equações macroscópicas e à soluçãode escoamentos em geometria simples, restrito em sua maioria a casosunidimensionais.

3.3.1 Propriedades da Colisão e Distribuição de Equilíbrio

O operador de colisão, cuja forma é dada pelo termo no ladodireito da Eq. (3.4), tem o papel de modi�car a distribuição f(x1, ξ1, t)levando em consideração a interação entre as partículas. Se tal in-teração consiste em colisões elásticas, a atuação deste operador nãodeve alterar a massa, a quantidade de movimento e a energia locais,limitando-se a afetar a forma como estas quantidades são transferidas.Denotando o operador como Ω[f ], estas a�rmações podem ser repre-sentadas na forma de equações, respectivamente, como∫

Ω[f ]dξ1 = 0, (3.8)∫Ω[f ]ξ1dξ1 = 0, (3.9)

42

∫Ω[f ]

ξ212dξ1 = 0. (3.10)

De especial interesse é a forma que a função distribuição assumeno equilíbrio, ou seja, quando Ω[f ] = 0. No caso do operador de colisãoda equação de Boltzmann, é possível demonstrar que esta relação ésatisfeita pela função

feq(ξ) = ρ

(1

2πRT

)D/2e−

(ξ−u)22RT (3.11)

a qual é denominada de função distribuição de Maxwell-Boltzmann ouMaxwelliana.

3.3.2 Limite Macroscópico

Um dos aspectos mais importantes de um modelo cinético é oseu comportamento em nível macroscópico. Como comportamento emnível macroscópico entende-se a dinâmica das variáveis ρ, u e T nolimite onde o comprimento do sistema é algumas ordens de grandezamaior que a distância média entre as partículas que constituem o gás.

A técnica mais utilizada para realizar a conexão entre as duasescalas é conhecida como análise de Chapman-Enskog[25]. Esta aná-lise consiste em uma expansão da função distribuição e dos gradientes,sendo que esta é realizada para um estado próximo do equilíbrio local.A aplicação deste procedimento à equação de Boltzmann resulta nasequações de Navier-Stokes-Fourier, ou seja, até uma dada aproximaçãoo gás se comporta como um �uido newtoniano que obedece a lei deFourier.

Este procedimento será descrito em maior detalhe no Capítulo6, onde a análise multiescala do modelo proposto é realizada.

3.3.3 Modelos de colisão

Com o intuito de simpli�car o operador de colisão, e consequen-temente a própria equação de Boltzmann, diversos autores propuseramversões alternativas para este termo. Dentre os modelos propostos naliteratura, o operador BGK encontra-se entre os de forma mais simples.Este operador foi proposto em 1954 por Bhatnagar, Gross e Krook[26]

e apresenta a forma

43

ΩBGK =1

τ(f − feq) . (3.12)

Nota-se a partir da expressão anterior que quanto mais distante o �uidoestiver do equilíbrio local maior será a taxa de variação da função distri-buição devido às colisões. O termo 1/τ é interpretado como a frequênciaentre colisões, sendo τ um parâmetro ajustável denominado tempo derelaxação.

O operador da Eq. (3.12) mantém as principais característicasdo operador de colisão originalmente proposto por Boltzmann. En-tre estas cita-se a conservação local da massa, da quantidade de movi-mento e da energia, assim como um equilíbrio descrito pela distribuiçãode Maxwell-Boltzmann. Como consequência destas características, asequações macroscópicas resultantes deste modelo também possuem aforma das equações de Navier-Stokes-Fourier, sendo os coe�cientes detransportes linearmente dependentes do tempo de relaxação. Comoprincipal limitação deste modelo, tem-se que tanto a viscosidade ηquanto a condutividade térmica k são proporcionais ao tempo de rela-xação, o que impossibilita o ajuste independente destes coe�cientes detransporte.

3.4 EQUAÇÃO DE ENSKOG

A aplicabilidade da equação de Boltzmann é limitada pelo con-junto de suposições realizadas ao derivá-la, ou seja, seu escopo de apli-cação �ca restrito a gases de baixa densidade, onde as hipóteses citadasanteriormente são válidas. A possibilidade de descrever sistemas densosutilizando uma versão generalizada da equação de Boltzmann foi alvode grande interesse nas últimas décadas. Porém, apesar desses esforços,nenhuma solução geral foi ainda encontrada[23].

Por este motivo, a equação proposta por Enskog em 1922 aindaé amplamente utilizada para descrição de gases densos. Enskog propôsmodi�car a equação de Boltzmann de modo a descrever um gás com-posto por esferas rígidas. Consequentemente, algumas das hipótesesrealizadas anteriormente por Boltzmann precisaram ser abandonadas.

A primeira delas é o conceito de colisão localizada. Nota-se naEq. (3.4) que todas as funções distribuições são avaliadas no ponto x1.Esta imposição é uma aproximação, onde as posições x2, x̄1 e x̄2 soba esfera de ação são substituídas pela posição x1, pois no limite emque raio de ação σ tende a zero todas as posições convergem para x1.

44

Esta suposição é aceitável para partículas pontuais, porém deixa de serválida quando o gás é composto de esferas rígidas.

A segunda hipótese que precisa ser descartada é a do caos mo-lecular. Em um gás denso composto por esferas rígidas, independen-temente da situação ser pré ou pós-colisional, nem todo o volume dogás está acessível em decorrência da ocupação do espaço pelas outraspartículas do gás. Logo, observa-se que haverá uma correlação entreas posições das partículas. Outro ponto que pode afetar a correlaçãoentre as posições das partículas é o encobrimento. Este efeito ocorrequando duas partículas se encontram próximas de tal modo que estasservem de escudo parcial contra a colisão de uma terceira partícula,diminuindo assim a probabilidade de ocorrência de uma colisão.

Com o relaxamento das hipóteses empregadas anteriormente naderivação da equação de Boltzmann, torna-se necessária a incorporaçãodo termo de correlação na distribuição de probabilidade conjunta. Comesta adição, obtém-se para a situação pré-colisional

F (2)(x1,x2, ξ1, ξ2) = χ

(x1 −

1

2σk̂

)F (1)(x1, ξ1)F

(1)(x1 − σk̂, ξ2).

Como ilustrado na Figura 3, k̂ é um versor que aponta do centro dapartícula 2 para o centro da partícula 1 e σ é o diâmetro da partícula,que neste caso coincide com o raio da esfera de ação. A função corre-lação χ, também denominada de função distribuição radial, é avaliadano ponto de encontro, ou seja, no ponto médio entre as duas partículas.Uma primeira aproximação para a função correlação pode ser obtidado aumento da probabilidade de se encontrar a partícula em um dadoponto causada pela ocupação do espaço vizinho por outras partículas

χ =1

1− ρ(4πσ3)/(3m)=

1

1− bρ,

onde b ≡ (4πσ3)/(3m), é o volume especí�co de uma molécula. Umfator de correlação que leve em consideração mais efeitos de interaçãoentre as partículas pode ser obtido utilizando as técnicas da mecânicaestatística[27].

Realizando os mesmos procedimentos que levam à equação deBoltzmann conjuntamente com as modi�cações descritas acima, obtém-se uma nova forma para o operador de colisão,

45

Figura 3: Colisão de duas eferas rígidas conforme o modelo de Enskog

Ω[f ] =1

m

∫|ξ1 − ξ|

[χ(x1 +

σ

2k̂)f(x1, ξ̄1, t)f(x1 + σk̂, ξ̄2, t)

−χ(x1 −

σ

2k̂)f(x1, ξ1, t)f(x1 − σk̂, ξ2, t)

]rdrdϵdξ2,

sendo a equação que utiliza este operador denominada de Equação deEnskog. A exemplo do que ocorre com o operador de colisão obtido porBoltzmann, é necessária a simpli�cação da expressão anterior quandoo objetivo é a aplicação a problemas realistas em mecânica dos �uidos.Isto pode ser alcançado através da expansão das funções f e χ em sériede Taylor ao redor do ponto x1, em conjunto com a suposição de quea função distribuição descreve um estado próximo do equilíbrio local.Estas hipóteses resultam no operador de colisão simpli�cado de Enskog,

Ω[f ] = Ωb[f ]− bρχf0{(ξ − u) ·

[∇ ln(ρ2χT )

+3

5

(|ξ − u|2

RT− D + 2

2

)∇ lnT

]+

2

D + 2

[2

RT(ξ − u)(ξ − u) : ∇u+

(|ξ − u|2

RT− D + 2

2

)∇ · u

]}(3.13)

na qual Ωb representa o operador de colisão obtido por Boltzmann. Arelação anterior constitui uma escolha mais adequada para representar

46

o efeito das interações de curta distância, visto que as hipóteses con-sideradas são menos restritivas que aquelas utilizadas por Boltzmann,tornando esta mais adequada àbi descrição de gases densos.

47

4 EQUAÇÃO LB

Neste capítulo, será realizada uma discussão dos conceitos eda estrutura básica do método de Boltzmann para redes ou métodoLattice Bolzmann (LBM). É realizada, inicialmente, uma revisão dasorigens históricas do método, começando pela sua primeira derivaçãopor meio dos modelos de gás em rede ou Lattice Gas Automata (LGA),passando então para os modelos com operador de colisão simpli�cado.Na segunda parte do capítulo são apresentadas as formulações maisrecentes sobre o método, onde este é interpretado como uma formadiscreta das equações da teoria cinética dos gases.

4.1 MODELOS DE GÁS EM REDE

Historicamente, o LBM surgiu através do aprimoramento de ou-tro método numérico, o LGA ou, como também é conhecido, métodobooleano. Este por sua vez tem origem no artigo publicado em 1986por Frisch, Hasslacher e Pomeau[28] que propôs um método baseado emautômatos celulares capaz de reproduzir o comportamento macroscó-pico de �uidos newtonianos.

A estrutura do modelo de gás em rede consiste em um mundomicroscópico simpli�cado, onde as partículas que constituem o �uidoestão limitadas a ocupar os vértices de uma rede regular, também co-nhecidos como sítios. A velocidade das partículas, denotada por ci, érestrita a um conjunto de b velocidades, sendo que o módulo e a direçãodestas são escolhidos de maneira que as partículas transladem até ossítios mais próximos em um passo de tempo.

Se a ocupação de um determinado sítio é limitada a uma partí-cula por direção, o estado de cada sítio pode ser descrito através de umconjunto de variáveis booleanas ni(x, t), a qual assume o valor 1 paraa presença e 0 para a ausência de uma partícula com velocidade ci naposição x no tempo t. O estado de um sítio para um gás em rede queutiliza um conjunto de velocidades de forma hexagonal é apresentadona Figura 4, juntamente com sua representação matemática.

A evolução do sistema é caracterizada por duas etapas distintas:a colisão e a propagação, como ilustrado na Figura 5. Na propagação,as partículas movem-se até os sítios vizinhos na direção de sua veloci-dade ci. Já na etapa de colisão, o estado sítio de cada sítio é alterado de

48

1

2 3

4

56

Figura 4: Modelo do gás em rede e representação da ocupação atravésde variáveis boolenanas.

acordo com uma regra que conserva o número de partículas e a quan-tidade de movimento no sítio. Esta etapa tem como intuito simular acolisão das partículas neste mundo microscópico simpli�cado.

Colisão Propagação

Figura 5: Etapas de colisão e propagação no método booleano emuma rede hexagonal.

No LGA utiliza-se, por simplicidade, a premissa de que as par-tículas tem massa unitária. Assim, a massa especí�ca1 pode ser iden-ti�cada como o número de partículas no sítio,

ϱ =

b∑i=0

ni(x, t), (4.1)

1Neste trabalho adotou-se duas notações diferentes para as grandezas macros-

cópicas. Os símbolos ϱ , u e Θ são utilizados, respectivamente, para representar amassa especí�ca, velocidade, e temperatura nas versões adimensionais empregadas

nos métodos LGA e LBM, comumente denominadas de variáveis em unidade de

rede. Os símbolos ρ, v e T são utilizados para representar estas mesmas grandezasna forma dimensional.

49

e a velocidade pode ser obtidade através da quantidade de movimentototal no sítio,

ϱu =b∑

i=0

nici. (4.2)

A utilização das relações anteriores e a repetida sucessão dos processosde colisão e propagação possibilitam determinar a evolução do campode massa especí�ca e velocidade do �uido.

A primeira di�culdade encontrada na utilização do LGA foi aobtenção de um comportamento macroscópico isotrópico. Este pontonegativo pode ser contornado através da escolha de um conjunto develocidades ci cujos tensores calculados pela soma

∑i cici . . . ci são

isotrópicos até quarta ordem. Quando esta condição de isotropia é sa-tisfeita, observa-se por meio de uma análise multiescala que a dinâmicado �uido segue em média a equação de Navier-Stokes para escoamentosincompressíveis.

O método LGA demonstrou ser uma ferramenta útil no estudode problemas em mecânica dos �uidos, por reproduzir uma dinâmicamacroscópica não linear e extremamente complexa utilizando um al-goritmo simples, de fácil implementação e baixo custo computacional.Adicionalmente, a interação do �uido com fronteiras sólidas pode serfacilmente modelada utilizando a re�exão das partículas que vão ao en-contro do contorno. Esta condição de contorno, conhecida como bounceback, garante velocidade nula no contorno e permite uma fácil aplicaçãodo LGA a geometrias complexas.

Apesar dos pontos positivos relacionados acima, este métodoapresentou várias de�ciências signi�cativas. Dentre estas ressalta-se:a) a dependência da pressão com a velocidade do �uido, problema quepode ser atenuado permitindo que mais partículas possam estar emrepouso nos sítios, b) a ausência de invariância galileana na equaçãomacroscópica para a evolução da quantidade de movimento, c) a exis-tência de um limite inferior para o valor da viscosidade, di�cultandoassim o tratamento de problemas com valores altos do número de Rey-nolds e d) a presença de ruídos ou �utuações.

A principal estratégia adotada para contornar a presença das �u-tuações foi a utilização de médias temporais, espaciais ou de ensemble.Como média de ensemble entende-se a média dos resultados de váriassimulações do mesmo problema com condições iniciais microscópicasdistintas, porém com condições macroscópicas idênticas. Esta aborda-gem acarreta em um aumento signi�cativo do custo computacional, o

50

que motivou propostas alternativas que deram origem ao método LB.

4.2 ORIGEM DA EQUAÇÃO LB

Visando eliminar os ruídos que ocorrem no Lattice Gas, Mc-Namara e Zanetti[29] propuseram um método alternativo que simuladiretamente a média de ensemble, ou seja, a média da ocupação de umsítio fi = ⟨ni⟩. O procedimento adotado pelos autores na obtenção daequação para fi a partir do LGA é similar ao realizado para dedução daequação de Boltzmann a partir das leis da mecânica newtoniana. Poreste motivo, o novo método foi denominado Lattice Boltzmann Method(LBM) ou Método de Boltzmann para redes.

Neste caso, a grandeza calculada, fi(x, t), é um número realque representa a fração de partículas que se encontram no sítio x comvelocidade ci no tempo t e é conhecida como a função distribuiçãodiscreta. A equação que rege a evolução do LBM, denominada de latticeBoltzmann equation (LBE) ou equação de Boltzmann para redes, podeser escrita na forma geral

fi(x+ hci, t+ δ) = fi(x, t) + δΩi. (4.3)

Nesta equação estão sintetizados os processos de colisão e pro-pagação herdados do LGA, onde o parâmetro h é a distância entre ossítios adjacentes e δ é o tempo que transcorre no problema físico du-rante um passo na simulação. O cálculo no lado direito da equaçãorepresenta a colisão, enquanto a propagação consiste na atribuição doresultado desse cálculo para o lado esquerdo. O termo Ωi é conhecidocomo operador de colisão, novamente em analogia com a equação deBoltzmann da teoria cinética.

No LBM, a massa especí�ca e a quantidade de movimento adi-mensionais são calculados, respectivamente, por meio das relações

ϱ =b∑

i=0

fi, (4.4)

ϱu =

b∑i=0

fici. (4.5)

A primeira versão do método de Boltzmann para redes herdoudiversas características do seu antecessor, tanto positivas, como a es-tabilidade numérica e a separação em etapas de colisão e propagação,

51

quanto negativas, como a dependência da pressão com a velocidade, aausência de invariância galileana e a limitação nos valores da viscosi-dade.

Porém, nesta transição do LGA para o LBM, perde-se um dosaspectos mais notáveis do modelo booleano, o baixo custo computa-cional. Isto se deve, em grande parte, a dois fatores: ao aumento daquantidade de memória necessária acarretado pelo uso de variáveis deponto �utuante e ao aumento do tempo de processamento requeridopela etapa de colisão, haja visto que este operador precisa computara contribuição de todos os possíveis cenários de colisão. Estes fatorestornaram praticamente inviável a aplicação desta versão inicial do LBMpara problemas tridimensionais.

O operador de colisão na forma inicialmente sugerida na Ref. [29]consiste em uma média do termo que trata das colisões no LGA, o queresulta em um operador não linear em fi. Por este motivo, a equação re-sultante �cou conhecida como LBE não linear. Um trabalho posterior,apresentado por Higuera e Zanetti[30], concentrou-se na simpli�caçãodeste operador. O caminho escolhido pelos autores para realizar talsimpli�cação foi expandir o operador de colisão em torno da distribui-ção encontrada no equilíbrio. Desprezando termos de segunda ordemnessa expansão, os autores chegaram a uma forma simpli�cada para ooperador de colisão,

Ωi =

b∑j=0

Aij(fj − feqj ) (4.6)

onde Aij permanece constante durante a simulação e pode ser deter-minado a partir do operador de colisão da LBE não linear. A equaçãoque resulta do uso desse operador foi chamada de LBE quase-linear,pois embora tenha uma aparência linear em fi, esta ainda apresentatermos não lineares, uma vez que feqj depende de potências de u, quepor sua vez dependem de fi. Mesmo com a permanência de termos nãolineares, este procedimento tornou a etapa de colisão mais e�ciente,permitindo o uso do LBM em problemas tridimensionais.

Motivados por este operador quase-linear, Chen et al.[31] suge-riram utilizar uma forma diagonal para a matriz Aij , resultando naseguinte expressão para o operador de colisão:

Ωi = −1

τ(fi − feqi ), (4.7)

onde τ é um parâmetro do modelo chamado de tempo de relaxação,

52

que permite controlar os coe�cientes de transporte. Este operador decolisão é o análogo discreto do modelo BGK encontrado na Eq. (3.12),e por este motivo a equação LB resultante do uso desse operador éconhecida como Lattice BGK Equation (LBGK).

Independentemente, Qian, d`Humières e Lallemand[32] realiza-ram uma abordagem semelhante. A contribuição especí�ca destes au-tores foi notar que a distribuição de equilíbrio utilizada no operadorBGK não precisa ser necessariamente aquela herdada do método LGA.Partindo desta constatação, estes sugeriram uma nova forma para adistribuição de equilíbrio

feqi = ϱwi

[1 + 3(ci · u) +

9

2(ci · u)2 −

3

2u2

]. (4.8)

Esta expressão foi obtida pelos autores utilizando uma distribuição deequilíbrio em série de potências até segunda ordem da variável u, sendoos coe�cientes dessa expansão escolhidos de modo que o comportamentomacroscópico corresponda aquele das equações de Navier-Stokes.

Como consequência dessa alteração, as redes formadas pelo con-junto de velocidades ci não �cam mais restritas aos requisitos de iso-tropia do LGA, sendo os pesos wi escolhidos de maneira a assegurar aisotropia das equações macroscópicas. A Eq. (4.8) pode ser utilizadacom redes unidimensionais, bidimensionais ou tridimensionais, sendoestas denotadas pela sigla DdQb, onde b é o número de velocidade e dé a dimensão. As principais redes deste grupo podem ser visualizadasna Figura 6, sendo que os pesos wi são encontrados no Apêndice B.

D3Q19D2Q9D1Q3

Figura 6: Redes utilizadas por Qian, d`Humières e Lallemand[32]

Esta nova abordagem trouxe vantagens como a eliminação dostermos que dão origem à ausência de invariância galileana e à depen-dência da pressão com a velocidade. Estes modelos se tonaram ampla-mente utilizados, em virtude principalmente da simplicidade do opera-

53

dor BGK e das vantagens aqui citadas.

4.3 CONEXÃO COM A TEORIA CINÉTICA

Conforme descrito na seção anterior, o LBM foi inicialmente pro-posto como um aprimoramento dos métodos de gás em rede. Ressalta-se ainda que nesta primeira abordagem não havia uma conexão diretaentre a equação de Boltzmann e a LBE. A similaridade na nomencla-tura destas equações se motivou por estas serem análogas: A primeiradescrevendo a evolução da média de essemble de um sistema compostopor partículas que seguem as leis de Newtown, e a segunda, a evoluçãodesta média para um mundo microscópico simpli�cado e discreto, cujadinâmica é ditada pelas regras do LGA.

Em 1997, He e Luo[33] demonstraram pela primeira vez que arelação entre a LBE e a equação de Boltzmann é mais profunda doque uma simples analogia na derivação de ambas. Estes autores con-seguiram obter a LBE prosta por Qian, d`Humières e Lallemand[32]

através da aplicação de técnicas numéricas especí�cas a equação deBoltzmann com o modelo BGK. A técnica proposta consistia em umprocesso de discretização combinada das variáveis tempo, posição e ve-locidade, além de uma expansão em série da função distribuição deequilíbrio. A partir do trabalho de He e Luo[33], determinadas formasda LBE passaram a ser interpretadas como formas discretas da equa-ção de Boltzmann com o operador BGK, o que gerou a questão se talprocedimento poderia ser aprimorado com o objetivo de obter novosmodelos numéricos.

Posteriormente, He, Chen e Doolen[34] apontaram uma falhanesta primeira derivação, onde nem todos os termos de primeira ordemno passo de tempo δ eram considerados, propondo uma correção no pro-cedimento original. Além desta discordância, na forma concebida porHe e Luo[33] o procedimento de discretização do espaço de velocidadesera restrito a obter as redes e distribuições de equilíbrio já existentes,sem apontar uma forma mais geral para realizar este procedimento quepermitisse melhorar sistematicamente a precisão da equação discreta.

De modo independente, Shan, Yuan e Chen[8] e Philippi et al.[7]

estenderam o processo de discretização do espaço de velocidades, permi-tindo melhorar de forma sistemática e consistente a ordem de descriçãoda distribuição de equilíbrio. No trabalho de Shan e colaboradores, ofoco é desenvolver modelos capazes de simular escoamentos com altosvalores do número de Knudsen, ou seja, para gases rarefeitos ou em

54

pequenas escalas de comprimento. Já o trabalho de Philippi et al temcomo foco a derivação de modelos para problemas não-isotérmicos demaneira consistente com a teoria cinética.

Por ser central para o presente trabalho, o procedimento de dis-cretização é apresentado em detalhe nesta seção. Vale ressaltar que, de-vido a evolução na forma e entendimento do processo de discretização,a abordagem utilizada aqui consiste no e