Aula 10
esquemas numéricos para a resolução dos
sistemas de equações de conservação
métodos de resolução
• diferenças finitas e volumes finitos
esquemas explícitos de 1ª ordem no tempo para a
resolução de
t  U   x  F  G
diferenças finitas, exemplo: explícito, 1ª ordem no espaço e no tempo)
Uin1

Uin


t n

Fi 1  Fin  tG in
x
n +1
n
i 1
i
i+1
funções U e F, e respectivas derivadas, discretizadas e transformadas em
funções de malha.
derivadas substituídas por acréscimos finitos: diferenças finitas centradas,
progressivas, regressivas ou combinações entre estas.
métodos de resolução
• diferenças finitas e volumes finitos
esquemas explícitos de 1ª ordem no tempo para a
resolução de
t  U   x  F  G
diferenças finitas, fórmula geral:




1 Uin1; t , x  0 Uin ; t , x  0
bibliografia: Richtmyer e Morton 1967, Cunge et al. 1980, Hirsh 1989,
1990.
métodos de resolução
• diferenças finitas e volumes finitos
esquemas explícitos de 1ª ordem no tempo para a
t  U   x  F  G
resolução de
volumes finitos
Uin1  Uin 


t *n
Fi  1  Fi*n1  tG in
2
2
x
fluxos avaliados nas fronteiras entre volumes centrados em i1, i, i+1...
*
Fi
1
*
Fi
1
2
2
U i 1
2
2
i3/2
i1
i
i1/2
i+1
i+3/2
i+1/2

Fi*n1  F Fin2 , Fin1 , Fin , Fin1,...

métodos de resolução
• diferenças finitas e volumes finitos
esquemas explícitos de 1ª ordem no tempo para a
resolução de
t  U   x  F  G
volumes finitos, em geral
Uin1  Uin 



t *n
Fi  1  Fi*n1  tG*i G in1 , G in , G in1 , G in11 , G in 1 , G in11
2
2
x

bibliografia: Leveque 1980, Hirsh 1990, Toro 1998, Leveque 2002.
métodos de resolução
• diferenças finitas e volumes finitos
conceitos fundamentais para a caracterização dos esquemas
• ordem (ver acetatos)
exemplo: discretização por diferenças progressivas
 xi1  xi 2
3
2
fi1  fi   x  f   xi 1  xi    x  f 
 O  xi 1  xi 
ì
ì
2!
f f
 x  f   i 1 i  O  x 
ì
x


erro de truncatura
exemplo: discretização por diferenças centradas
fi1  fi   x  f 
 x  xi    x  f  ì
ì i 1
fi1  fi   x  f 
xi  xi 1    2 x  f 

ì
ì
2
 
fi 1  fi 1  2  x  f   O x3
ì
 xi1  xi 2
2!
 xi  xi1 2


 O  x  x  
 O  xi 1  xi 
3
i
2!
x  f  
ì
i 1
3
 
fi 1  fi 1
 O x 2
2x
“segunda ordem”: erro
de truncatura
proporcional a x2
métodos de resolução
• diferenças finitas e volumes finitos
conceitos fundamentais para a caracterização dos esquemas
• consistência (ver acetatos)
exemplo: esquema upwind não conservativo aplicado a t    x   g
in1  in  i
como
t n
(i  in1 )  tgi
x
 
in1  in  t  i t  O t 2
 
in   t  i t  O t 2  in  i
e
 
in  in1   x  ì x  O x2
 
t n
(i 1   x  ì x  O x 2  in1 )  tgi
x
t  i  O  t   i  x  ì  gi  O  x 
e, quando x → 0
t    x   g
ou seja, obtém-se a equação inicial.
logo, é consistente.
métodos de resolução
• diferenças finitas e volumes finitos
conceitos fundamentais para a caracterização dos esquemas
• estabilidade (ver acetatos)
exemplo: esquema upwind não conservativo aplicado a t    x   0
in1  in  i

t n
i 1  in
x

verificar a frequência angular e o factor de ampliação associados às perturbações
propagadas no domínio de cálculo
sendo
I ( n 1) t ix 
in 1  ˆ e 
I ( n 1) t  (i 1) x 
n 1  ˆ e 
e
ˆ e I  nt ix 
in  
I nt  (i 1) x 
n  ˆ e 
i 1
i 1
obtém-se

ˆ eI (n1)t 
ˆ eI ix  
ˆ eI nt 
ˆ eI ix  Cr 
ˆ eI nt 
ˆ eI (i1)x  
ˆ eI nt 
ˆ eI ix

t
em que
Cr   i
x

métodos de resolução
• diferenças finitas e volumes finitos
• estabilidade
exemplo: esquema in 1  in   i
ˆ eI nt 
ˆ eI ix obtém-se
dividindo por 
eI t
 1  Cr

t n
i 1  in
x

eIx  1
a frequência angular numérica poderá ter componente imaginária (ao
contrário das frequências angulares do sistema hiperbólico)
  R  I I
assim
eI  t e t
R
I
 1  Cr
eIx  1
 cos  R t   I sin  R t   e t  1  Cr   cos  x   I sin  x   1
I
dividindo a equação nas suas partes real e imaginária
 cos   t  eI t  1  Cr  cos  x   1
R


 sin  R t  eI t  Cr sin  x 

métodos de resolução
• diferenças finitas e volumes finitos
• estabilidade
exemplo: esquema in 1  in   i
combinando as equações



t n
i 1  in
x

 1  Cr 2 sin 2 x e2I t e2I t  1  Cr cos x  1
 
   



Cr sin  x 
 tan   t  
R

1  Cr 1  cos  x  
2
 I t
e

1

Cr
cos

x

1
 Cr 2 sin 2  x 







Cr sin  x 
 tan  R t  

1  Cr 1  cos  x  



2

métodos de resolução
• diferenças finitas e volumes finitos
• estabilidade
exemplo: esquema in 1  in   i

t n
i 1  in
x

 eI t  1  2Cr  Cr  1  cos  x   1



Cr sin  x 
tan


t




R
1  Cr 1  cos  x  

t
x
se Cr  i
 1  i 
 i
x
t

I t
1
não há ampliação/atenuação
 e


  t  x  R  1
não há dispersão, celeridade numérica não
R


depende do número de onda
assim, se
Cr  1
o esquema reproduz a forma de propagação hiperbólica:
sem dissipação e sem dispersão.
métodos de resolução
• diferenças finitas e volumes finitos
exemplo: esquema in 1  in   i
• estabilidade

t n
i 1  in
x

para quaisquer outros valores de Cr é necessário traçar os retratos de amplitude...
e t
I

 2x  
 1  2Cr  Cr  1  cos 
  1
L

 

1.4
L : comprimento de onda das
factor ampliação (-)
1.3
perturbações para quaisquer
outros valores de Cr é necessário
traçar os retratos de amplitude
1.2
Cr = 1.1
1.1
Cr = 1.0
Cr = 0.9
1
0.9
0.8
Cr = 0.5
0.7
0.6
0.5
0.4
1
10
100
L /x (-)
1000
métodos de resolução
• diferenças finitas e volumes finitos
• estabilidade
exemplo: esquema in 1  in   i
... e de fase...
celeridade numérica adimensional (-)


 x 
Cr sin  2 


R
L
L




atan 


  2x
 x   
 1  Cr 1  cos  2   
L  



1.4
1.3
1.2
Cr = 1.1
1.1
Cr = 1.0
1
0.9
Cr = 0.9
0.8
0.7
0.6
Cr = 0.5
0.5
0.4
1
10
100
L /x (-)
1000

t n
i 1  in
x

métodos de resolução
• diferenças finitas e volumes finitos
conceitos fundamentais para a caracterização dos esquemas
• estabilidade (ver acetatos)
condição de Courant–Friedrichs–Lewis
Cr  
t
x
1   
 i
x
t
n +1


n
i 1
i
i+1
métodos de resolução
• diferenças finitas e volumes finitos
conceitos fundamentais para a caracterização dos esquemas
• convergência (ver acetatos)
teorema de equivalência de Lax: convergente sse é consistente, estável e bem
condicionado
métodos de resolução
• esquemas de discretização por diferenças finitas
• tipo upwind
upwind não conservativo
1990’s, actual com complemento de 2ª
ordem
• tipo Lax-Wendroff
Lax-Wendroff 2ª ordem, MacCormack
finais 1980’s, 1990’s com TVD, actual...
• tipo “box”
Preissmann
anos 70, 80 e início 90’s
métodos de resolução
• esquema upwind
in1
com

in
 i 1
2
t n
n
 t
(i  i 1 )   i  1
(in1  in )  tgi
2 x
x
  1 (   )
2
  1 (   )
2
0
0
n +1
n +1
n
n
i1
t    x   g
(1ª ordem, explícito) para
i
i 1
i+1
estabilidade condicionada a
Cr  
i
i+1
t
t


 1  notar que Cr  
 0
x
x


métodos de resolução
• esquema Lax-Wendroff
(2ª ordem, explícito) para
t    x   g
in1
2
2



t
n
n
n
n
n
n
i
i  t 
 i 
(i 1  i 1 ) 
(


2



i 1
i
i 1 )  tgi
 
2 x
2  x 
esquema instável
n +1
n
i1
i
i+1
necessita viscosidade artificial ou correcção TVD (total variation diminishing)
estabilidade condicionada a
Cr  
t
1
x
métodos de resolução
• esquema tipo box (Preissmann, implícitio)
t    x   g
para









 n1
1   n 1
i 1  in1 
i  in
t
t

 n1
1  n
i 1  in1 
i 1  in
x
x
 in11  in1  (1  ) in 1  in
t () 
t
 in11  in 1  (1  ) in1  in
 x () 
x








ni 1  ni1  (1  )ni
2
  n1  
  n1
1 
n
n






1
1

 i

 i 1 
i 2
i 2

t

x

t

x





(1  ) n
1  n
 in1 
i  in 1
i 1  in  gin 1
2 x
2
t
t


métodos de resolução
• esquema tipo box (Preissmann, implícitio)
t    x   g
para
  n1  
  n1
1 
n
n






1
1

 i

 i 1 
i 2
i 2
x 
x 
 t
 t

(1  ) n
1  n
 in1 
i  in 1
i 1  in  gin 1
2 x
2
t
t


n +1

valores usuais:

  12
  0.7
n
i
i+1
estabilidade condicionada a
 1


2   1  0
2
Cr
métodos de resolução
• notas:
- o esquema upwind é muito difusivo porque o erro de
truncatura é da ordem de x2;
- os esquemas de Lax-Wendroff são muito dispersivos e,
sendo de 2ª ordem, são oscilatórios; necessitam sempre
correcções TVD ou de viscosidade artificial;
- os esquemas implícitos do tipo Preissmann são dispersivos
e difusivos (a disfusão esconde a dispersão); sendo implícitos
apresentam bom desempenho no que diz respeito ao tempo
de cálculo (t pode ser independente de );
métodos de resolução
• esquema upwind
Uin1

Uin

Λ i 1
2
(1ª ordem, explícito) para t  U  M x U  G
t n
n
 t
(Ui  Ui 1 )  Λ i  1
(Uin1  Uin )  tG i
2 x
x
com
 
Λ i 1
2
 (1) 
  i 1
2




0

0




 

1
i 2 

 
(k )
 
1(k )1
2

 
1(k )1
2
ver: diagonalização da matriz Λ = S1A1BS



 12  1(k )1  1(k )1 
2 
 2


 12  1(k )1  1(k )1 
2 
 2
métodos de resolução
• esquema upwind
Uin1

Uin

Λ i 1
2
(1ª ordem, explícito) para t  U  M x U  G
t n
n
 t
(Ui  Ui 1 )  Λ i  1
(Uin1  Uin )  tG i
2 x
x
( k )  0
( k )  0
n +1
n +1
n
n
i1
i
i 1
i+1
estabilidade condicionada a
Cr  max i( k )
k ,i
t
1
x
i
i+1
métodos de resolução
• esquema MacCormack
para t  U  M x U  G




t n
Fi  Fin1  tG in
x
t p
Uic  Uin 
Fi 1  Fip  tG ip
x
1
Uin1  Uip  Uic
2
n +1
Uip  Uin 

(2ª ordem, explícito)
nota: o desempenho do esquema é
melhorado se se alternar a aplicação das
diferenças progressivas e regressivas

p
n
estabilidade condicionada a
Cr  max
k ,i
i( k )
t
1
x
i1
i
i+1
métodos de resolução
• esquemas de discretização por volumes finitos
• tipo Godunov
HLLC (com Riemann solvers)
final década 1990, actualmente
• quasi-2ª ordem
Ying (2000)
proposto em 2000
métodos de resolução
• esquema HLLC
(Harten-Lax-van Leer + Contact wave, 1ª
ordem, explícito) para t  U   x  F   G


t *n
F i  1  F*in 1  tG in
2
2
x
Uin1  Uin 
stencil de Godunov
problema de Riemann
local
t
L
valores médios nos
volumes de cálculo
R
U*L
UL
i
5
2
i
3
2
i  12 i  12 i 
3
2
x
t
SL
S*
SR
U*R
UR
x
estrutura da solução para o
problema de Riemann local
métodos de resolução
• esquema HLLC
(Harten-Lax-van Leer + Contact wave, 1ª
ordem, explícito) para t  U   x  F   G
t
SL
S*
U*L
SR
U*R
UL
Riemann solvers para o esquema HLLC (podem
ser aproximados)
SL  uin 1  cin 1
2
S* 
UR
x
u*in 1 
2
h*in 1
2


1 n
1
ui 1  uin 
2
2


 
uin 1
2

9.8hin 
1
1
 hin1  hin  uin1  uin
2
4


SR  uin 1  cin 1
2
2
u*in 1
2
hin uin

hin1 uin1
hin  hin1
hin1  hin 

9.8hin1 
 hin1  hin 
 hin1  hin 
ghin1  ghin

cin 1 
2
2
   
2

9.8 12  hin1  hin

hin  hin1
nota: não há
descontinuidade em u e h
através da onda de
contacto
2



métodos de resolução
• esquema HLLC
(Harten-Lax-van Leer + Contact wave, 1ª
ordem, explícito) para t  U   x  F   G
t
SL
S*
U*L
SR
U*R
UL
expressão alternativa para as velocidades das
ondas
SL  uin  cin  Li  1
UR
2
S*  u*in 1
2
x
 Li  1
2
R i 1
2
1
se h*ni  1  hin
2


n 2
 1 hin h*ni  1  hin
hi
se h*ni  1  hin
2
2
2

 1
se h*ni  1  hin1
2


n 2
n
n
 1 hin1 h*ni  1  hin1
h
se
h
1  hi 1
i

1
*
i

2
2
 2
 


 
SR  uin  cin  R i  1
2
métodos de resolução
• esquema HLLC
(Harten-Lax-van Leer + Contact wave, 1ª
ordem, explícito) para t  U   x  F   G
t
SL
U*L
UL
S*
SR
U*R
Yb*L in

1
2

1 n
Ybi 1  Ybin
2
Yb*R in
UR
x
C*L in
1
2
 Cin
fluxos:
F*i  1
2
if 0  SL
 FL
 F F S U U
 *L
L
L  *L
L  if S L  0  S*

 F*R  FR  SR  U*R  U R  if S*  0  SR

if 0  SR
FR

1
2


1 n
Ybi 1  Ybin
2
 C*R in
1
2
 Cin1

métodos de resolução
• esquema HLLC
(Harten-Lax-van Leer + Contact wave, 1ª
ordem, explícito) para t  U   x  F   G
t
SL
U*L
UL
fluxos, expressões alternativas:
S*
SR
U*R
k = 3:
 
UR
x
F*i  1
2
k = 1, 2:
F*i  1
2
 Fin


 Fin SR  SL Fin1  SR SL Uin1  Uin

n

 Fi 1


n
 n
F
C
 1i *L i  1
2




n
n
 F1i 1 C*R 1
i 2


 
S*  0
S*  0
if 0  SL
 SR  SL 
if SL  0 and SR  0
if SL  0 and SR  0
métodos de resolução
• esquema de Ying (quasi 2ª ordem, quasi-implícito)
para
t  U   x F  S
Uin1  Uin 


t *n
F i  1  F*in 1  tH in 1  tG in
2
2
x
 uh 


nota: o vector de fluxo contém apenas os termos de fluxo físicos, eg. F  u 2 h 
 Cuh 


o vector H contém o declive da superfície livre, i.e. a soma do declive do fundo e
da altura do escoamento
métodos de resolução
esquemas de diferenças e volumes finitos, notas:
- os diferenças finitas upwind e volumes finitos Godunov simples são
equivalentes; são ambos altamente difusivos, não dispersivos e possuem
grande aptência para capturar choques;
- o esquema MacCormack é do tipo Lax-Wendroff; é de 2ª ordem no espaço e,
logo, oscilatório; só se obtém soluções de qualidade se se controlar as
oscilações com viscosidade artificial ou com algoritmos TVD com limitadores
de fluxo;
métodos de resolução
esquemas de diferenças e volumes finitos, notas:
- o esquema HLLC é adequado para escoamentos que desenvolvem choques
fortes; a sua aplicação a problemas quasi-estacionários tem vindo a ser
tentada com resultados encorajadores
- o esquema de Ying parece simular satisfatoriamente escoamentos
transitórios e escoamentos permanentes; todavia é bastante difusivo e pode
gerar oscilações na presença de choques fortes.
métodos de resolução
condições de fronteira
(ver acetatos)
- método das características (ver acetatos)
- difícil de aplicar em problemas com leito móvel;
- implementável com trechos fictícios a montante e jusante;
- células “fantasma”
- resolução das equações nos nós de fronteira (requer nós fictícios a
montante e a jusante);
- fácil de implementar; matematicamente, pode representar um
problema mal condicionado;
métodos de resolução
condições de fronteira
(ver acetatos)
- modelos em equilíbrio, condição de fronteira para as equações de
conservação da massa de sedimentos
- não desejável:
Yb1n1  Yb1n 


t
B
qs (t )  qs1n
x
1 p
- preferível, integrar na primeira célula de cálculo:
x

0
t
t t
(1  p)t Yb   x qs  dtdx  0