Aplicação de Métodos de Quarta Ordem para Resolver as
Equações de Navier-Stokes em um Canal com uma Oclusão
Katia Prado Fernandes
Programa de Pós Graduação em Modelagem Computacional, LNCC
25651-075, Petrópolis, RJ
E-mail: [email protected],
Paulo F. A. Mancera
Departamento de Bioestatı́stica, IBB, UNESP
18618-000, Botucatu, SP
E-mail: [email protected].
Resumo: Métodos numéricos de diferenças finitas de alta ordem para as equações de NavierStokes nas formulações função corrente e função corrente-vorticidade são classificados como
compactos ou não compactos. Neste trabalho, constrói-se um método compacto e um não compacto de quarta ordem para resolver numericamente as equações de Navier-Stokes na formulação
função corrente em uma malha uniforme. Apresenta-se os erros RMS e máximo e verifica-se que
o resultado numérico para a combinação de métodos de segunda ordem com métodos de quarta
ordem são mais precisos do que o método de diferença central de segunda ordem.
Palavras-chave: Equações de Navier-Stokes, Diferenças Finitas, Métodos de Alta Ordem
1
Introdução
Nas últimas décadas nota-se o desenvolvimento de vários métodos numéricos de diferenças
finitas para resolver numericamente as equações de Navier-Stokes, tanto em variáveis primitivas,
quanto na formulação função corrente-vorticidade. Na maioria das vezes os métodos usados são
de segunda ordem, resultando em uma molécula computacional de cinco pontos.
Pode-se obter métodos numéricos de alta ordem para as equações de Navier-Stokes usando
diferenças centrais de quarta ordem, porém os mesmos resultam em moléculas computacionais
enormes (ver Figura 1.a), ditos métodos não compactos. Entretanto, é possı́vel obter métodos
de diferenças finitas de alta ordem, com precisão de quarta ordem, que são computacionalmente
eficientes e estáveis e tem como atrativo alto grau de precisão além de uma molécula computacional pequena (ver Figura 1.b). Esses métodos são conhecidos como compactos (ver Gupta [1],
Mancera & Hunt [3, 4], Pandit et al. [6, 7] e Spotz & Carey [8]).
Nr
r
r
r
O r
r
r
r
r
r
r
r
r
r
sr
r
r
r
r
r
r
r
r
r
r
r rL
r
r
r
r
r
r
r
r
r
r
r
r
r
r
rs
r
r
r
r
r
r
r
r
r
r
r
r
S
a)
b)
Figura 1: Molécula computacional com: a) 29 pontos; b) 25 pontos.
Em diversos trabalhos são desenvolvidos métodos compactos 3×3 para as equações de Navier-
585
Stokes na formulação função corrente-vorticidade (ver Gupta [1]). Contudo, neste trabalho é
apresentado um método compacto e um não compacto, de quarta ordem 5 × 5, na formulação
função corrente. A escolha por essa formulação tem como atrativo a atribuição de condições de
fronteira apenas para a função corrente, sem o uso da vorticidade.
2
Objetivos
Os objetivos deste trabalho são calcular os erros RMS (Root Mean Square) e máximo para a
combinação de métodos de segunda e quarta ordem usando duas diferentes malhas, uma com o
dobro de pontos da outra e compará-los com os erros do método de diferença central de segunda
ordem.
3
Metodologia
No decorrer do trabalho são apresentados os procedimentos para a obtenção de um método
compacto e um não compacto, ambos de quarta ordem, para as equações de Navier-Stokes estacionárias, bidimensionais, na formulação função corrente. Inicialmente descreve-se a motivação
biológica do problema a ser estudado.
3.1
Motivação Biológica
Durante décadas muitos pesquisadores descreveram a aterosclerose como sendo o acúmulo
gradativo de gordura nas paredes das artérias, entretanto, sabe-se que a aterosclerose é um processo inflamatório crônico, identificado como uma doença na parede interna da artéria (ı́ntima)
e um processo alimentado pela inflamação (ver McKay et al. [5]). Com o passar do tempo,
a inflamação na parede da artéria começa a impedir a passagem do fluxo sanguı́neo levando a
obstrução.
O estudo de métodos numéricos para um canal com uma oclusão retangular (ver Figura
2), tem como uma de suas motivações o problema da obstrução de artérias, entretanto, devido a dificuldade em se trabalhar com um canal em 3D, simplificou-se o modelo utilizando
uma representação bidimensional para um problema tipo estenose. De tal situação decorre o
desenvolvimento e estudo do método numérico.
-
Figura 2: Canal com uma oclusão.
3.2
Equações de Navier-Stokes
As equações de Navier-Stokes que modelam escoamentos de fluidos incompressı́veis podem
ser escritas em variáveis não primitivas (vorticidade ζ e função corrente ψ), dadas por
ζ=
∂v
∂u
−
,
∂x ∂y
u=
∂ψ
,
∂y
v=−
∂ψ
,
∂x
(1)
em que u e v são as componentes do vetor velocidade nas direções x e y, respectivamente.
Substituindo as velocidades u e v na expressão da vorticidade ζ, tem-se
∂2ψ ∂2ψ
+
= −ζ.
∂x2
∂y 2
(2)
586
A formulação função corrente-vorticidade das equações de Navier-Stokes é dada pela expressão (2) e pela equação para fluidos estacionários expressa por
∂2ζ
∂ψ ∂ζ
∂ψ ∂ζ
∂2ζ
+
= Re
,
−
∂x2 ∂y 2
∂y ∂x
∂x ∂y
(3)
em que Re é o número de Reynolds.
Substituindo a equação (2) em (3) tem-se
∂4ψ
∂4ψ
∂4ψ
+
2
+
∂x4
∂x2 ∂y 2
∂y 4
∂3ψ
∂3ψ
+
∂x3
∂x∂y 2
∂ψ
= Re
∂y
!
∂ψ
−
∂x
∂3ψ
∂3ψ
+
∂x2 ∂y
∂y 3
!!
,
(4)
a qual é dita a formulação função corrente das equações de Navier-Stokes.
3.3
Obtenção dos Métodos de Quarta Ordem
Um método não compacto de quarta ordem para as equações de Navier-Stokes estacionárias é
obtido aproximando todas as derivadas da equação (4) por diferenças centrais de quarta ordem,
o que resulta em uma molécula computacional de 29 pontos. A molécula computacional é
mostrada na Figura 1.a, tendo os seguintes pontos cardeais
N=
Re DX
−1
Re DX
−1
−
, S=
+
,
4
3
4
6h
8h
6h
8 h3
−1
Re DY
L=
+
,
4
6h
8 h3
(5)
−1
Re DY
O=
−
,
4
6h
8 h3
em que h, DX e DY , são respectivamente, o espaçamento da malha e aproximações de quarta
∂ψ
ordem para ∂ψ
∂x e ∂y .
Quando o método não compacto é aplicado próximo a uma fronteira rı́gida, faz-se necessário
uma condição de fronteira extra para obter a solução procurada (ver Figura 3).
r
r
r r
r
r
r
r
r
r
r
r
r
r
r
r
r
r
r
r
r
r
r
r
r fronteira rı́gida
r r
r
r
Figura 3: Molécula próxima de uma fronteira sólida.
Assume-se que o fluido próximo de fronteiras rı́gidas tem a mesma velocidade da parede, ou
seja, as velocidades u e v são nulas. Da aplicação da condição de contorno de não deslizamento
tem-se duas condições de fronteira para a função corrente, ψ = 1 e ∂ψ
∂η = 0 com η = x ou y.
Entretanto, são requeridas três condições.
Para as fronteiras rı́gidas horizontais, tem-se que a condição ∂ψ
∂x = 0 e então para o ajuste
do método é necessário substituı́-la na equação (4), logo
∂4ψ
= 0,
∂y 4
(6)
fornece a condição de fronteira extra utilizada no método não compacto de quarta ordem para o
problema do canal com uma oclusão local. Condição similar é obtida para as fronteiras rı́gidas
verticais.
Após ilustrar a construção de um método não compacto de quarta ordem, apresenta-se os
procedimentos para a obtenção de um método compacto de quarta ordem para a equação (4).
Observa-se que o método compacto será obtido com as exclusões dos pontos cardeais. Seguindo
Wittkopf [9]
587
1. As derivadas parciais de primeira e segunda ordens de (4) são calculadas em relação a x
e a y. Estas expressões são denotadas por N Sx, N Sxx, N Sy e N Syy. Em seguida, para
estas expressões, todas as derivadas de ordem superior a dois são aproximadas usando
diferença central de segunda ordem, e as outras por diferença central de quarta ordem.
Estas aproximações são denotadas por [N Sx], [N Sxx], [N Sy] e [N Syy].
2. Na equação (4) representada por N S4 todas as derivadas são aproximadas por diferenças
centrais de quarta ordem. Esta aproximação é representada por [N S4 ] e resulta numa
molécula como já exibida na Figura 1.b, com os pontos cardeais dados em (5).
A equação
2
[N S4] + h
DX[N Sy] − DY [N Sx]
[N Sxx] + [N Syy]
+ Re
6
12
= 0,
(7)
elimina os pontos cardeais resultando em um método compacto 5 × 5 (ver Figura 1.a) para a
equação (4), com precisão O(Re2 h4 ).
3.4
Geometria e Resolução do Problema
Considera-se um canal com duas quinas reentrantes, tendo paredes sólidas em y = ±1 para
1
1
x < m e x > n, y = ± para m < x < n e ≤ |y| ≤ 1 para x = m e x = n, em que m e n são
2
2
inteiros, m < n (ver Figura 4). Devido a simetria na geometria o problema é resolvido apenas
para y ≥ 0 (ver Mancera & Hunt [4] e Pandit et al. [6]).
ψ=1
ψ=1
3y y 3
−
ψ=
2
2
ψ=1?
ψ=
3y y 3
−
2
2
ψ=0
Figura 4: Canal com uma oclusão.
As condições de fronteira para o problema são dadas por
=
1,
∂ψ
∂y
= 0 em y = 1, x ≤ m, x ≥ n, e y = 21 , m < x < n,
=
1,
∂ψ
∂x
= 0 em x = m e x = n, 21 ≤ y ≤ 1,
=
0,
∂2ψ
∂y 2
→
3
2y
− 21 y 3 ,
ζ → 3y,
quando x → −∞,
ψ →
3
2y
− 21 y 3 ,
ζ → 3y,
quando x → +∞.


ψ









ψ






ψ









ψ







= 0 em y = 0,
(8)
com fluxo de Poiseuille tanto a montante (sentido contrário a corrente), quanto a jusante
(seguindo o sentido da corrente).
Conforme verificado na Figura 3 para aplicação do método não compacto nos pontos interiores ao domı́nio computacional são necessários dois pontos fantasmas (fora do domı́nio computacional) para cada fronteira. No intuito de simplificar a quantidade de linhas de pontos
fantasmas, considera-se uma molécula computacional deslocada duas unidades de uma fronteira
sólida (linha preta), ilustrada na Figura 5.
588
r
r
rr
r
r
r
r
r
r
r
r
rr
rr
tr
rr
rr
r
r
r
rr
r
r
Figura 5: Mólecula próxima a fronteira sólida.
Os procedimentos adotados para implementar a combinação de métodos de segunda e quarta
ordem em uma geometria como a exibida na Figura 5 são dados por
• Aplica-se na primeira linha de pontos internos ao domı́nio computacional o método de
diferença central de segunda ordem (linha ). Por exemplo,
1
∂ψ
=
(ψi+1,j − ψi−1,j ) + O(h2 ).
∂x
2h
• Aplica-se o método de diferença central de quarta ordem (compacto/não compacto) nas
demais linhas internas do domı́nio computacional.
• Para as fronteiras rı́gidas (linha sólida) tem-se as seguintes condições
ψ = 1,
∂ψ
=0
∂x
ou
∂ψ
= 0.
∂y
(9)
• Para o eixo de simetria tem-se como condição ψ = 0. Verifica-se que para essa condição
os pontos são simétricos, então, ψi,j−1 = ψi,j+1 .
• Para o comportamento do fluido tanto na entrada quanto na saı́da do canal, assume-se a
seguinte condição
ψ=
1
3
y − y3.
2
2
(10)
• Na linha de pontos fantasmas aplica-se diferenças finitas ascendentes/descendentes de
). Por exemplo,
quarta ordem ou diferença central de segunda ordem (linha
∂ψ
1
=
(3ψi,j + 10ψi,j−1 − 18ψi,j−2 + 6ψi,j−3 − ψi,j−4 ) + O(h4 ).
∂y
12h
4
Resultados Numéricos
As simulações envolveram dois métodos numéricos de quarta ordem (compacto e não compacto) e o método de diferença central de segunda ordem para a formulação função corrente (ver
equação (4)). Obteve-se os resultados para duas malhas com diferentes tamanhos, sendo uma
com 256 × 64 e a outra com 512 × 128.
Para o problema do canal com um ressalto, tem-se uma malha uniforme, com o comprimento
do canal sendo quatro, e espaçamentos hx = hy = 1/ny em ambas as direções. A geometria
considerada para as simulações foi
• antes da oclusão: 0 ≤ x ≤ 1 e 0 ≤ y ≤ 1;
• oclusão: 1 ≤ x ≤ 2 e 0 ≤ y ≤ 0.5;
589
• depois da oclusão: 2 ≤ x ≤ 4 e 0 ≤ y ≤ 1.
Nas Tabelas 1 e 2 a , são apresentados os erros RMS e máximo para os seguintes métodos na
formulação função corrente: diferença central de segunda ordem (SO), método não compacto e
aproximações de quarta ordem na linha de pontos fantasmas (MNC 4h), método não compacto
e aproximações de segunda ordem na linha de pontos fantasmas (MNC 2h), método compacto
e aproximações de segunda ordem na linha de pontos fantasmas (MC 2h) e método compacto e
aproximações de quarta ordem na linha de pontos fantasmas (MC 4h).
Re
0
10
50
100
250
500
SO
1.15(-4)
1.20(-4)
1.69(-4)
2.06(-4)
2.93(-4)
4.28(-4)
Erros ψ - Métodos
MNC 4h MNC 2h MC 2h
2.63(-5)
2.49(-5) 2.48(-5)
2.72(-5)
2.57(-5) 2.42(-5)
3.25(-5)
3.11(-5) 4.46(-5)
3.44(-5)
3.13(-5) 7.03(-5)
3.75(-5)
3.08(-5)
**
5.41(-5)
3.91(-5)
**
MC 4h
2.62(-5)
2.51(-5)
4.58(-5)
7.47(-5)
**
**
Tabela 1: Erros RMS
Re
0
10
50
100
250
500
SO
4.95(-4)
7.14(-4)
1.13(-3)
1.47(-3)
1.84(-3)
2.11(-3)
Erros ψ - Métodos
MNC 4h MNC 2h MC 2h
1.08(-4)
1.06(-4) 1.05(-4)
1.47(-4)
1.49(-4) 1.51(-4)
2.00(-4)
2.06(-4) 2.20(-4)
2.17(-4)
2.23(-4) 2.65(-4)
2.45(-4)
2.24(-4)
**
3.63(-4)
3.02(-4)
**
MC 4h
1.08(-4)
1.49(-4)
2.20(-4)
2.66(-4)
**
**
Tabela 2: Erros máximo
Como pode ser visto nas tabelas, os erros RMS variam de ∼ 10−5 a ∼ 10−4 e os erros
máximo de ∼ 10−4 a ∼ 10−3 , respectivamente, para 0 ≤ Re ≤ 500. O sı́mbolo ** indica que
o método de Newton não convergiu para a dada precisão estipulada depois de 10 iterações.
Também é observado que os erros RMS para as combinações dos métodos mistos MNC 2h e
MC 2h apresentaram erros menores quando comparados com os métodos mistos MNC 4h e MC
4h. Nota-se que para todos os erros a ordem de grandeza é a mesma, porém a precisão de
alta ordem, que era esperada para os métodos nominados MNC 4h e MC 4h, não foi obtida.
Esse comportamento deve-se as duas quinas reentrantes (no afunilamento e na expansão do
canal), aos pontos singulares e ao uso de malha uniforme. O fato é reforçado com os resultados
apresentados por Hunt [2] para o problema do degrau fazendo uso de malha não uniforme.
Para o problema estudado, verificou-se que os métodos mistos MNC apresentam erros menores
quando comparados com os erros dos métodos mistos compacto MC, e isto pode ser explicado
pelo fato que o método não compacto tem erro de truncamento de ordem O(Re h4 ), enquanto
os erros de truncamento para o método compacto são de ordem O(Re2 h4 ).
Na Figura 6 são apresentadas as linhas de corrente e para a construção dessas linhas foi
realizado os cálculos em uma malha 512 × 128. Nas Figuras (6(a))-(6(b)) são exibidas as linhas
de corrente para o método MNC 2h, figuras essas semelhantes para os métodos MNC 4h, MNC
2h, MC 4h e MC 2h. Para Re = 500, os métodos mistos MC 4h e MC 2h não convergiram. A
Figura 6(c) apresenta linhas de corrente para o método MNC 4h.
Nota-se na Figura 6(a) a simetria (antes e após a oclusão) que só ocorre em Re = 0; na
Figura 6(b) tem-se uma pequena recirculação na oclusão e após a oclusão a recirculação é bem
a
Em que a notação p(−q) é equivalente a p × 10−q .
590
acentuada, e na Figura 6(c) há duas regiões de recirculação, sendo respectivamente, uma na
oclusão e outra mais intensa após a oclusão.
0.8
0.8
0.8
0.6
0.6
0.6
y
1
y
1
y
1
0.4
0.4
0.4
0.2
0.2
0.2
0
0
1
2
3
4
0
0
1
2
3
4
0
0
1
2
x
x
x
(a) Re = 0
(b) Re = 250
(c) Re = 500
3
4
Figura 6: Linhas de corrente.
5
Conclusão
Apresentou-se resultados para o problema do canal com um ressalto utilizando a combinação
de métodos de quarta ordem com os métodos de segunda ordem. Escolheu-se essa combinação
com o intuito de diminuir as dificuldades computacionais. Conclui-se que essa combinação além
da facilidade de implementação do método para a geometria considera, ainda produz resultados
mais precisos quando comparados com o método de segunda ordem.
Referências
[1] Gupta, M. M., High accuracy solutions of incompressible Navier-Stokes equations, J. Comput. Phys., 93 (1991) 343-359.
[2] Hunt, R., The numerical solution of the laminar flow in a constricted channel at moderately
high Reynolds number using Newton iteration, Int. J. for Numer. Meth Fl., 11 (1990) 247259.
[3] Mancera, P. F. A. and Hunt, R., Some experiment with high order compact methods using
a computer algebra software - Part I, Appl. Math. Comput., 174 (2006) 775-794.
[4] Mancera, P. F. A. and Hunt, R., Fourth-order method for solving the Navier-Stokes equations in a constricting channel, Int. J. Numer. Meth. Fluids, 25 (1997) 1119-1135.
[5] McKay, C., McKee S., Mottram, N., Mulholland, T. and Wilson, S., Towards a model of
atherosclerosis, University of Strathclyde, Glasgow, Scotland, (2005) 1-29.
[6] Pandit, S. K., Kalita, J. C. and Dalal, D. C., A transient higher order compact scheme
for incompressible viscous flows on geometries beyond rectangular, J. Comput. Phys., 126
(2007) 1110-1124.
[7] Pandit, S. K., Kalita, J. C. and Dalal, D. C., A fourth-order accurate compact scheme
for the solution of steady Navier-Stokes equations on non-uniform grids, Computers and
Fluids, 37 (2008) 121-134.
[8] Spotz, W. F. and Carey, G. F., High-order compact scheme for the steady stream-function
vorticity equations, Int.J. Numer. Meth. Engng., 38 (1995) 3497-3512.
[9] Wittkopf. A. “High order wide and compact schemes for the steady incompressible NavierStokes equations”, Master’s thesis, Simon Fraser University, 1994.
591