Tópicos Complementares

Disciplina: Métodos Numéricos em Termofluidos

Professor: Adriano Possebon Rosa

Departamento de Engenharia Mecânica

Faculdade de Tecnologia
Universidade de Brasília

Sumário
1 Introdução 2

2 Problema da Convecção Natural 2

2.1 Equação da Energia . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Aproximação de Boussinesq . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Equações Governantes . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.4 Adimensionalização . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.5 Implementação Numérica: Método de Projeção . . . . . . . . . . . . 9

3 Entrada e Saída de Fluido 10

4 Problema do Degrau 12

5 Método de Projeção de Primeira Ordem Implícito 12

6 Método de Projeção de Segunda Ordem 18

7 Formulação Vorticidade-Função de Corrente 20

7.1 Equações Governantes . . . . . . . . . . . . . . . . . . . . . . . . . . 20
7.2 Aproximação Usando Diferenças Finitas . . . . . . . . . . . . . . . . 21
7.3 Restrições . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

8 Cálculo da Ordem do Método 23

1
1 Introdução
Nesta aula são discutidos alguns tópicos complementares que serão a base para a
elaboração do artigo científico. Os tópicos são abordados de forma introdutória e
em cada um deles são apresentadas referências para aqueles leitores que queiram se
aprofundar no tema.

2 Problema da Convecção Natural

O objetivo aqui é apresentar o problema de convecção natural bidimensional
para o caso de uma cavidade quadrada de lado L, representada na figura (1). A
parede esquerda é mantida a uma temperatura fixa TH e a parede direita é mantida
a TC , com TH > TC , durante todo o processo. As paredes de cima e de baixo
são isoladas termicamente, ou seja, Q̇ = 0 nessas paredes. Resolver esse problema
significa obter as distribuições de velocidade e temperatura dentro da cavidade,
assim como obter o número de Nusselt. Vamos começar com a dedução das equações
governantes.

y
Q̇ = 0

L TH TC

Q̇ = 0 x

L
Figura 1: Convecção Natural em uma Cavidade Quadrada.

2
2.1 Equação da Energia
A energia de um sistema é dada por

m(u · u)
E = Eint + , (1)
2
em que m é a massa, u é a velocidade e

Eint = mcT (2)

representa a energia interna, onde T é a temperatura e c é o calor específico do
material. A energia específica, ou energia por unidade de massa, é dada por

u·u u2
e = cT + = cT + . (3)
2 2
Primeira lei da Termodinâmica para um sistema: a taxa de variação de energia do
sistema é igual ao calor que chega ao sistema somado ao trabalho realizado pela
vizinhança sobre o sistema. Na forma de equação, temos

dE
= Q̇ + Ẇ , (4)
dt
em que Q̇ é a quantidade de calor que chega ao sistema por unidade de tempo e Ẇ
é o trabalho realizado pela vizinhança sobre o sistema. Assim, temos
Z Z
Q̇ = − (q̇ · n̂) dS = − (∇ · q̇) dV . (5)
S V

Na equação (5) foi usado o teorema da divergência par obter a segunda igualdade.
Aqui, q̇ é o vetor fluxo de calor. Usando a lei de Fourier

q̇ = −k∇T (6)

na equação (5) e considerando a condutividade térmica do fluido, k, constante,

resulta:
Z Z
Q̇ = ∇ · (k∇T ) dV = k∇2 T dV . (7)
V V

O trabalho, por sua vez, é dado por

Z Z
Ẇ = (ρg) · u dV + (n̂ · σ) · u dS . (8)
V S

σ é o tensor de tensões, e representa as forças superficiais que são aplicadas no

sistema, por unidade de área. A componente σij atua no plano i e na direção j.
Para um fluido Newtoniano,

σ = −pI + 2µD , (9)

3
em que p é a pressão, I é o tensor identidade, µ é a viscosidade do fluido (constante)
e
1h i
D= ∇u + (∇u)T (10)
2
é o tensor taxa de deformação. Usando o teorema da divergência:
Z Z
n̂ · (σ · u) dS = ∇ · (σ · u) dV . (11)
S V

Então, o trabalho pode ser reescrito como:

Z Z
Ẇ = (ρg) · u dV + ∇ · (σ · u) dV . (12)
V V

A energia do sistema é:
Z
E= ρe dV . (13)
V

Assim, a primeira lei da termodinâmica pode ser escrita como:

Z
dE d
= ρe dV = Q̇ + Ẇ . (14)
dt dt V

Usando o Teorema do Transporte de Reynolds:

Z Z   Z
d ∂e De
ρe dV = ρ + ρu · ∇e dV = ρ dV . (15)
dt V V ∂t V Dt
Assim:
Z Z Z Z
De
ρ dV = 2
k∇ T dV + (ρg) · u dV + ∇ · (σ · u) dV . (16)
V Dt V V V

Ou:
Z  
De
ρ 2
− k∇ T − (ρg) · u − ∇ · (σ · u) dV = 0 . (17)
V Dt
Como a integral na equação (17) é zero para qualquer limite de integração, podemos
concluir que o integrando é nulo sempre, pelo Teorema da Localização. Assim:

De
ρ = k∇2 T + (ρg) · u + ∇ · (σ · u)
Dt (18)
= k∇2 T + (ρg) · u + u · (∇ · σ) + σ : (∇u) .

Substituindo a expressão para a energia específica:

DT ρ Du2
ρc + = k∇2 T + (ρg) · u + u · (∇ · σ) + σ : (∇u) . (19)
Dt 2 Dt

4
A equação de Cauchy é dada por:

Du
ρ = ρg + ∇ · σ . (20)
Dt
Fazendo o produto escalar da equação de Cauchy por u resulta:

ρ Du2
= u · ρg + u · (∇ · σ) . (21)
2 Dt
Substituindo (21) em (19):

DT
ρc + u · ρg + u · (∇ · σ) = k∇2 T + (ρg) · u + u · (∇ · σ) + σ : (∇u) . (22)
Dt
Assim:

DT
ρc = k∇2 T + σ : (∇u) . (23)
Dt
Considerando um fluido Newtoniano, temos

DT
ρc = k∇2 T + (−pI + 2µD) : (∇u) . (24)
Dt
Ou:

DT
ρc = k∇2 T − p(∇ · u) + 2µD : ∇u . (25)
Dt
Mas D : ∇u = D : D. Considerando um fluido incompressível (∇ · u = 0), temos:

DT
ρc = k∇2 T + 2µ (D : D) . (26)
Dt
O último termo dessa equação,

Φ = 2µDij Dij , (27)

é chamado de função de dissipação viscosa. Em muitos escoamentos, Φ é muito
menor do que os outros termos da equação (neste trabalho vamos assumir essa
hipótese). Portanto, a equação para a energia é dada por:
 
∂T
+ u · ∇T = α∇2 T . (28)
∂t
Aqui, α = k/(ρc) é a difusividade térmica e tem unidade de m2 /s no SI.
Comentários:

• a velocidade interfere no campo de temperatura, pois aparece no termo advec-

tivo da equação (28);

5
 • o trabalho total do tensor de tensões é dado por:

∇ · (σ · u) = u · (∇ · σ) + σ : (∇u) (29)

O primeiro termo do lado direito da equação contribui para uma variação na

energia cinética do elemento, como pode ser observado na equação (21). O
segundo termo contribui para um aumento na energia interna, ou seja, para
um aumento de temperatura, como pode ser visto na equação (23).

Número de Nusselt. A lei de resfriamento por convecção de Newton, para o caso

estudado no presente trabalho, é dada por:

Q̇ = h (TH − TC ) , (30)
em que Q̇ é a taxa de transferência de calor por convecção por unidade de área e h
é o coeficiente de transferência de calor por convecção. O número de Nusselt N u é
o coeficiente de transferência de calor adimensionalizado:

hLc
Nu = . (31)
k
Como nas paredes o fluido está em repouso, o calor é transmitido apenas por con-
dução, ou seja,

∂T
Q̇ = −k , (32)
∂x x=0

na parede esquerda, por exemplo. Q̇ = 0 nas paredes de cima e de baixo, já que estas
são termicamente isoladas. Esse calor que chega da parede esquerda por condução
é transferido por convecção no fluido, de maneira que temos a igualdade:

∂T
Q̇ = −k = h (TH − TC ) . (33)
∂x x=0

Isolando h:

k ∂T
∂x
h=− x=0
. (34)
(TH − TC )
Substituindo (34) em (31) resulta:

k ∂T
∂x Lc
Nu = − x=0
. (35)
(TH − TC ) k
No caso em que o fluido está em repouso, ou seja, temos transferência apenas por
condução pelo meio, resulta que N u = 1.

6
2.2 Aproximação de Boussinesq
A variação de temperatura leva a uma variação de massa específica que faz aparecer
uma força de empuxo. Precisamos relacionar a variação de temperatura T com a
variação de massa específica ρ. A propriedade que nos dá essa relação é o coeficiente
de expansão volumétrica β, definido como:
 
1 ∂ρ
β=− . (36)
ρ ∂T P

β é uma propriedade do fluido e o subíndice P indica que essa definição é válida para
pressão constante. Fazendo uma aproximação da equação (36), podemos reescrevê-la
como:

1 ∆ρ 1 (ρ − ρo )
β≈− =− . (37)
ρo ∆T ρo (T − To )
ρo e To representam a massa específica e a temperatura do fluido longe da parede (ou
em alguma configuração de referência), respectivamente. Assim, ρ e T são pequenos
desvios de ρo e To . Isolando ρ em (37) temos:

ρ = ρo − βρo ∆T . (38)
A equação de Navier-Stokes é:
 
∂u
ρo + u · ∇u = −∇p + µ∇2 u + ρg . (39)
∂t
Substituindo a relação (38) na equação (39), resulta:
 
∂u
ρo + u · ∇u = −∇p + µ∇2 u + ρo g − ρo gβ (T − To ) . (40)
∂t
Considerando que a aceleração gravitacional é dada por

g = −gêy , (41)
e definindo uma nova pressão como sendo

p∗ = p + ρo gy , (42)
temos
 
∂u
ρo + u · ∇u = −∇p + µ∇2 u + ρo gβ (T − To ) êy . (43)
∂t
Nessa equação já foi usada a nova pressão (sem o asterisco). Note que a força extra
que aparece atua apenas na direção y , que é a direção da aceleração gravitacional.
Ao conjunto de passos que leva à equação (43) é dado o nome de aproximação de
Boussinesq [1, 2].

7
2.3 Equações Governantes
As equações governantes desse problema são a equação da conservação de massa
(fluido incompressível),

∇·u=0 , (44)
a equação de Navier-Stokes com a correção de Boussinesq,
 
∂u
ρo + u · ∇u = −∇p + µ∇2 u + ρo gβ (T − To ) êy , (45)
∂t
e a equação da energia (temperatura),
 
∂T
+ u · ∇T = α∇2 T . (46)
∂t
Nas equações, α = k/(ρo c), To é uma temperatura de referência e ρo é a massa
específica do fluido nessa temperatura. No nosso problema específico, a temperatura
de referência é TC .

2.4 Adimensionalização
Para adimensionalizar as equações governantes do problema, (44), (45) e (46), va-
mos utilizar L como um comprimento característico, α/L como uma velocidade
característica, L2 /α como um tempo característico e ρo α2 /L2 como uma pressão
característica. Assim, as grandezas características usadas para adimensionalizar as
equações são:

α L2 α2
Lc = L , uc = , tc = e pc = ρ o . (47)
L α L2
As variáveis adimensionais são:
u L α pL2
T − TC
∇∗ = L∇ u∗ = =u , t∗ = t p∗ = p , . θ=
uc α L2 ρo α 2
TH − TC
(48)
Note que θ é a temperatura adimensional. Substituindo essas variáveis de volta nas
equações governantes resulta:

∇∗ · u∗ = 0 , (49)

∂u∗
+ u∗ · ∇∗ u∗ = −∇∗ p∗ + P r∇∗2 u∗ + RaP rθêy e (50)
∂t∗

∂θ
∗
+ u∗ · ∇∗ θ = ∇∗2 θ . (51)
∂t
Na equação (50),

8
 ν µ/ρo µc
Pr = = = (52)
α k/(ρo c) k
é o número de Prandtl e

gβ (TH − TC ) L3c
Ra = (53)
να
é o número de Rayleigh.

2.5 Implementação Numérica: Método de Projeção

As equações que devem ser resolvidas são (já na forma adimensional)

∇·u=0 , (54)

∂u
+ u · ∇u = −∇p + P r∇2 u + RaP rθêy e (55)
∂t

∂θ
+ u · ∇θ = ∇2 θ . (56)
∂t
A figura (2) mostra o problema proposto na forma adimensional. A velocidade possui
valor zero em todas as paredes. Para a temperatura, as condições de contorno são:

• θ = 1 para a parede esquerda;

• θ = 0 para a parede direita;

• Q̇ = ∂θ
∂y
= 0 paras as paredes de cima e de baixo.

Uma maneira de resolver numericamente essas equações governantes é por meio do

método de projeção, acoplado com o método de diferenças finitas para aproximar
as derivadas. Os passos para resolver as equações são:

−u
 u ∆t
 ∗ k
+ uk · ∇uk = P r∇2 uk + RaP rθk êy
PASSO 1 (57)
u = 0 em ∂Ω
∗


 θk+1 −θk

 ∆t
+ uk · ∇θk = ∇2 θk



PASSO 2 θ = 1 em x = 0 ; θ = 0 em x = 1 (58)



= 0 em y = 0 e y = 1

 ∂θ
∂y
 2 k+1
 ∇p = ∆t1
∇ · u∗
PASSO 3 (59)
∇pk+1 · n̂ = 0 em ∂Ω


9
 y
Q̇ = 0

1 1 0

Q̇ = 0 x

1
Figura 2: Convecção Natural em uma Cavidade Quadrada. Problema adimensional.

PASSO 4 uk+1 = u∗ − ∆t∇pk+1 (60)



Posição da temperatura: usando uma malha defasada, a temperatura pode ser

posicionada no centro do elemento, no mesmo local onde se encontra a pressão.
Número de Nusselt. Na forma adimensional, a equação (35) é dada por:

∂θ
Nu = − . (61)
∂x x=0

Esse é o número de Nusselt local, ou seja, para cada valor y da parede esquerda. O
número de Nusselt médio N um é dado por:
Z 1
∂θ
N um = − dy . (62)
0 ∂x x=0

3 Entrada e Saída de Fluido

Muitos problemas em fluidos envolvem a entrada e a saída de fluido de um domínio.
Para entendermos como essas condições de contorno são implementadas, vamos con-
siderar o escoamento de um fluido entre duas placas planas paralelas, como mostra a

10
figura (3). Vamos considerar como base a metodologia desenvolvida para o problema
da cavidade cisalhante. As paredes de cima e de baixo são sólidas e estão paradas.
As condições de contorno aqui serão iguais às da cavidade, sendo necessário apenas
zerar a velocidade da tampa. Ou seja, para as paredes sólidas, temos

u=v=0

e
∂p
=0.
∂x

∂p
u = 1; v = 0; ∂x =0

y
p = 0; ∂u
∂x = 0
Entrada Saı́da

1
x
L
Figura 3: Escoamento entre placas planas, com entrada e saída de fluido.

O que era a parede esquerda da cavidade agora se torna a entrada de fluido.

Na entrada nós vamos especificar a velocidade e “liberar” a pressão. Dessa forma,
vamos impor
u = 1 ,v = 0
e
∂p
=0.
∂x
Assim, a condição de contorno para a pressão na entrada é a mesma da parede sólida.
Assim, não é preciso alterar a equação da pressão nessa parede. Para v também não
é necessário alterar nada, já que continua zero. É preciso alterar apenas as condições
para u, já que agora u = 1.
Na saída, nós vamos fazer o oposto: vamos “liberar” as velocidades e vamos fixar
a pressão. Assim, na saída as condições de contorno são
∂u ∂v
= =0
∂x ∂x
e
p=0.
Agora é necessário modificar as equações das 3 incógnitas nas proximidades da saída.
Muito cuidado com a pressão: é necessário alterar as equações do método SOR.

11
 Para testar, faça simulações com Re = 100 e L = 10. Use 100 pontos na direção
x e 10 na direção y, com um ∆t de 0.01. Você vai observar a região de entrada
e depois uma convergência do escoamento para o perfil parabólico. Compare seus
resultados com a solução analítica.

4 Problema do Degrau
Um problema com entrada e saída bem interessante é o problema do degrau, re-
presentado na figura (4). Agora, na superfície da esquerda, temos metade sendo
entrada e metade sendo parede sólida. Nesse caso você deve aplicar as condições
de contorno de parede para 0 ≤ y ≤ 0.5 e as condições de contorno de entrada de
fluido para 0.5 < y ≤ 1.
O resultado é bem complexo, com o aparecimento de uma recirculação em baixo,
próximo da entrada e, para Re maiores, uma segunda recirculação em cima, a jusante
da entrada. Um valor de interesse é o tamanho da primeira recirculação que se
forma após o degrau. Você pode comparar os seus resultados com outros resultados
numéricos da literatura [3, 4] ou com resultados experimentais [5]. Aqui L deve ser
pelo menos 15 para que você tenha bons resultados. As figuras (5) e (6) mostram o
resultado para Re = 400, sendo a segunda um gráfico de contorno da pressão.
∂p
u = 1; v = 0; ∂x =0

y
p = 0; ∂u
∂x = 0
Entrada Saı́da
0.5
0.5
x
L
Figura 4: Escoamento com degrau, com entrada e saída de fluido.

5 Método de Projeção de Primeira Ordem Implícito

Uma desvantagem do Método de Projeção de Primeira Ordem Explícito, que foi o
que vimos no estudo do Problema da Cavidade Cisalhante, é que há uma restrição
no passo de tempo ∆t causada pelo termo difusivo. Uma das maneiras de contornar
esse obstáculo é implementando a versão implícita desse método:

12
 2

1.5

1
Y

0.5

-0.5
xr
-1
0 5 10 15 20
X

Figura 5: Linhas de corrente do escoamento com degrau, para Re = 400. Foram

usados 400 pontos na direção x e 50 na direção y.
1.5

1
Y

0.5

0
xr

-0.5
0 5 10 15 20
X

Figura 6: Campo de pressão do escoamento com degrau, para Re = 400. Foram

usados 400 pontos na direção x e 50 na direção y.

 u ∆t
−u
 ∗ k
+ uk · ∇uk = 1
Re
∇ 2 u∗
PASSO 1 (63)
u = ub em ∂Ω
∗

 2 k+1
 ∇p = 1
∆t
∇ · u∗
PASSO 2 (64)
∇p k+1
· n̂ = 0 em ∂Ω


PASSO 3 uk+1 = u∗ − ∆t∇pk+1 (65)



Agora, porém, faz-se necessário resolver 3 sistemas de equações lineares em cada

passo de tempo: um para ustar , um para vstar e um para a pressão p. Nesta seção
vamos ver uma metodologia para implementar esse método, usando uma malha
defasada.
Para ustar temos:

13
    
ustar [i, j] − u[i, j] u[i + 1, j] − u[i − 1, j] u[i, j + 1] − u[i, j − 1]
+ u[i, j] + C1 =
∆t 2∆x 2∆y
"
1 ustar [i + 1, j] − 2ustar [i, j] + ustar [i − 1, j]
+
Re ∆x2
#
ustar [i, j + 1] − 2ustar [i, j] + ustar [i, j − 1]
,
∆y 2
(66)
em que
1h i
C1 = v[i, j + 1] + v[i − 1, j + 1] + v[i, j] + v[i − 1, j] . (67)
4
Reorganizando:
"
∆t ustar [i + 1, j] − 2ustar [i, j] + ustar [i − 1, j]
ustar [i, j]− +
Re ∆x2
#
ustar [i, j + 1] − 2ustar [i, j] + ustar [i, j − 1]
=
∆y 2
(   (68)
u[i + 1, j] − u[i − 1, j]
= u[i, j] − ∆t u[i, j]
2∆x
  )
u[i, j + 1] − u[i, j − 1]
+ C1 .
2∆y

Ou:
 
2∆t 2∆t
ustar [i, j] 1 + + =
Re∆x2 Re∆y 2
(    )
u[i + 1, j] − u[i − 1, j] u[i, j + 1] − u[i, j − 1]
= u[i, j] − ∆t u[i, j] + C1
2∆x 2∆y
" #
∆t ustar [i + 1, j] + ustar [i − 1, j] ustar [i, j + 1] + ustar [i, j − 1]
+ + .
Re ∆x2 ∆y 2
(69)

Com isso, conseguimos isolar ustar [i, j]. Usando o método iterativo de Gauss-
Seidel para resolver um sistema de equações lineares, podemos escrever, após algu-
mas operações de soma e subtração de ustar no lado direito da equação (69),

ustar [i, j] ← ustar [i, j] + R , (70)

em que R é o incremento em cada iteração para encontrar o ustar correto. A partir

14
da equação (69), R é dado por:
("    #
u[i + 1, j] − u[i − 1, j] u[i, j + 1] − u[i, j − 1]
R = λ u[i, j] − ∆tu[i, j] − ∆tC1
2∆x 2∆y
"
∆t ustar [i + 1, j] − 2ustar [i, j] + ustar [i − 1, j]
− ustar [i, j] − +
Re ∆x2
!#)
ustar [i, j + 1] − 2ustar [i, j] + ustar [i, j − 1]
,
∆y 2
(71)
em que  −1
2∆t 2∆t
λ= 1+ 2
+ . (72)
Re∆x Re∆y 2
As condições de contorno são aplicadas diretamente à ustar . Assim, ustar [i, −1] =
−ustar [i, 0] e ustar [i, Ny ] = 2U [i] − ustar [i, Ny − 1]. Essas condições devem ser substi-
tuídas na equação (71), com atenção ao fato de que o valor de λ também se altera.
Os detalhes podem ser encontrados no algoritmo (1).

O procedimento para a obtenção de vstar é análogo. A equação, em sua versão

implícita, é dada por:
   
vstar [i, j] − v[i, j] v[i + 1, j] − v[i − 1, j] v[i, j + 1] − v[i, j − 1]
+ C2 + v[i, j] =
∆t 2∆x 2∆y
"
1 vstar [i + 1, j] − 2vstar [i, j] + vstar [i − 1, j]
+
Re ∆x2
#
vstar [i, j + 1] − 2vstar [i, j] + vstar [i, j − 1]
.
∆y 2
(73)
Nessa equação, C2 corresponde ao valor de u interpolado:
1h i
C2 = u[i + 1, j] + u[i, j] + u[i + 1, j − 1] + u[i, j − 1] . (74)
4
Reescrevendo:
"
∆t vstar [i + 1, j] − 2vstar [i, j] + vstar [i − 1, j]
vstar [i, j]− +
Re ∆x2
#
vstar [i, j + 1] − 2vstar [i, j] + vstar [i, j − 1]
=
∆y 2
(    )
v[i + 1, j] − v[i − 1, j] v[i, j + 1] − v[i, j − 1]
= v[i, j] − ∆t C2 + v[i, j] .
2∆x 2∆y
(75)

15
Ou:
 
2∆t 2∆t
vstar [i, j] 1 + + =
Re∆x2 Re∆y 2
(    )
v[i + 1, j] − v[i − 1, j] v[i, j + 1] − v[i, j − 1]
= v[i, j] − ∆t C2 + v[i, j]
2∆x 2∆y
" #
∆t vstar [i + 1, j] + vstar [i − 1, j] vstar [i, j + 1] + vstar [i, j − 1]
+ + .
Re ∆x2 ∆y 2
(76)

Com isso, conseguimos isolar vstar [i, j]. Usando o método iterativo de Gauss-
Seidel para resolver um sistema de equações lineares, podemos escrever, após algu-
mas operações de soma e subtração de vstar no lado direito da equação (76),

vstar [i, j] ← vstar [i, j] + R , (77)

em que R é o incremento em cada iteração para encontrar o vstar correto. A partir

da equação (76), R é dado por:
("    #
v[i + 1, j] − v[i − 1, j] v[i, j + 1] − v[i, j − 1]
R = λ v[i, j] − ∆tC2 − ∆tv[i, j]
2∆x 2∆y
"
∆t vstar [i + 1, j] − 2vstar [i, j] + vstar [i − 1, j]
− vstar [i, j] − +
Re ∆x2
!#)
vstar [i, j + 1] − 2vstar [i, j] + vstar [i, j − 1]
,
∆y 2
(78)
em que  −1
2∆t 2∆t
λ= 1+ 2
+ . (79)
Re∆x Re∆y 2
As condições de contorno são aplicadas diretamente à vstar . Assim, vstar [−1, j] =
−vstar [0, j] e vstar [Nx , j] = −vstar Nx − 1, j]. Essas condições devem ser substituídas
na equação (78), com atenção ao fato de que o valor de λ também se altera. Os
detalhes podem ser encontrados no algoritmo (2).
No caso do método implícito apresentado nesta seção, o algoritmo completo para
resolver o problema da cavidade é igual ao algoritmo (??), com a diferença de que
agora ustar é dado pelo algoritmo (1) e vstar pelo algoritmo (2).

16
Algorithm 1 Calculate ustar
INPUT: ustar , u, v, Nx , Ny , ∆x, ∆y, ∆t, tol
OUTPUT: ustar
error = 100
while error > tol do
R-max = 0
for i = 1, Nx - 1 do
for j = 0, hNy - 1 do i
1
C1 = v[i, j + 1] + v[i − 1, j + 1] + v[i, j] + v[i − 1, j]
4
if j = 0hthen i −1
2∆t 3∆t
λ = 1 + Re∆x 2 + Re∆y 2
(" #
   
u[i+1,j]−u[i−1,j] u[i,j+1]−u[i,j−1]
R=λ u[i, j]−∆tu[i, j] 2∆x
−∆tC1 2∆y
"
∆t ustar [i+1,j]−2ustar [i,j]+ustar [i−1,j]
− ustar [i, j] − Re ∆x2
+
!#)
+ ustar [i,j+1]−3u
∆y 2
star [i,j]

else if j = Ny − 1 then
h i−1
2∆t 3∆t
λ = 1 + Re∆x2 + Re∆y2
(" #
   
R = λ u[i, j]−∆tu[i, j] u[i+1,j]−u[i−1,j]
2∆x
−∆tC1 u[i,j+1]−u[i,j−1]
2∆y
"
∆t ustar [i+1,j]−2ustar [i,j]+ustar [i−1,j]
− ustar [i, j] − Re ∆x2
+
!#)
+ 2U [i]−3ustar [i,j]+u
∆y 2
star [i,j−1]

else h i−1
2∆t 2∆t
λ = 1 + Re∆x 2 + Re∆y 2
(" #
   
u[i+1,j]−u[i−1,j] u[i,j+1]−u[i,j−1]
R=λ u[i, j]−∆tu[i, j] 2∆x
−∆tC1 2∆y
"
∆t ustar [i+1,j]−2ustar [i,j]+ustar [i−1,j]
− ustar [i, j] − Re ∆x2
+
!#)
+ ustar [i,j+1]−2ustar
∆y 2
[i,j]+ustar [i,j−1]

end if
ustar [i, j] ← ustar [i, j] + R
if |R| > R-max then
R-max ← |R|
end if
end for
end for
error ← R-max
end while 17
 BOUNDARY UPDATE:
for j = 0, Ny -1 do
ustar [0, j] = 0
ustar [Nx , j] = 0
end for
for i = 0, Nx do
ustar [i, −1] = −ustar [i, 0]
ustar [i, Ny ] = 2U [i] − ustar [i, Ny − 1]
end for

6 Método de Projeção de Segunda Ordem

O erro associado aos métodos de projeção de primeira ordem, explícito e implícito,
são de primeira ordem no tempo e de segunda ordem no espaço. Nesta seção é
apresentado um método de projeção de segunda ordem tanto no tempo quanto no
espaço. Este método foi desenvolvido por Bell, Colella e Glaz [6] e modificado por
Brown, Cortez e Minion [7]. Os passos para a implementação desse método são:

 u ∆t
−u
 ∗ k
1
u∗ + uk − ∇pk−1/2



+ uk+1/2 · ∇ uk+1/2 = 2Re
∇2
PASSO 1 (80)
u = ub em ∂Ω
∗


 2 k+1
 ∇φ = 1
∆t
∇ · u∗
PASSO 2 (81)
∇φ k+1
· n̂ = 0 em ∂Ω


PASSO 3 uk+1 = u∗ − ∆t∇φk+1 (82)



PASSO 4 ∆t
(83)

pk+1/2 = pk−1/2 + φk+1 − 2Re
∇2 φk+1
Note que agora existe uma função φ, que tem um comportamento parecido com a
pressão p. Os termos advectivos no passo 1 são obtidos por meio de uma interpolação
no tempo dada por:
3 1
f k+1/2 = f k − f k−1 . (84)
2 2
Dessa forma, apenas o termo difusivo fica implícito na equação (80) e o sistema de
equações algébricas resultante é linear. Três sistemas de equações lineares devem
ser resolvidos em cada passo de tempo: um para ustar , um para vstar e um para φ.
A pressão é obtida em um instante de tempo intermediário entre k e k + 1.

18
Algorithm 2 Calculate vstar
INPUT: u, v, Nx , Ny , ∆x, ∆y, ∆t, tol
OUTPUT: vstar
error = 100
while error > tol do
R-max = 0
for i = 0, Nx - 1 do
for j = 1, hNy - 1 do i
1
C2 = u[i + 1, j] + u[i, j] + u[i + 1, j − 1] + u[i, j − 1]
4
if i = 0 hthen i−1
3∆t 2∆t
λ = 1 + Re∆x 2 + Re∆y 2
(" #
   
R = λ v[i, j] − ∆tC2 v[i+1,j]−v[i−1,j]
2∆x
− ∆tv[i, j] v[i,j+1]−v[i,j−1]
2∆y
"
∆t vstar [i+1,j]−3vstar [i,j]
− vstar [i, j] − Re ∆x2
!#)
vstar [i,j+1]−2vstar [i,j]+vstar [i,j−1]
∆y 2

else if ih= Nx − 1 then i

−1
3∆t 2∆t
λ = 1 + Re∆x 2 + Re∆y 2
(" #
   
R = λ v[i, j] − ∆tC2 v[i+1,j]−v[i−1,j]
2∆x
− ∆tv[i, j] v[i,j+1]−v[i,j−1]
2∆y
"
∆t −3vstar [i,j]+vstar [i−1,j]
− vstar [i, j] − Re ∆x2
!#)
vstar [i,j+1]−2vstar [i,j]+vstar [i,j−1]
∆y 2

else h i−1
2∆t 2∆t
λ = 1 + Re∆x 2 + Re∆y 2
(" #
   
v[i+1,j]−v[i−1,j] v[i,j+1]−v[i,j−1]
R=λ v[i, j] − ∆tC2 2∆x
− ∆tv[i, j] 2∆y
"
∆t vstar [i+1,j]−2vstar [i,j]+vstar [i−1,j]
− vstar [i, j] − Re ∆x2
!#)
vstar [i,j+1]−2vstar [i,j]+vstar [i,j−1]
∆y 2

end if
vstar [i, j] ← vstar [i, j] + R
if |R| > R-max then
R-max ← |R|
end if
end for
end for
error ← R-max
end while 19
 BOUNDARY UPDATE:
for j = 0, Ny do
vstar [−1, j] = −vstar [0, j]
vstar [Nx , j] = −vstar [Nx − 1, j]
end for
for i = 0, Nx − 1 do
vstar [i, 0] = 0
vstar [i, Ny ] = 0
end for

7 Formulação Vorticidade-Função de Corrente

Vamos abordar o problema da cavidade cisalhante quadrada, de lado L e tampa se
movendo com velocidade U0 , usando a formulação vorticidade-função de corrente.

7.1 Equações Governantes

As equações governantes para esse problema são a equação da continuidade
∂u ∂v
+ =0 (85)
∂x ∂y
e a equação de Navier-Stokes
Du
ρ = −∇p + µ∇2 u . (86)
Dt

Usando a intensidade de velocidade U0 e o tamanho L da cavidade, as equações

acima podem ser adimensionalizadas. O resultado é (todas as variáveis já estão na
forma adimensional):
∂u ∂v
+ =0 (87)
∂x ∂y
e
Du 1 2
= −∇p + ∇u, (88)
Dt Re
em que
ρU0 L
Re = (89)
µ
é o número de Reynolds. O domínio para essa forma adimensional é 0 ≤ x ≤ 1 e
0 ≤ y ≤ 1 e a velocidade da parede superior se torna U (x) = 1.
Em vez de resolvermos diretamente u, v e p, que são chamadas de variáveis pri-
mitivas, vamos resolver o problema nas variáveis vorticidade ω e função de corrente
Ψ. A equação da continuidade é satisfeita automaticamente se definirmos Ψ como:
∂Ψ ∂Ψ
u= e v=− . (90)
∂y ∂x

20
A vorticidade, no caso bidimensional, por sua vez, é dada por:

∂v ∂u ∂ 2Ψ ∂ 2Ψ
ω= − =− 2 − 2
= −∇2 Ψ . (91)
∂x ∂y ∂x ∂y
O rotacional da equação de Navier-Stokes fornece uma equação para a vorticidade:

1 ∂ 2ω ∂ 2ω
 
∂ω ∂ω ∂ω 1 2
+u +v = ∇ω= + 2 . (92)
∂t ∂x ∂y Re Re ∂x2 ∂y

Usando a definição de Ψ na equação acima, resulta:

∂ 2ω ∂ 2ω
 
∂ω ∂Ψ ∂ω ∂Ψ ∂ω 1 2 1
+ − = ∇ω= + 2 . (93)
∂t ∂y ∂x ∂x ∂y Re Re ∂x2 ∂y

Assim, temos duas equações (91 e 93) e duas incógnitas (ω e Ψ). Esse é o novo
sistemas de equações que iremos resolver usando diferenças finitas.

7.2 Aproximação Usando Diferenças Finitas

Para resolver numericamente as equações (91) e (93), o domínio é discretizado com
N espaços na direção x e N espaços na direção y. Os incrementos espaciais são
∆x = ∆y = 1/N . A função de corrente Ψ(x, y, t) é aproximada pela sua equivalente
numérica Ψ(i∆x, j∆y, k∆t) = Ψki,j , com 0 ≤ i, j ≤ N e k ≥ 0. O mesmo vale para
a vorticidade ω. Note que neste caso a malha não é deslocada: as variáveis ω e Ψ
estão no mesmo lugar do domínio numérico, que corresponde às quinas das células
computacionais.
Usando a aproximação de diferenças finitas, a equação da vorticidade se torna:

! !
k+1 k
ωi,j − ωi,j Ψki,j+1 − Ψki,j−1 k
ωi+1,j k
− ωi−1,j
+
∆t 2∆y 2∆x
! !
Ψki+1,j − Ψki−1,j k
ωi,j+1 k
− ωi,j−1
− = (94)
2∆x 2∆y
!
k k k k k k
1 ωi+1,j − 2ωi,j + ωi−1,j ωi,j+1 − 2ωi,j + ωi,j−1
= + .
Re ∆x2 ∆y 2

A equação que relaciona a função de corrente e a vorticidade, por sua vez, pode
ser escrita como:

Ψk+1 k+1 k+1

i+1,j − 2Ψi,j + Ψi−1,j Ψk+1 k+1 k+1
i,j+1 − 2Ψi,j + Ψi,j−1
+ k+1
= −ωi,j . (95)
∆x2 ∆y 2
Portanto, em cada passo tempo, primeiro calculamos diretamente o novo ω, que é o
k+1
ωi,j . Em seguida usamos a equação (95) e o novo valor de ω para calcular Ψk+1
i,j .

21
 A condição de contorno para Ψ é: Ψi,j = 0 para todos os pontos que estejam
em alguma das paredes. Assim, Ψi,j = 0 se i = 0, ou se i = N , ou se j = 0, ou se
j = N.
A condição inicial para Ψ é Ψ0i,j = 0 para todos os pontos do domínio.
Já a condição de contorno para ω é atualizada em todo passo de tempo e é
específica para cada parede:

Parede inferior (j = 0):

2 Ψki,0 − Ψki,1
k
ωi,0 = (96)
∆y 2

Parede superior (j = N ):


2 Ψki,N − Ψki,N −1 2Ui
k
ωi,N = 2
− (97)
∆y ∆y

Parede esquerda (i = 0):

2 Ψk0,j − Ψk1,j
k
ω0,j = (98)
∆x2

Parede direita (i = N ):

2 ΨkN,j − ΨkN −1,j
k
ωN,j = (99)
∆x2

Essas condições são válidas para todos os instantes de tempo k, para k ≥ 0. Assim,
a condição inicial para ω é:

se j =

0 6 N;
ωij0 = (100)
− 2U
∆y
i
se j = N.

7.3 Restrições
Para evitar instabilidades numéricas, devemos obedecer às seguintes restrições quanto
aos tamanhos do passo de tempo ∆t e do incremento espacial ∆x:

∆t < ∆x ,
1
∆x <
Re
e
1
∆t < Re∆x2 . (101)
4
22
8 Cálculo da Ordem do Método
É sempre bom checar os resultados que estamos obtendo numericamente. Não pode-
mos simplesmente confiar nos números que o computador exporta apenas porque é
um computador. É importante ter um bom conhecimento sobre o problema que está
sendo investigado, sobre como seria a solução esperada e também sobre quais se-
riam os tipos de resultados absurdos (concentração negativa, temperatura absoluta
negativa, velocidade muito alta, etc).
Um resultado numérico tem erros causados pela precisão de máquina e pela
aproximação discreta, além de poder ter erros humanos, gerados no momento da
implementação (erro de programação ou erro de conta, por exemplo). Por isso
precisamos apurar se o código implementado fornece resultados confiáveis.
Para saber se os resultados numéricos são válidos, podemos compará-los com
resultados analíticos, numéricos ou experimentais. Apenas por uma questão de
nomenclatura, quando comparamos com resultados analíticos (para casos simples)
ou numéricos (já consagrados na literatura), estamos verificando o nosso código. Já
quando comparamos com resultados experimentais estamos validando.
Uma forma de investigar os resultados e ver se eles estão se comportando bem é
por meio de uma análise de convergência de malha. Isso é feito a partir da resolução
de um mesmo problema com malhas cada vez mais finas: a solução deve convergir
para um valor constante à medida em que ∆t (ou ∆x) tende a zero.
Outra forma é a partir de uma análise da ordem do método implementado. A
ordem é a taxa na qual o erro do método numérico diminui na medida em que o
passo de tempo (ou o tamanho da malha) se aproxima de zero. Vimos, ao longo
do curso, diferentes métodos numéricos. Alguns eram O(∆t), outros O(∆t2 ). No
espaço, vimos métodos O(∆x) e O(∆x2 ). Na prática, podemos testar a ordem para
ver se a implementação foi feita de forma correta e se o código se comporta realmente
na ordem esperada.
Vamos ver como a ordem é calculada. Considere h como sendo o parâmetro
relacionado à discretização do domínio temporal ou espacial. Assim, h pode ser o
tamanho da malha (∆x ou ∆y) ou o tamanho do passo de tempo (∆t). A solução
exata para um dado problema é f , enquanto a solução numérica aproximada, usando
h como parâmetro, é dada por fh .
O erro é a diferença entre o valor numérico e o exato. Se o método é de ordem
p, podemos dizer que
|fh − f | ≈ Chp ,
em que C é uma constante. Dividindo h por dois e resolvendo o mesmo problema,
encontramos um outro resultado, que chamaremos de fh/2 . Assim,
 p
h Chp
|fh/2 − f | ≈ C = p .
2 2
Da mesma forma, para h/4 temos
 p
h Chp
|fh/4 − f | ≈ C = 2p .
4 2

23
 Nas raras situações em que conhecemos a solução exata do problema, podemos
calcular a ordem do método comparando as soluções usando h e h/2. Isso é feito da
seguinte forma:
fh − f Chp
≈ p p
= 2p . (102)
fh/2 − f Ch /2
Nas situações mais comuns, quando não temos uma solução exata, podemos
calcular a ordem do método usando os resultados das malhas h, h/2 e h/4:
p p
fh − fh/2 Chp − Ch Chp − Ch
≈ Chp Ch 2p
p = 1 Chp
2p

Chp
= 2p . (103)
fh/2 − fh/4 2p
− 22p 2 2p
− 22p

Assim, é possível calcular a ordem do método usando os resultados de 3 simulações

diferentes. Uma observação importante é que muitas vezes a ordem teórica só é
alcançada quando h é muito pequeno.
Todo o material desta seção foi retirado do livro de diferenças finitas do Leveque
[8].

24
Referências
[1] A. Bejan. Convection Heat Transfer. John Wiley and Sons, Inc., New Jersey,
4th edition, 2013.

[2] Y. Cengel and A. Ghajar. Heat and Mass Transfer. McGraw Hill Education,
New York, 5th edition, 2015.

[3] J. Kim and P. Moin. Application of a fractional-step method to incompressible

navier-stokes equations. Journal of Computational Physics, 59:308–323, 1985.

[4] G. Biswas, M. Breuer, and F. Durst. Backward-facing step flows for various
expansion ratios at low and moderate reynolds numbers. Journal of Fluid Engi-
neering, 126:362–374, 2004.

[5] B. F. Armaly, F. Durst, J. C. F. Pereira, and B. Schonung. Experimental and

theoretical investigation of backward-facing step flow. Journal of Fluid Mecha-
nics, 127:473–496, 1983.

[6] J. B. Bell, P. Colella, and H. M. Glaz. A Second-Order Projection Method for

the Incompressible Navier-Stokes Equations. Journal of Computational Physics,
85:257–283, 1989.

[7] D. L. Brown, R. Cortez, and M. L. Minion. Accurate Projection Methods for

the Incompressible Navier-Stokes Equations. Journal of Computational Physics,
168:464–499, 2001.

[8] R. J. Leveque. Finite Difference Methods for Ordinary and Partial Differential
Equations. SIAM, Philadelphia, 1st edition, 2007.

25