MMC44 - Modelagem e Simulação
Computacional em Recursos Hídricos
Propagação de cheias em rios
Prof. Benedito C. Silva
IRN/UNIFEI
Objetivo

Qual é o hidrograma em um local B a
jusante se o hidrograma em um local
A a montante é conhecido?
magnitude da vazão
 Níveis máximos
 tempo de ocorrência dos picos
 tempo de chegada da onda

?
Exemplo rio Uruguai
distância: 192 Km
Propagação

O que acontece com uma onde de
cheia enquanto viaja ao longo do
rio?
translação
 amortecimento
 contribuição de afluentes
 efeitos de jusante

Translação
A
B
Q
Hidrograma em A
Hidrograma em B
t
Amortecimento
A
B
Q
Hidrograma em A
Hidrograma em B
t
Efeitos de jusante
A
h em B (maré)
B
Q
Hidrograma em A
Hidrograma em B
t
Efeitos de jusante


Exemplo rio Amazonas
Entre Óbidos e Macapá



Óbidos (efeito ausente ou imperceptível)
Macapá (efeito da maré é
preponderante)
Aramanduba (ponto intermediário –
efeitos mistos)
26/12/97
11/12/97
26/11/97
11/11/97
27/10/97
12/10/97
27/9/97
12/9/97
28/8/97
13/8/97
29/7/97
14/7/97
29/6/97
14/6/97
30/5/97
15/5/97
30/4/97
15/4/97
31/3/97
16/3/97
1/3/97
14/2/97
30/1/97
15/1/97
31/12/96
Cota(cm)
Cotas no posto Obidos (17050000) no ano de 1997
(último posto c/ medidas de vazão no rio Amazonas)
900
800
700
600
500
400
300
200
100
0
-100
300
250
200
Cota(cm)
Cotas no posto Macapá (19500000) no ano de 1997
(último posto fluviométrico do rio Amazonas)
500
450
400
350
150
100
50
0
26/12/97
11/12/97
26/11/97
11/11/97
27/10/97
12/10/97
27/9/97
12/9/97
28/8/97
13/8/97
29/7/97
14/7/97
29/6/97
14/6/97
30/5/97
15/5/97
30/4/97
15/4/97
31/3/97
16/3/97
1/3/97
14/2/97
30/1/97
15/1/97
31/12/96
26/12/38
11/12/38
26/11/38
11/11/38
27/10/38
12/10/38
27/9/38
12/9/38
28/8/38
13/8/38
29/7/38
14/7/38
29/6/38
14/6/38
30/5/38
15/5/38
30/4/38
15/4/38
31/3/38
16/3/38
1/3/38
14/2/38
30/1/38
15/1/38
31/12/37
Cota(cm)
Cotas no posto Aramanduba (18400000) no ano de 1938
(Rio Amazonas, montante da confluência com rio Xingu, jusante de Óbidos)
600
500
400
300
200
100
0
Velocidade de propagação de
ondas de cheia


Ondas de cheia se propagam para
jusante com uma velocidade que é
maior do que a própria velocidade
média da água.
Assim, a velocidade de propagação
da onde de cheia em um rio cuja
velocidade média, durante uma
cheia, é de 1 m.s-1, é superior a 1
m.s-1, podendo chegar a 1,6 m.s-1,
por exemplo.
Celeridade

A velocidade de propagação das
ondas de cheia em rios pode ser
estimada pela celeridade
cinemática, que pode ser obtida
com base nas características médias
das seções transversais do rio e de
sua declividade.
Velocidade da onda cinemática

Uma abordagem mais intuitiva...
Ideia mais intuitiva sobre celeridade
de onda cinemática

Considerando que:

Então:
Q2
S f  S0
Q  f  A
Q1
Ideia mais intuitiva sobre celeridade
de onda cinemática

Imaginando uma onda de cheia gerada
por um aumento de vazão de Q1 para
Q2...
Frente da onda
Q2
Q1
Ideia mais intuitiva sobre celeridade
de onda cinemática

Depois de um pequeno intervalo de
tempo a onda deve ter se deslocado
para jusante
Frente da onda
t
Q2
t  t
Q1
Ideia mais intuitiva sobre celeridade
de onda cinemática

O volume adicional necessário para o
avanço da onda é:
Volume Q2  Q1  t
Frente da onda
t
Q2
t  t
Q1
Ideia mais intuitiva sobre celeridade
de onda cinemática
Ou,
Volume  A2  A1  x
Considerando que A1 é a área molhada da seção quando Q = Q1;
e que A2 é a área molhada da seção quando Q = Q2;
t
Q2
x
t  t
Q1
Ideia mais intuitiva sobre celeridade
de onda cinemática
Os dois volumes devem ser iguais, então:
Q2  Q1  t  A2  A1  x
t
Q2
x
t  t
Q1
Ideia mais intuitiva sobre celeridade
de onda cinemática
Reorganizando:
x Q2  Q1 

t  A2  A1 
Mas, x/t é a velocidade de avanço da onda...
t
Q2
x
t  t
Q1
Ideia mais intuitiva sobre celeridade
de onda cinemática
Se pensarmos numa variação de vazão bem pequena
x Q

t A
t
Q2
x
t  t
Q1
Assim...

A celeridade da onda cinemática
pode ser estimada por:
dx Q
c

dt A
Portanto a velocidade com que se propaga a onda de cheia
depende da relação entre vazão e área molhada, que é uma
característica da seção transversal, da declividade e da
rugosidade do rio ou canal.
Celeridade
c
combinando:
dQ
dA
Rh  S
Q  u  A  A
n
2
3
1
2
pode-se obter:
c
5
u
3
O que significa que a celeridade
(velocidade de propagação da onda de cheia)
é superior à velocidade média da água.
Celeridade


A velocidade de propagação das
ondas de cheia é maior do que
a própria velocidade média da
água
Além disso, a velocidade de
propagação das cheias tende a
ser maior para cheias maiores,
porque o nível da água e a
velocidade média tendem a ser
maiores.
5
c  u
3
Celeridade

Por outro lado, em rios com grandes
planícies de inundação, a velocidade
de propagação das ondas de cheia
tende a diminuir drasticamente no
momento em que o rio começa a
transbordar.
Celeridade diminui
Celeridade aumenta
Celeridade
Evidências experimentais
Murrumbidgee river - Wang e Laurenson, 1983 Water Resources Research
Tempo de viagem
192 km
mais ou menos 2 dias
Cálculos de propagação


Modelos hidrodinâmicos
Modelos simplificados
Equações de Saint-Venant

As equações utilizadas para
descrever o comportamento do
escoamento em rios e canais foram
inicialmente derivadas no século
XIX
Hipóteses assumidas
1.
O escoamento é unidimensional; a
velocidade é uniforme e igual à
média; a nível da água é
horizontal na seção transversal.



Escoamento em meandros e transições
é fortemente tridimensional
Velocidade é maior no centro da seção
Em curvas o nível da água pode não ser
horizontal
Hipóteses assumidas
2.
Pressão é hidrostática (depende
apenas da profundidade)
Variações de forma da seção devem ser
relativamente suaves.
Hipóteses assumidas
3.
4.
É possível usar fórmulas para
perda de carga semelhantes às
usadas em escoamento
permanente (como Manning).
A declividade do canal é pequena,
o cosseno do ângulo é quase igual
a 1.
Continuidade ou conservação de
massa
Considere um volume de controle entre as seções x=x1 e x=x2
e ao longo de um intervalo de tempo de t=t1 a t=t2
x1
x2
A
A
Continuidade ou conservação de
massa de água:
x2
 A  A dx
t2
t1
t2
 uA  uA dt
=
x1
x1
x2
t1
considerando que Q=u.A
e que a massa específica da água é constante:
x2
t2
  A   A dx   Q
t2
x1
t1
x2
t1
 Q x1 dt  0
Forças agindo sobre o volume de
controle



Elevation View


Plan View
Fg = Força da
gravidade
Ff = Força de atrito
com o fundo e
margens
Fe = Força devida
às contrações e
expansões da seção
Fw = força de atrito
do vento na
superfície
Fp = diferença de
pressão nos limites
de montante e
jusante do volume
de controle
Equações na forma integral
 uA  uA dx   u
x2
t2
t2
t1
x1
2
 
t1
x2
t2
  A   A dx   Q
t2
x1




 
 A x1  u  A x 2 dt  g  I1 x1  I1 x 2 dt  g   I 2 dxdt  g   AS0  S f dxdt
2
t1
x2
t2
t 2 x2
t 2 x2
t1
t1 x1
t1 x1
 Q x1 dt  0
t1
Equações estabelecidas sem exigir que A; Q; u fossem
variáveis contínuas e diferenciáveis.
Também não exigem que x2-x1 seja uma distância
infinitesimalmente pequena.
Alguns esquemas numéricos estão baseados na forma
integral das equações, em vez de usar a forma diferencial.
Método dos volumes finitos se baseia em equações na
forma integral.
Equações na forma integral

Considerando
variáveis contínuas
e diferenciáveis e
usando expansão
em série de Taylor
pode se chegar a
outra forma:
 dtdx   g 
 Q  u 2 A
x1 t1  t  x
x 2t 2


 Q A 
x1 t1  x  t dtdx  0
x 2t 2
 I1




I

A
S

S
0
f dtdx
x1 t1  x 2

x 2t 2
Equações na forma diferencial

Em um volume
infinitamente
pequeno as
equações anteriores
também devem
valer, assim:
Q A

0
x t
 
Q  u 2 A
I1

g
 gAS0  S f   gI2
t
x
x
Equações na forma diferencial


considerando
que
Q=u.A
Q A

0
x t

Q   Q 2
 
 gI1   gAS0  S f   gI2
t x  A

Equações na forma diferencial

combinando os termos
com I1 e I2
Q A

0
x t
Q   Q 2 
 dh

    gA  S0   gASf  0
t x  A 
 dx

que não é mais exatamente uma equação de conservação de quantidade de
Movimento, muitas vezes chamada de “equação dinâmica”
h é a profundidade
Equações de Saint-Venant na forma
mais usual atualmente
A Q

0
t x
Q   Q 2 
h
  g  A   g  A  S f  0
 
t x  A 
x
Sf 
Q Q  n
A R
2
2
3
Vazão e nível da água
Solução das equações de
Saint-Venant



Não existem soluções analíticas para as
equações de Saint-Venant na maior parte
das aplicações úteis.
Somente nas décadas mais recentes é que
os métodos numéricos e os computadores
digitais permitiram a solução das equações
completas de Saint-Venant.
Atualmente existem diversos programas
computacionais de modelos matemáticos
que resolvem as equações de Saint-Venant
numericamente para resolver problemas de
propagação de vazão em rios e canais.

Zij1  Zij11  Zij  Zij1
Z

t
2  t
Aplicando este
esquema:




Zij11  Zij1
Zij1  Zij
Z
 
 1   
x
xi
xi

Z
Z  
j 1
i 1 
Zij1
2
  1    Z
j
i 1 
Zij


2
Q A

q  0
x
t
A estas equações:
Q     Q2 
 h


 g A 
 Sf   0
t x  A 
 x

Sf 
n2  Q  Q
2
A R
4
3
Equações de diferenças
numéricas
continuidade
Q

j1
j1

Q
i1
i
xi
  1  Q
j
j

Q
i1
i
Q
x 
i
Dinâmica
xi
  A
j1
j1
j
i  Qi1  Qi
2  t
j1
j1
j

A

A
i
i1
i
 Aij1
2  t
 Qij1
 q
j1
j

q
i
i
2
0

j1
j1


2
2



Q


Q
 
  g  A j1  h j1  h j1  x  S j1 
  
i
f
i1
i
 A 

 A 




i1
i




j
j


2
   Q2 


Q
 
  g  A j  h j  h j  x  S j   0
1    
i
f 
i1
i

 A 
A
i1 
i



A  0.5  Ai  Ai1
R  0.5  Ri  Ri1
Sf 
n2  Q  Q
A2  R
4
3

Incógnitas – não linear
continuidade
Q

j1
j1

Q
i1
i
xi
  1  Q
j
j

Q
i1
i
Q
x 
i
Dinâmica
xi
  A
j1
j1
j
i  Qi1  Qi
2  t
j1
j1
j

A

A
i
i1
i
 Aij1
2  t
 Qij1
 q
j1
j

q
i
i
2
0

j1
j1


2
2



Q


Q
 
  g  A j1  h j1  h j1  x  S j1 
  
i
f
i1
i
 A 

 A 




i1
i




j
j


2
   Q2 


Q
 
  g  A j  h j  h j  x  S j   0
1    
i
f 
i1
i

 A 
A
i1 
i



A  0.5  Ai  Ai1
R  0.5  Ri  Ri1
Sf 
n2  Q  Q
A2  R
4
3

Esquema de Preissmann

Cada trecho entre
duas seções define
duas equações:



continuidade
dinâmica
Cada seção tem
duas incógnitas:

h e Q no tempo
futuro
Modelos hidrodinâmicos

Atualmente existem programas de
modelos como o HEC-RAS que
podem ser utilizados para os
cálculos de propagação de cheias
em rios
Simulação hidrodinâmica
Modelos simplificados


Em função da dificuldade que havia para
resolver as equações de Saint-Venant, um
grande número de métodos simplificados
foi criado para calcular os efeitos que
ocorrem em ondas de cheia à medida que
se propagam ao longo de rios
Estes métodos utilizam a equação da
continuidade mas simplificam ao máximo a
equação da conservação da quantidade de
movimento
Escoamento em rios
Modelo Muskingum

dS
 I Q
dt
S  f ( I , Q)
Criado na década
de 1930 por
McCarthy para
representar a
propagação de
vazão ao longo do
rio Muskingum.
Supõe que S
(armazenamento) está
relacionado a I (vazão de
entrada) e Q (vazão de saída)
Escoamento em rios: Muskingum

Continuidade

Relação
dS
 IQ
dt
S = K [X I +(1-X) Q]
A vazão (Q) na seção de jusante é
dada por:
Qt 1  C1It 1  C2 I t  C3Qt
t
t
t
KX 
K (1  X) 
2 ; C 
2 ; C 
2
C1 
2
3
t
t
t
K (1  X) 
K (1  X) 
K (1  X) 
2
2
2
C1+C2+C3=1
K é o tempo
médio de
deslocamento
da onda
X é um
ponderador
entre as
vazões de
entrada e
saída
 KX 
52
Intervalo de tempo

Para que os coeficientes da equação sejam positivos
t
2  0 e 2KX  t
C1 
t
K(1  X) 
2
 KX 
2KX  t  2K(1  X)
0  X  0,5
t é o intervalo de tempo para
simulação da propagação
t
2  0 e 2K(1- X)  t
C3 
t
K (1  X) 
2
K (1  X) 
t
2X 
 2(1  X)
K
t / K
2
1
Região válida
0
0
0,5
X
53
KeX

O método Muskingum tem dois
parâmetros de cálculo (K e X) que
devem ser definidos antes dos
cálculos.
Definir valor de X


O parâmetro X é um ponderador
adimensional cujo valor deve estar
entre 0 e 0,5, mas na maior parte
dos rios e canais naturais seu valor é
próximo a 0,3.
Dependendo do valor de X ocorre
mais ou menos amortecimento da
onda de cheia.


Para um valor de X igual a 0,5 não ocorre
amortecimento.
Quando X é igual a zero o amortecimento
é máximo.
Efeito de X
Efeito de K
Definir K



O parâmetro K têm unidades de tempo e
deve ser expresso nas mesmas unidades de
t.
O valor de K pode ser estimado pelo tempo
de viagem do pico da cheia do início ao
final do trecho de rio, ou seja, a distância
dividida pela celeridade.
Quanto maior o valor de K, mais afastados
no tempo ficam os picos de vazão na
entrada e saída do trecho de canal.
Estimativa de K e X
Tradicional método da laçada
S/∆t
X=X1
X= Xn
Variar o valor de X até
que se crie uma laçada,
com forma mais
próxima possível de
uma reta
Ajustar uma linha de
tendência linear
S/∆t = a. QI + b
K=a
K será igual ao
coeficiente angular da
reta
QI
𝑆𝑡+1
𝑆𝑡
= 0,5 𝐼𝑡+1 + 𝐼𝑡 − 𝑄𝑡+1 + 𝑄𝑡 +
∆𝑡
∆𝑡
Método da laçada - Exemplo
A tabela abaixo apresenta os hidrogramas de vazões
medidos nas seções de entrada e saída de um trecho de
rio. Determine os valores dos parâmetros K e X
Tempo
(h)
0
6
12
18
24
30
36
42
48
54
60
66
72
78
84
90
Entrada
(m3/s)
30
120
286
412
373
306
246
198
165
141
123
108
93
81
72
63
Saída
(m3/s)
30
39
45
93
181
237
264
261
246
225
202
184
174
153
135
117
60
Muskingum-Cunge
Adaptado para estimativa com base em
parâmetros físicos do trecho
Qo
X  0,5(1 
)
bo So co x
0,3
0, 4
5 S 0 Q0
co 
3 bo 0, 4 .n 0, 6
Ou:
x
K
co
Q0
0 ,8
x 
 0,8.c0 .t  x 0, 2
b0 .S0 .c0
1/2




c0 .t 
Q0
x 
1  1  1,5.
 
2
2  
b0 .S0 .t.c 0  


Roteiro de Ajuste
1) Fixar ∆t = tp/5 ou outro valor para ∆t ≤ tp/5
2) Adotar valor de Qo = 2/3 da vazão máxima do
hidrograma de entrada
3) Calcular co
4) Calcular ∆x por processo iterativo
2,5Qo
5) A primeira estimativa de ∆x pode ser obtida por xo 
bo So co
6) Calcular K e X e verifique se está dentro da
faixa de validade
7) Caso contrário modifique ∆x
62
Exemplo


Determine o hidrograma 18 km a jusante de uma seção de um rio.
As características do trecho são: largura=30m, declividade=0,0007
m/m; rugosidade de Manning n=0,045.
o tempo tp = 240 min e =240/5=48 min, ∆t=40min. A vazão
máxima de montante é 130 m3/s, Qo=87m3/s
5 So0,3Qo0, 4
co 
 1,86m / s
0, 6
0, 4
3 n bo
x 
2,5.87
 5.568m
30x0,0007x1,86
Por convergência
x  6018 m
x  6000 m adotado
X=0,31
K = 1,34
t
Tempo
(40min)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
vazão de entrada
m3 / s
20
30
60
90
100
130
115
95
80
60
40
20
20
20
20
vazão de saída
m3 / s
20
20
20
20
21,1
27,0
42,2
63,9
85,9
103,0
102,4
92,4
77,2
59,4
41,9
63
Muskingum Cunge não linear



A celeridade não é constante
Os parâmetros do método de
Muskingum Cunge deveriam variar
Celeridade varia com o nível da
água ou com a vazão
Celeridade diminui
Celeridade aumenta
Muskingum Cunge não linear



Substituir K e X (C1, C2 e C3)
constantes por variáveis
A cada passo de tempo é necessário
recalcular o valor de K e X (C1, C2
e C3)
Só o que não muda é o x
Muskingum Cunge não linear



Qual vazão usar como referência?
Criar tabela Q x C a partir de tabela
hxAxQ
Solução não-linear

Cálculo de X e K em cada célula de
cálculo
t
t+1
It+1
Calcular K e X
com base em:
Qt+1
(1) Qt
t
It
Qt
(2) Qt, It e It+1
(3) todos.
i
i+1
x
Qot 
I t  Qt  I t 1
3
67
Exemplo Jacuí

Linear x Não-linear
Evento
1
2
3
4
Linear
0,91
0,83
0,92
0,88
Não-Linear
0,97
0,94
0,96
0,98
68