Universidade da
Beira Interior
AULA 11
MODELAÇÃO BIDIMENSIONAL DA
PROPAGAÇÃO DE ONDAS DE CHEIA COM
FRENTE ABRUPTA
João Leal (UBI)
IST, Lisboa
2006
Universidade da
Beira Interior
MODELO CONCEPTUAL
(morfodinâmico)
h
hw
hc
uw
clear water
uc
sheet-flow
zb
immobile bed
Depth average
theory

C  Cc hc uc ( x )  uc ( y )
2
  w 1   s  1 C 
2

w h w
Cc
ub = 0
u  uc hc  uwhw  h
C w= 0
C b= 1 - p
h u( x )  u ( y )
2
2

w h w + c h c
NOTA:
 u x  

u
u y  
Universidade da
Beira Interior
EQUAÇÕES DE CONSERVAÇÃO
U F  U  G  U 
MORFODINÂMICO
t
Variáveis dependentes
 zs 
 u h 
( x) 
U
u( y ) h 


 ze 

x

S
y
Vectores de fluxo
hu( y )
hu( x )








wuw( x ) uw( y ) hw  c uc ( x ) uc ( y ) hc 
 wuw( x ) 2 hw  c uc ( x ) 2 hc   
F U  

 G U  
2
2

u
u
h


u
u
h
c c( x) c( y ) c 
 wuw( y ) hw  c uc ( y ) hc   
 w w( x ) w ( y ) w




Chu( y )
Chu( x )





  g w hw2  2w hw hc  c hc 2
Termos de fonte
0




  g   w hw  c hc  zb  bc ( x ) 


x
S

z
  g   w hw  c hc  b  bc ( y ) 
y




0



2
Variáveis dependentes
cota da sup. livre
caudal mássico por
unidade de largura
zs  h  zb
uh
ze  1  p  zb  Cc hc
cota do fundo equivalente aos
sedimentos acumulados na coluna de água
Universidade da
Beira Interior
ESQUEMA DE MACCORMACK
 O esquema (MacCormack, 1969) foi desenvolvido e implementado no
âmbito da dinâmica de gases. A sua facilidade de implementação
aliada a uma precisão de segunda ordem faz com que seja um dos
esquemas mais utilizados na modelação numérica de ondas de cheia.
 O esquema numérico de MacCormack é explícito e constitui a variante
de dois passos mais popular do esquema de Lax-Wendroff. Pertence à
classe dos métodos de passo fraccionado (“fractional-step”) e garante
uma aproximação de segunda ordem no tempo e no espaço.
Universidade da
Beira Interior
ESQUEMA DE MACCORMACK 1D
x  x
Diferenças Finitas
U F

0
t x
t  t
Uin1 

1 p
Ui  Uic
2
Previsão
 n t n
n
U

F

F
i
i

1
i

x
p
Ui  
U n  t F n  F n
i 1
 i x i

Correcção




 n t p
p
U

F

F
i 1
se n  2, 4,...
 i x i
c
Ui  
U n  t F p  F p
se n  3,5,...
 i x i 1 i


se n  2, 4,...


se n  3,5,...
Aplicação alternada de diferenças progressivas e regressivas,
respectivamente, nos passos de previsão e de correcção
Universidade da
Beira Interior
tempo
(t )
n+ 1
correcção
(progressivas)
n
previsão
(regressivas)
n 1
n
i 1
i
i +1
FRONT. DE JUSANTE
n+ 1
FRONT. DE MONTANTE
x
barragem
ESQUEMA DE MACCORMACK 1D
1
0
1
2
CONDIÇÕES INICIAIS
3
m c 1 m c m c +1
m x 1 m x
espaço
(x )
Universidade da
Beira Interior
ESQUEMA DE MACCORMACK 2D
x  x
Diferenças Finitas
t  t
y  y
U F G


0
t x y
Uin,j1 

1 p
Ui , j  Uic, j
2
Previsão
Uip, j
t
 n
U

i
,
j

x

 U n  t
 i , j x

 U n  t
 i , j x

t
 Uin, j 
x

F
n
i 1, j
F
n
i, j

 Fin1, j

F
 Fin, j

F
 Fin1, j

n
i 1, j
n
i, j
t
x
t

x
t

x
t

x
 Fin, j 

Correcção
G
n
i , j 1
 G in, j

se n  2, 6,...
G
n
i , j 1
 G in, j

se n  3, 7,...
G
n
i, j
 G in, j 1

se n  4,8,...
G
n
i, j
 G in, j 1

se n  5,9,...
Uic, j
t
 n
 Ui , j  x

 U n  t
 i , j x

 U n  t
 i , j x

t
Uin, j 
x

F
p
i, j

F
 Fi ,pj

F
 Fi p1, j

F
 Fi ,pj

p
i 1, j
p
i, j
t
x
t

x
t

x
t

x
 Fi p1, j 
p
i 1, j
G
p
i, j
 G ip, j 1

se n  2, 6,...
G
p
i, j
 G ip, j 1

se n  3, 7,...
G
p
i , j 1
 G ip, j

se n  4,8,...
G
p
i , j 1
 G ip, j

se n  5,9,...
Aplicação alternada de diferenças progressivas e regressivas,
respectivamente, nos passos de previsão e de correcção e também nos
fluxos segundo x (F) e segundo y (G)
Universidade da
Beira Interior
ESQUEMA DE MACCORMACK 2D
y
i ,j +1
p
y
i ,j +1
c
c
i +1,j
x
i ,j
i -1,j
p
i -1,j
c
c
p
i ,j -1
Tempo de cálculo n
Tempo de cálculo n +1
y
i ,j +1
y
i ,j +1
p
c
i +1,j
x
i ,j
i -1,j
c
i +1,j
x
i ,j
i ,j -1
p
p
c
i ,j -1
Tempo de cálculo n +2
c
i +1,j
x
i ,j
i -1,j
p
p
i ,j -1
Tempo de cálculo n +3
Universidade da
Beira Interior
CORRECÇÃO DE ALTA RESOLUÇÃO TVD
 De acordo com o teorema de Godunov, a utilização de um esquema
de ordem superior à primeira produz oscilações espúrias na
presença de descontinuidades.

a)
b)
SOLUÇÃO DO PROBLEMA DE RIEMANN
 3,0
3,0
2,5
2,5
2,0
2,0
1,5
1,5
1,0
-1,0
1,0
-1,0
-0,5
0,0
1,0
0,5
distância
esquema primeira ordem (Godunov)
simula mal (“adoça”) as descontinuidades
c)
-0,5
0,0
1,0
0,5
distância
esquema de segunda ordem (MacCormack)
Simula bem as descontinuidades mas
apresenta oscilações numéricas
Universidade da
Beira Interior
a)
b)
CORRECÇÃO DE ALTA
RESOLUÇÃO TVD

3,0

3,0
dos
2,5
 A tentativa de eliminar as2,5oscilações numéricas resultantes
termos de ordem superior
que
2,0 deu origem aos esquemas TVD
2,0
garantem que a variação1,5
total das sucessivas soluções
1,5 numéricas
n 1
n
não aumenta no tempo, ou
1,0 seja, TV  U   TV  U 
1,0
-1,0
-0,5
1,0
0,5
geramdistância
novos
0,0
 Esta propriedade garante que não se
-1,0
-0,5
0,0
mínimos ou
máximos locais e que os mínimos e máximos locais existentes não
são decrescentesb)ou crescentes, respectivamente.c)

1,0
0,5
distância
3,0

3,0
2,5
2,5
2,0
2,0
1,5
1,5
1,0
-1,0
-0,5
0,0
1,0
0,5
distância
esquema de MacCormack sem TVD
1,0
-1,0
-0,5
0,0
1,0
0,5
distância
esquema de MacCormack com TVD
1
0,5
distânc
Universidade da
Beira Interior
CORRECÇÃO DE ALTA RESOLUÇÃO TVD
 Porém, o teorema não garante a unicidade da solução fraca obtida,
podendo obter-se soluções sem significado físico (espúrias). Para
que tal não aconteça é necessário garantir a satisfação da condição
de entropia:
 k  Uesquerda   Sdes   k  U direita 
0,5
esquema de MacCormack sem condição de
entropia (choque não físico)
Y [m]
0,4
0,3
0,2
0,1
0
-8
-6
-4
-2
0
2
x [m]
4
6
8
10
12
Universidade da
Beira Interior
CORRECÇÃO DE ALTA RESOLUÇÃO TVD
 A condição de entropia desenvolvida por Harten e Hyman (1983)
escreve-se
1,3
i 1 2, j
 a1,3
 i 1 2, j

1,3
 i 1 2, j

se
1,3
ai1,3
1 2, j  i 1 2, j
se
1,3
ai1,3
1 2, j  i 1 2, j

1,3
i , j 1 2

1,3
1,3
1,3
1,3


1,3
i 1 2, j  max 0, ai 1 2, j  ai , j , ai 1, j  ai 1 2, j 
 a1,3
 i , j 1 2

1,3
 i , j 1 2

se
1,3
ai1,3
, j 1 2  i , j 1 2
se
1,3
ai1,3
, j 1 2  i , j 1 2


1,3
1,3
1,3
1,3


1,3
i , j 1 2  max 0, ai , j 1 2  ai , j , ai , j 1  ai , j 1 2 
Universidade da
Beira Interior
ESTABILIDADE NUMÉRICA
 O esquema de MacCormack, tal como todos os esquemas explícitos,
tem que verificar a condição de estabilidade de Courant-FriedrichsLewy.
 Esta condição impõe o tamanho dos passos de cálculo da seguinte
forma:
t  Cr 
xy
 max
x 2  y 2
Passo de tempo
Número de Courant
Valor máximo das características no
instante de tempo anterior
Universidade da
Beira Interior
ESCOLHA DAS VARIÁVEIS DEPENDENTES
 Na resolução numérica de sistemas de equações às derivadas parciais
em que a solução apresente descontinuidades (choques) é necessário
ter em atenção que as variáveis dependentes devem ser escolhidas de
forma ao sistema ser fisicamente conservativo.
 Note-se que para o sistema ser matematicamente conservativo bastará
U F  U  G  U 


S
t
x
y
Porém isso, por si só, não garante que o sistema é fisicamente
conservativo.
escreve-lo na forma conservativa
Universidade da
Beira Interior
Exemplo: Eqs. de Saint-Venant 1D para canais horizontais
prismáticos com secção rectangular e sem atrito
variáveis
primitivas
variáveis
conservativas
h  hu

0
t
x
h q

0
t x


q  q2 h  gh2 2

0
t
x


u  u2 2  gh

0
t
x
NOTA: ambos os sistemas são matematicamente conservativos
mas o segundo não é fisicamente conservativo.
U FU

0
t
x
Formulações conservativas VS. não-conservativas
Universidade da
Beira Interior
Admitindo que a solução dos sistemas contém um choque
Condições de Rankine-Hugoniot: são aplicáveis a soluções
descontínuas de sistemas de leis de conservação
hiperbólicos
F  Si U

Si
F(UR )  F(UL )  Si UR  UL 
com
UR
U FU

0
t
x
UL
x
Formulações conservativas VS. não-conservativas
Universidade da
Beira Interior
Condições de Rankine-Hugoniot aplicadas ao sistema
com variáveis conservativas
qR  qL


hR  hL 
q 2 h  gh 2 2  q 2 h  gh 2 2  Scons q  q 
 R

 R
R
R
L
L
L
L

Scons
ghL
 uR 
hR  hL 
2hR
Formulações conservativas VS. não-conservativas
Universidade da
Beira Interior
Condições de Rankine-Hugoniot aplicadas ao sistema
com variáveis primitivas
hRuR  hLuL


hR  hL 
u 2 2  gh  u 2 2  gh   Sncons u  u 
 R
 R
R
L
L
L

Sncons
hL2
 uR  2g
hR  hL
hR hL
Formulações conservativas VS. não-conservativas
Velocidade do
choque S
Universidade da
Beira Interior
facilmente se demonstra que Scons  Sncons
50
40
30
20
10
0
Scons
Sncons
Si
UR
UL
0
5
10
15
20
Força do choque hR/hL
x
CONCLUSÃO:
Soluções numéricas baseadas nas variáveis primitivas dão velocidades dos choques
erradas (menores do que a real), tanto mais erradas quanto maior for a força do
choque
Universidade da
Beira Interior
Exemplo:
0,5
Y [m]
0,4
0,3
0,2
0,1
0
-8
-6
-4
-2
8
6
4
2
x [m]
solução de Ritter
RUBATS - variáveis dependentes h e q
RUBATS - variáveis dependentes h e u
0
10
12
Universidade da
Beira Interior
MODELO NUMÉRICO
conservação da massa da
mistura
Conservação da quantidade
de movimento da mistura nas
direcções x e y
conservação da massa
de sedimentos
correção TVD
aproximações de Roe
condição de entropia de Harten e Hymen
Limitador de fluxo de Van Leer
viscosidade artificial de Jameson
Universidade da
Beira Interior
VISCOSIDADE ARTIFICIAL
 A introdução da viscosidade artificial no esquema de MacCormack
traduz-se, tal como na metodologia TVD, na introdução de um termo de
correcção que pode ser visto como um termo de dissipação artificial,
mas ao contrário do TVD não é auto-adaptativo.
Uin,j1 
1 p
t n
t n
Ui , j  Uic, j 
Di 1 2, j  Din1 2, j 
Di , j 1 2  Din, j 1 2
2
x
y






Universidade da
Beira Interior
VISCOSIDADE ARTIFICIAL
viscosidade artificial de Jameson (1981)



Din1 2, j  i 1 2, j Ui 1, j  Ui , j  i1 2, j Ui 2, j  3Ui 1, j  3Ui , j  Ui 1, j
2

 2
 2
 2
i 1 2, j  max i , j , i 1, j
 4
i 1 2, j
4

 2
i , j  
 2
 

i 1 2, j
4
 max 0,    
 
u x i 1 2, j  ci 1 2, j
 


 2
 u 
x i, j




 ci , j
zsi 1, j  2 zsi , j  zsi 1, j
z
si 1, j
 2 zsi , j  zsi 1, j

2
2
2
i , j1 2  max i , j , i, j1
 4
i , j 1 2

4

 2
i , j  
 2
2
 

i , j1 2
4


 max 0,   

u y i , j 1 2  ci , j 1 2  


 
4
   1 256
2
   1 4
Din, j 1 2  i , j1 2 Ui , j 1  Ui , j  i, j1 2 Ui , j 2  3Ui , j 1  3Ui , j  Ui , j 1
2
 u 
y i, j
 ci , j


zsi , j 1  2 zsi , j  zsi , j 1
z
si , j 1
 2 zsi , j  zsi , j 1
Universidade da
Beira Interior
TRATAMENTO DOS TERMOS DE FONTE
 O tratamento dos termos de fonte é um aspecto fundamental na
resolução numérica de equações.
 Existem duas formas de tratar os termos de fonte: i) aplicação de um
esquema de divisão (“splitting”ou “pointwise approach”); ii) aplicação
de um esquema “upwind”.
 No caso do esquema MacCormack a primeira alternativa é mais fácil de
implementar, dado que a segunda é facilmente implementada em
esquemas que usam discretização upwind do vector do fluxo, mas a
sua aplicação a esquemas centrados não é trivial.
Universidade da
Beira Interior
TRATAMENTO DOS TERMOS DE FONTE
Discretização de derivadas no termo de fonte
 Na discretização dos termos de fonte existe ainda outro problema
relacionado com a existência de derivadas parciais, como sejam os
declives do fundo zb/x e zb/y, que necessitam ser discretizadas
 Propõe-se uma discretização “upwind” de acordo com a do vector fluxo

 x

 x
 zbi , j 1  zbi , j  y

y  
 zbi , j  zbi , j 1  y
 zbi 1, j  zbi , j

zb x  
 zbi , j  zbi 1, j
zb

 x
F x   Fi , j  Fi 1, j  x
G y   G i , j 1  G i , j  y
G y   G i , j  G i , j 1  y
se F x  Fi 1, j  Fi , j
se
se
se
Universidade da
Beira Interior
CONDIÇÕES INICIAIS E DE FRONTEIRA
 A resolução de um sistema de equações às derivadas parciais através
da aplicação de um esquema numérico a um domínio de cálculo finito
exige a formulação de condições iniciais e de fronteira, que definem o
contorno desse domínio.
Universidade da
Beira Interior
CONDIÇÕES INICIAIS E DE FRONTEIRA
 O número de condições físicas iniciais ou de fronteira a ser
especificado em cada ponto do contorno do domínio de cálculo deve
ser igual ao número de características que “entram” nesse ponto
 O esquema MacCormak-TVD é “upwind” no tempo, pelo que em cada
ponto do contorno referente às condições iniciais entram todas as
características do sistema de equações. Assim, não são necessárias
condições iniciais numéricas, adoptando-se condições iniciais físicas
do tipo:
hm
h  x, y , t  0   
 h j
hsm
zb  x, y, t  0   
 hs j
se x  0
se x  0
se x  0
se x  0
q x   x, y, t  0   0
q y   x, y, t  0   0
Universidade da
Beira Interior
CONDIÇÕES INICIAIS E DE FRONTEIRA
 No que respeita às condições físicas em cada fronteira, a situação não
é tão simples como a das condições iniciais, pois o número de
características que entra em cada ponto do contorno depende do facto
de a fronteira se situar a montante ou a jusante, à esquerda ou à direita.
 As condições físicas devem ainda, em cada ponto do contorno, ser
estabelecidas de acordo com o tipo de informação (fase sólida ou fase
líquida) propagada pelas características que entram nesse ponto.
Assim, tem que se ter em conta a variação das características com o
número de Froude
Universidade da
Beira Interior
CONDIÇÕES INICIAIS E DE FRONTEIRA
Parede de montante:
três condições físicas referentes às
características positivas: 1 e 4,
associadas à fase líquida, e 3,
associada à fase sólida
q x  1, y, t   0
q y  1, y, t   0 C 1, y, t   0
uma condição numérica associada a 2
Parede de jusante:
uma condição física referentes à
característica negativa: 2,
associada à fase líquida
q  mx , t   0, 4 2 g h  mx , t   h j 
32
três condições numéricas associadas a 1 , 3 e 4
C  x,0, t   C  x,2, t 
Universidade da
Beira Interior
CONDIÇÕES INICIAIS E DE FRONTEIRA
Parede lateral esquerda:
três condições físicas referentes às
características positivas: 1 e 4,
associadas à fase líquida, e 3,
associada à fase sólida
q x   x, 0, t   q x   x, 2, t 
q y   x, 0, t   q y   x, 2, t 
C  x,0, t   C  x, 2, t 
uma condição numérica associada a 2
Parede lateral direita:
uma condição física referentes à
característica negativa: 2,
associada à fase líquida



q y  x, my  1, t  q y  x, my  1, t

três condições numéricas associadas a 1 , 3 e 4
C  x,0, t   C  x,2, t 
Universidade da
Beira Interior
CONDIÇÕES INICIAIS E DE FRONTEIRA
Parede que constitui o alargamento:
três condições físicas referentes às
características positivas: 1 e 4,
associadas à fase líquida, e 3,
associada à fase sólida
q x   mc  1, y, t   q x   mc  1, y, t  com
q y   mc  1, y, t   q y   mc  1, y, t  com
y  bm
y  bm
C  mc  1, y, t   C  mc  1, y, t  com y  bm
uma condição numérica associada a 2
Universidade da
Beira Interior
Instalação Experimental
9.5
1.0
1.0 0.8
1.0
1.0
lift-gate
0.5
downstream
2.0
0.5 0.5 0.8
perspex
wall
2.0
0.2
0.5
2.0
9.7
upstream
LEGEND:
lift mechanism
Planta do canal
localizado no LNEC
video camera
filming area
pressure tranducer
in the lateral wall
pressure transducer
in the bottom
Universidade da
Beira Interior
Condições iniciais
Lift-gate
foram realizados 2 tipos de
ensaios: fundo fixo e fundo
móvel (areia e pedra-pomes)
water
movable bed
0.07 m
0.40 m
Upstream
Downstream
Fundo de areia
(baixa mobilidade)
diâmetro mediano, d = 0.8 mm
densidade, s = 2.65
Fundo de pedra-pomes
(elevada mobilidade)
diâmetro mediano, d = 1.3 mm
densidade, s = 1.40
Universidade da
Beira Interior
Resultados experimentais
AREIA
Vista frontal
Universidade da
Beira Interior
Resultados experimentais
PEDRA-POMES
Vista frontal
Universidade da
Beira Interior
Resultados numéricos: Cota da sup. livre, zs
Componentes da velocidade, u(x) and u(y)
AREIA
zs
Vista lateral
u(x)
Planta
u(y)
Planta
Universidade da
Beira Interior
Resultados numéricos: Cota da sup. livre, zs
Componentes da velocidade, u(x) and u(y)
PEDRA-POMES
zs
Vista lateral
u(x)
Planta
u(y)
Planta
Universidade da
Beira Interior
Resultados numéricos
Cota da sup. livre, zs
T = 20
Areia
Pedra-pomes
Os resultados numéricos ajustam-se bem aos observados experimentalmente:
A reflexão na parede lateral origina um ressalto hidráulico.
O ensaio com pedra-pomes apresenta uma propagação longitudinal
mais lenta
Universidade da
Beira Interior
Resultados numéricos
Componentes da velocidade, u(x) e u(y)
Areia
T = 20
Pedra-pomes
u(x)
u(y)
Existe uma zona de recirculação, mais nítida no ensaio com pedra-pomes.
A propagação transversal é mais rápida no caso da pedra-pomes.