Equação de Águas Rasas

Em construção Grupo: Gabriel Schmökel, Julia Remus e Pedro Inocêncio Rodrigues Terra

Introdução

Tsunami é um fenômeno da natureza caracterizado por uma sucessão de ondas marinhas, que devido ao seu grande volume e alta velocidade, podem se tornar catastróficas ao atingir a costa. Sismos, erupções vulcânicas, deslizamentos de terra, impactos e outros movimentos submarinos são a causa para a formação deste evento, sendo a grande maioria provocado pelos movimentos das placas tectônicas.

Formação de um Tsunami

Vamos analisar a sequência de passos da formação de uma Tsunami formada a partir de um abalo sísmico:

I. A convergência das placas tectônicas, devido as correntes de convecção, faz com que existam forças de tensão entre as placas.

IMAGEM

A tensão entre as placas eventualmente ultrapassa o limite máximo, o que provoca o deslizamento brusco de uma das placas sobre a outra, gerando um grande deslocamento de volume de água na vertical. Como a tsunami ocorre em grandes profundidades, ela pode passar despercebida para um barco que navega nas proximidades, uma vez que amplitude da onda é menor.

IMAGEM

II. A onda gerada se propaga ao longo de todas as direções do plano da água.

IMAGEM

III. A medida que a onda se aproxima da superfície ela diminui sua velocidade e aumenta sua amplitude

IMAGEM

Temos o interesse de descrever fisicamente a propagação da Tsunami de acordo com a topografia da água e do mar, por essa razão não iremos estudar o efeito físico que causou o deslocamento do volume de água.

Teoria

Derivação das EQs. de Águas Rasas

Para obter as equações de águas rasas devemos partir da equação da continuidade e das equações da quantidade de movimento de Navier-Stokes:

${\frac {d\rho }{dt}}+\nabla .(\rho \mathbf {u} )=0\qquad (3)$

${\frac {D\mathbf {u} }{Dt}}+{\frac {1}{\rho }}\nabla p+{\frac {1}{\rho }}\nabla .{\boldsymbol {\tau }}+\mathbf {g} =0\qquad (4)$

$\rho$ é a densidade; p é a pressão; $\mathbf {u} =(u,v,w)$ é o vetor velocidade do fluído, onde u,v e w são as velocidades das partículas que compõe o fluído nas direções x,y,z; $\mathbf {g}$ é o vetor aceleração da gravidade; ${\boldsymbol {\tau }}$ é o tensor tensão, onde as componentes deste tensor são as tensões normais e tangenciais de cisalhamento, expressas por $\tau _{ij}$, no qual $i$ indica a direção e $j$ o plano normal.

Introduzindo as condições de contorno para a superfície $z(x,y,t)$ e para a profundidade do oceano $h(x,y)$:

${\frac {D\eta }{Dt}}={\frac {\partial \eta }{\partial t}}+\mathbf {v} .\nabla \eta =w$ , onde $z=\eta (x,y,t)\qquad (4)$

$\mathbf {u} .\nabla (z+h(x,y))=0$ , onde $z=-h(x,y)\qquad (5)$

$\eta$ é o deslocamento vertical da água sobre a superfície em repouso, $\mathbf {v} =(x,y,0)$ é o vetor velocidade do fluído nas direções horizontais x e y.

A equação da continuidade em (3) pode ser simplificada, já que a densidade do fluído no oceano $\rho$ não varia significativamente com o tempo e a posição.

$\nabla .\mathbf {u} =0\qquad (6)$

Integrando a expressão da continuidade em (6), utilizando a regra da integral de Leibniz , com os limites indo de $-h(x,y)$ até $\eta (x,y,t)$ chegamos na seguinte expressão:

$\int _{-h}^{\eta }\nabla .\mathbf {u} =\int _{-h}^{\eta }{\Big (}{\frac {\partial u}{\partial x}}+{\frac {\partial v}{\partial y}}+{\frac {\partial w}{\partial z}}{\Big )}dz={\frac {\partial }{\partial x}}\int _{-h}^{\eta }udz+{\frac {\partial }{\partial y}}\int _{-h}^{\eta }vdz+w{\Big |}_{-h}^{\eta }+\mathbf {u} .\nabla (z+h(x,y)){\Big |}_{-h}^{\eta }-u{\Big |}_{-h}^{\eta }{\frac {\partial \eta }{\partial x}}-v{\Big |}_{-h}^{\eta }{\frac {\partial \eta }{\partial y}}\qquad (7)$

Teorema de Leibniz:

${\frac {d}{dx}}\left(\int _{a(x)}^{b(x)}f(x,t)\,dt\right)=f{\big (}x,b(x){\big )}\cdot {\frac {d}{dx}}b(x)-f{\big (}x,a(x){\big )}\cdot {\frac {d}{dx}}a(x)+\int _{a(x)}^{b(x)}{\frac {\partial }{\partial x}}f(x,t)\,dt\qquad (8)$

Substituindo as condições de contorno da profundidade (5) em (7) obtemos:

${\frac {\partial }{\partial x}}\int _{-h}^{\eta }udz+{\frac {\partial }{\partial y}}\int _{-h}^{\eta }vdz-w{\Big |}_{eta}-\mathbf {v} .\nabla \eta =0\qquad (9)$

Substituindo a condição de contorno da superfície (4) em (9):

${\frac {\partial }{\partial x}}\int _{-h}^{\eta }udz+{\frac {\partial }{\partial y}}\int _{-h}^{\eta }vdz+{\frac {\partial \eta }{\partial t}}={\frac {\partial u(\eta +h)}{\partial x}}+{\frac {\partial v(\eta +h)}{\partial y}}+{\frac {\partial \eta }{\partial t}}={\frac {\partial uD}{\partial x}}+{\frac {\partial vD}{\partial y}}+{\frac {\partial \eta }{\partial t}}$

$\Rightarrow {\frac {\partial \eta }{\partial t}}+{\frac {\partial uD}{\partial x}}+{\frac {\partial vD}{\partial y}}=0\qquad (10)$

(10) é a primeira das equações das águas rasas que obtemos, onde $D$ é o comprimento da água total do fundo do oceano até a amplitude da onda. Podemos expressar (10) através do fluxo de descarga nas direções x e y, estás quantidades estão relacionadas com as velocidades da seguinte forma :

$M={\frac {\partial }{\partial x}}\int _{-h}^{\eta }udz=uD\qquad (11)$

$N={\frac {\partial }{\partial y}}\int _{-h}^{\eta }vdz=vD\qquad (12)$

Substituindo (11) e (12) em (10) chegamos na representação do fluxo de descarga para uma das equações de águas rasas.

$\Rightarrow {\frac {\partial \eta }{\partial t}}+{\frac {\partial M}{\partial x}}+{\frac {\partial N}{\partial y}}=0\qquad (10)$

Escrevendo as quantidades de movimento de Navier-Stokes nas componentes x,y e z:

Falhou ao verificar gramática (erro de sintaxe): {\displaystyle \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y}          + w\frac{\partial u}{\partial z} +\frac{1}{\rho}\frac{\partial P}{\partial x} +g_x = 0 \qquad (14) }

Falhou ao verificar gramática (erro de sintaxe): {\displaystyle \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y}          + w\frac{\partial v}{\partial z} +\frac{1}{\rho}\frac{\partial P}{\partial x} +g_y = 0 \qquad (15) }

${\frac {1}{\rho }}{\frac {\partial P}{\partial x}}+g_{z}=0\qquad (16)$

Na componente z em (15) negligenciamos a aceleração das partículas, pois a aceleração da gravidade é muito maior. Também tomamos como nulos as componentes $g_{x}$ e $g_{y}$ em (14) e passamos a definir $g_{z}=g$.

Resolvendo equação diferencial da componente z em (16) podemos obter a pressão, a qual é hidrostática.

$\partial P=\rho g\partial z\Rightarrow P=\rho g(\eta -z)\qquad (17)$

Substituindo a pressão em (14):

${\frac {\partial u}{\partial t}}+u{\frac {\partial u}{\partial x}}+v{\frac {\partial u}{\partial y}}+w{\frac {\partial u}{\partial z}}+g{\frac {\partial \eta }{\partial x}}=0\qquad (18)$

Integrando a equação (18) em relação a componente z com os limites indo

Integrando a expressão (18), utilizando a regra da integral de Leibniz [1] e as condições de contorno (4) e (5), com os limites indo de $-h(x,y)$ até $\eta (x,y,t)$ chegamos em outra das equações de águas rasas:

${\frac {\partial uD}{\partial t}}+{\frac {\partial u^{2}D}{\partial x}}+{\frac {\partial uvD}{\partial y}}+{\frac {g}{2}}{\frac {\partial D^{2}}{\partial x}}=0\qquad (18)$

Generalizando a equação (18), para a componente y, obtemos a última das equações de águas rasas:

${\frac {\partial vD}{\partial t}}+{\frac {\partial v^{2}D}{\partial y}}+{\frac {\partial uvD}{\partial x}}+{\frac {g}{2}}{\frac {\partial D^{2}}{\partial y}}=0\qquad (19)$

Na representação de fluxo de cargas as expressões (18) e (19) são apresentadas respectivamente como:

${\frac {\partial M}{\partial t}}+{\frac {\partial }{\partial x}}{\Big (}{\frac {M^{2}}{D}}{\Big )}+{\frac {\partial }{\partial y}}{\Big (}{\frac {MN}{D}}{\Big )}+gD{\frac {\partial \eta }{\partial x}}=0\qquad (20)$

${\frac {\partial N}{\partial t}}+{\frac {\partial }{\partial y}}{\Big (}{\frac {N^{2}}{D}}{\Big )}+{\frac {\partial }{\partial x}}{\Big (}{\frac {MN}{D}}{\Big )}+gD{\frac {\partial \eta }{\partial x}}=0\qquad (21)$

Iremos escrever as equações das águas rasas considerando o tensor de estresse ${\boldsymbol {\tau }}$. Os elementos deste tensor são responsáveis por causar nas partículas tensões tangenciais e perpendiculares, onde as tensões tangenciais são representadas por elementos $\tau _{ij}$ onde $i\neq j$, e as perpendiculares por elementos $\tau _{ij}$ onde $i=j$

Decompondo nas componentes x,y, e z de ${\frac {1}{\rho }}\nabla .{\boldsymbol {\tau }}$ presente em (4):

${\frac {1}{\rho }}{\Big (}{\frac {\partial }{\partial x}}\tau _{xx}+{\frac {\partial }{\partial y}}\tau _{xy}+{\frac {\partial }{\partial z}}\tau _{xz}{\Big )}\qquad (22)$

${\frac {1}{\rho }}{\Big (}{\frac {\partial }{\partial x}}\tau _{yx}+{\frac {\partial }{\partial y}}\tau _{yy}+{\frac {\partial }{\partial z}}\tau _{yz}{\Big )}\qquad (23)$

${\frac {1}{\rho }}{\Big (}{\frac {\partial }{\partial x}}\tau _{zx}+{\frac {\partial }{\partial y}}\tau _{zy}+{\frac {\partial }{\partial z}}\tau _{zz}{\Big )}\qquad (24)$

Considerando o fluído Newtoniano, então as tensões de cisalhamento serão proporcionais a uma taxa de deformação, onde a constante de deformidade é a viscosidade.

$\tau _{xy}=\tau _{yx}=\mu {\Big (}{\frac {\partial u}{\partial y}}+{\frac {\partial v}{\partial x}}{\Big )}\qquad (25)$

$\tau _{xz}=\tau _{zx}=\mu {\Big (}{\frac {\partial u}{\partial z}}+{\frac {\partial w}{\partial x}}{\Big )}\qquad (26)$

$\tau _{yz}=\tau _{zy}=\mu {\Big (}{\frac {\partial w}{\partial y}}+{\frac {\partial v}{\partial z}}{\Big )}\qquad (27)$

Substituindo (25),(26) em (25), integrando em relação a componente z, utilizando a regra de Leibnz e as condições de contorno (3) e (4), obtemos:

$\int _{-h}^{\eta }{\frac {1}{\rho }}{\Big (}{\frac {\partial }{\partial x}}\tau _{xx}+{\frac {\partial }{\partial y}}\tau _{xy}+{\frac {\partial }{\partial z}}\tau _{xz}{\Big )}={\frac {\tau _{x}}{\rho }}-A{\Big (}{\frac {\partial ^{2}}{\partial x^{2}}}M{\frac {\partial ^{2}}{\partial x^{2}}}M{\Big )}\qquad (28)$

Onde $A$ é a constante de viscosidade turbulenta, $\tau _{x}$ é uma força de resistividade relativa ao movimento do fluído com o fundo do oceano na direção x. Podemos negligenciar a constante de turbulência na situação em que não temos inclinações abrutas no fundo do mar. .

Considerando que o fluído é uniforme, então a expressão para Falhou ao verificar gramática (erro de sintaxe): {\displaystyle \frac{\tau_x}{\rho} é } é:

${\frac {\tau _{x}}{\rho }}={\frac {fM}{2D^{2}}}(M^{2}+N^{2})^{1/2}\qquad (29)$

$f$ é o coeficiente de fricção, porém o coeficiente de rugosidade de Manning $n$ é mais usado, alguns valores deste coeficiente são:

• Cimento puro e metal liso $n=0,010$
• Terra lisa $n=0,017$
• Pedras, ervas daninhas $n=0,035$
• Péssimo relevo de canal $n=0,060$
• Bom relevo de canal $n=0,025$

O coeficiente de fricção $f$ e o de rugosidade de Meanning $n$ estão relacionados por:

$f={\frac {2gn^{2}}{D^{1/3}}}\qquad (30)$

Substituindo (30) em (29) obtemos:

${\frac {\tau _{x}}{\rho }}={\frac {gn^{2}}{D^{7/3}}}M(M^{2}+N^{2})^{1/2}\qquad (31)$

Generalizando a expressão (31) para a componente y. ${\frac {\tau _{y}}{\rho }}={\frac {gn^{2}}{D^{7/3}}}N(M^{2}+N^{2})^{1/2}\qquad (32)$

Adicionando, repectivamente, (31) e (32) nas expressões (20) e (21), obtemos as equações de águas rasas considerando as forças de fricção do fundo do oceano.

${\frac {\partial M}{\partial t}}+{\frac {\partial }{\partial x}}{\Big (}{\frac {M^{2}}{D}}{\Big )}+{\frac {\partial }{\partial y}}{\Big (}{\frac {MN}{D}}{\Big )}+gD{\frac {\partial \eta }{\partial x}}{\frac {gn^{2}}{D^{7/3}}}M(M^{2}+N^{2})^{1/2}=0\qquad (33)$

${\frac {\partial N}{\partial t}}+{\frac {\partial }{\partial y}}{\Big (}{\frac {N^{2}}{D}}{\Big )}+{\frac {\partial }{\partial x}}{\Big (}{\frac {MN}{D}}{\Big )}+gD{\frac {\partial \eta }{\partial x}}{\frac {gn^{2}}{D^{7/3}}}N(M^{2}+N^{2})^{1/2}=0\qquad (34)$

Forma Conservativa

Um modelo mais simples - desconsiderando a fricção, viscosidade do líquido e as forças de Coriolis sobre ele - pode ser obtido . Para desenvolvê-lo são necessárias algumas premissas:

• O comprimento da onda é muito maior que as contribuições na direção ${\vec {z}}$
• A aceleração na direção da velocidade na direção ${\vec {z}}$ é zero
• As componentes das velocidades em ${\vec {x}}$ e em ${\vec {y}}$ (${\vec {u}}$ e ${\vec {v}}$) não variam em ${\vec {z}}$

O sistema então pode ser descrito pelas seguintes equações:

${\dfrac {\partial h}{\partial t}}+{\dfrac {\partial hu}{\partial x}}+{\dfrac {hv}{\partial y}}=0$

${\dfrac {\partial hu}{\partial t}}+{\dfrac {\partial \left(hu^{2}+{\dfrac {1}{2}}gh^{2}\right)}{\partial x}}+{\dfrac {huv}{\partial y}}=-gh{\dfrac {\partial b}{\partial x}}$

${\dfrac {\partial hv}{\partial t}}+{\dfrac {huv}{\partial x}}+{\dfrac {\partial \left(hv^{2}+{\dfrac {1}{2}}gh^{2}\right)}{\partial y}}=-gh{\dfrac {\partial b}{\partial y}}$

Onde $h$ é a altura do fluido desde a base, ${\vec {u}},{\vec {v}}$ são as velocidades médias na direções ${\vec {x}},{\vec {y}}$, $g$ é a constante gravitacional e $b(x,y)$ é função que descreve a superfície onde acontece o movimento .

Forma dissipativa

As equações de águas rasas na forma não conservativa são dadas por (10),(33) e (34). Para descrever numericamente o fenômeno foi utilizado discretização por diferenças finitas, onde realizamos derivadas centradas na região espacial, e para frente no região temporal. O erro de truncamento é de ordem $\Delta {\mathcal {O}}(x^{2})$ na região espacial, enquanto na temporal é de ordem $\Delta {\mathcal {O}}(t^{1})$. O método é conhecido como leap-frog method devido a discretização central na região espacial.

Discretizando a expressão (10) pelo leap-frog method:

${\frac {\partial \eta }{\partial t}}+{\frac {\partial uD}{\partial x}}+{\frac {\partial vD}{\partial y}}=0$

$\Rightarrow {\frac {\eta (x,y,t+\Delta t)-\eta (x,y,-t)}{\Delta t}}=-{\frac {M(x+\Delta x,y,t)-M(x-\Delta x,y,t)}{2\Delta x}}-{\frac {N(x,y+\Delta x,t)-N(x,y-\Delta y,t)}{2\Delta y}}$

$n_{i,j}^{n+1}=n_{i,j}^{n}-dt{\Bigg (}{\frac {M_{i+1,j}^{n}-M_{i-1,j}^{n}}{2dx}}+{\frac {M_{i,j+1}^{n}-M_{i,j-1}^{n}}{2dy}}{\Bigg )}$

Discretizando a expressão (33) pelo leap-frog method:

$\Rightarrow {\frac {\partial M}{\partial t}}=-{\Bigg (}{\frac {\partial }{\partial x}}{\Big (}{\frac {M^{2}}{D}}{\Big )}+{\frac {\partial }{\partial y}}{\Big (}{\frac {MN}{D}}{\Big )}+gD{\frac {\partial \eta }{\partial x}}+{\frac {gn^{2}}{D^{7/3}}}M(M^{2}+N^{2})^{1/2}{\Bigg )}\qquad$

Definindo as quantidades:

$M2\equiv {\frac {M^{2}(x,y,t)}{D(x,y,t)}}$

$MN\equiv {\frac {M(x,y,t)N(x,y,t)}{D(x,y,t)}}$

$f(x,y,t)\equiv {\frac {gn^{2}}{D^{7/3}}}M(M^{2}+N^{2})^{1/2}$

Das quantidades definidades e da derivada parcial do fluxo de descarga em relação ao tempo temos a respectiva avanço temporal para M:

$M_{i,j}^{n+1}=M_{i,j}^{n}-dt{\Bigg (}{\frac {M2_{i+1,j}^{n}-M2_{i-1,j}^{n}}{2x}}+{\frac {(MN)_{i,j+1}^{n}-(MN)_{i,j-1}^{n}}{2y}}+gD_{i,j}^{n}{\frac {\eta _{i+1,j}^{n}-\eta _{i-1,j}^{n}}{2x}}+f_{i,j}^{n}{\Bigg )}$

Generalizando a expressão anterior para o fluxo de descarga N temos:

$N_{i,j}^{n+1}=N_{i,j}^{n}-dt{\Bigg (}{\frac {N2_{i,j+1}^{n}-M2_{i,j-1}^{n}}{2y}}+{\frac {(MN)_{i+1,j}^{n}-(MN)_{i-1,j}^{n}}{2y}}+gD_{i,j}^{n}{\frac {\eta _{i+1,j}^{n}-\eta _{i-1,j}^{n}}{2y}}+f_{i,j}^{n}{\Bigg )}$

Desenvolvimento do cálculo

Forma conservativa 2D

Para descrever numericamente o fenômeno foi utilizado discretização por diferenças finitas e o método pra frente no tempo e no espaço (FTCS). As equações discretizadas podem ser observadas abaixo.

${\dfrac {h_{i,j}^{t+\Delta t}-h_{i,j}^{t}}{\Delta t}}+\left[{\dfrac {(hu)_{i+1,j}^{t}-(hu)_{i-1,j}^{t}}{2\Delta x}}\right]+\left[{\dfrac {(hv)_{i,j+1}^{t}-(hv)_{i,j-1}^{t}}{2\Delta y}}\right]=0$

${\dfrac {hu)_{i,j}^{t+\Delta t}-(hu)_{i,j}^{t}}{\Delta t}}+\left[{\dfrac {(hu^{2}+{\cfrac {1}{2}}gh^{2})_{i+1,j}^{t}-(hu^{2}+{\cfrac {1}{2}}gh^{2})_{i-1,j}^{t}}{2\Delta x}}\right]+\left[{\dfrac {(huv)_{i,j+1}^{t}-(huv)_{i,j-1}^{t}}{2\Delta y}}\right]=-gh_{i,j}^{t}b_{x.i,j}$

${\dfrac {(hv)_{i,j}^{t+\Delta t}-(hv)_{i,j}^{t}}{\Delta t}}+\left[{\dfrac {(huv)_{i+1,j}^{t}-(huv)_{i-1,j}^{t}}{2\Delta x}}\right]+\left[{\dfrac {(hv^{2}+1/2gh^{2})_{i,j+1}^{t}-(hv^{2}+1/2gh^{2})_{i,j-1}^{t}}{2\Delta y}}\right]=-gh_{i,j}^{t}b_{y.i,j}$

Para os contornos foi utilizado que:

• $u(x_{f},y)=-u(x_{f}-\Delta x,y)$

• $u(x_{i},y)=-u(x_{i}+\Delta x,y)$

• $v(x,y_{f})=-v(x,y_{f}-\Delta y)$

• $v(x,y_{i})=-v(x,y_{i}+\Delta y)$

• Nos contornos de x: ${\tfrac {\partial h}{\partial x}}=0$, discretizando essa derivada temos que: $h_{i(+/-)1,j}=h_{i,j}$

• Nos contornos de y: ${\tfrac {\partial h}{\partial y}}=0$, discretizando essa derivada temos que: $h_{i,j(+/-)1}=h_{i,j}$

No desenvolvimento do programa não foi conseguido alcançar os resultados esperados, pois o sistema não converge: a velocidade aumenta e diminui infinitamente, fazendo com que a altura da onda aumente indefinidamente. Mesmo assim, vamos apresentar o código escrito na linguagem Python, pois não conseguimos entender o erro.

Foi cosiderada uma superfície constante e, desta forma, suas derivadas são nulas e não interferem no cálculo.

#%% Bibliotecas 
import numpy as np

import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation
import matplotlib.patches as mpatches
from IPython.display import HTML, Image

#%% Parametros

L_xf = 10.0  # m
L_x0 = -L_xf

NX = 100

dx = (L_xf - L_x0) / NX


L_yf = 10.0  # m
L_y0 = -L_yf

NY = 100

dy = (L_yf - L_y0) / NY

N_INNER = (NX - 2) * (NY - 2)

g = 9.8 # m /s^2

# Tempo
dt = 0.002
Nt = 1000

time_interval = np.arange(0, Nt*dt, dt)

#%% Discretização do espaço x-y

x = np.linspace(L_x0, L_xf, NX-1)  
y = np.linspace(L_y0, L_yf, NY-1)
X, Y = np.meshgrid(x, y) 

#%% Condições iniciais, em distribuição gaussiana

sigma = 0.6
sigma_v = 0.6

h_2d = 0.8 * np.exp(-(((X)**2 / 2*(sigma**2)) + ((Y)**2 / 2*(sigma**2))))
u_2d = 0.1 * np.exp(-(((X)**2 / 2*(sigma_v**2)) + ((Y)**2 / 2*(sigma_v**2))))
v_2d = 0.1 * np.exp(-(((X)**2 / 2*(sigma_v**2)) + ((Y)**2 / 2*(sigma_v**2))))

#%% Vetores das variáveis anteriores e historico das variaveis

h_ant = np.copy(h_2d)
v_ant = np.copy(v_2d)
u_ant = np.copy(u_2d)

# Inicilização das listas para armazenar os valores

hist_h, hist_u, hist_v = [], [], []

Função para resolução das equações diferenciais com FTCS.

#%% Equação diferencial

# Fator de multiplicacao 
fator_x = (dt / (2*(dx)))
fator_y = (dt / (2*(dy)))

def resolve_pdes(h, vx, vy):

    """
    Função que retorna os valores de profundidade, velocidade em x e em y 
    no próximo tempo
    
    Parametros
    -----------
    h : array
                  profundidade no tempo t 
    vx : array
        velocidade em x no tempo t
      
    vy : array
        velocidade em y no tempo t
    
    Retorna
    -----------
    prox_h : array
                  profundidade no tempo t + dt
    prox_u : array
        velocidade em x no tempo t + dt
      
    prox_v : array
        velocidade em y no tempo t + dt
    
    """

    # Inicializa os vetores para armazenarem os resultados calculados
    prox_h = np.ones(shape = (np.shape(h)), dtype = np.float64)
    prox_u = np.ones(shape = (np.shape(vx)), dtype = np.float64)
    prox_v = np.ones(shape = (np.shape(vy)), dtype = np.float64)
    
        
    # Loop nos pontos discretizados 
    
    for i in range(1, NX - 1):
        for j in range(1, NY - 1):
        
            # Alturas e velocidades conforme a posicao do ponto:
            # _l : ponto a esquerda, _r: ponto a direita, _up: ponto acima, _d: abaixo
        
            # Condicao a esquerda ------------------
            if i == 1:  # primeiro x interno
                        
                h_l = h[i, j]
                u_l = -vx[i, j]
                # u_l = 0
                v_l = vy[i, j]
        
            else:
                h_l = h[i-1, j]
                u_l = vx[i-1, j]
                v_l = vy[i-1, j]
            # --------------------------------------
        
            # Condicao a direita -------------------
            if i == NX - 2:  # ultimo x interno

                h_r = h[i, j]
                u_r = -vx[i, j]
                # u_r = 0
                v_r = vy[i, j]
            
            else:
                h_r = h[i+1, j]
                u_r = vx[i+1, j]
                v_r = vy[i+1, j]
            # --------------------------------------
        
            # Condicao abaixo  ----------------------
            if j == 1:  # primeiro y interno 
                h_d = h[i, j]
                u_d = vx[i, j]
                v_d = - vy[i, j]
                # v_d = 0
            
            else:
                h_d = h[i, j - 1]
                v_d = vy[i, j - 1]
                u_d = vx[i, j - 1]
            # --------------------------------------
        
            # Condicao acima  ----------------------
            if j == NY - 2:  # utlimo y interno
            
            
                h_up = h[i, j]
                u_up = vx[i, j]
                # v_up = 0
                v_up = - vy[i, j]
            
            else:
                
                h_up = h[i, j + 1]
                v_up = vy[i, j + 1]
                u_up = vx[i, j + 1]
            # --------------------------------------


            ## Primeira Equação
        
            h_ij = h[i, j] - \
                  (h_r  * u_r  - h_l  * u_l) * fator_x - \
                  (h_up * v_up - h_d  * v_d) * fator_y 
            
            prox_h[i, j] = h_ij
        
            # ## Segunda equação
        
            hu_ij = (h[i, j] * vx[i, j]) - \
                    fator_x * (
                        ((h_r  * (u_r ** 2))  + (1/2 * g * (h_r ** 2))) -\
                        ((h_l  * (u_l ** 2))  + (1/2 * g * (h_l ** 2)))
                    ) - fator_y * (
                        (h_up * u_up * v_up ) - (h_d  * u_d  * v_d)
                    )
     
            prox_u[i, j] = hu_ij / h_ij

            # # ## Terceira Equação
        
            hv_ij = (h[i, j] * vy[i, j]) - \
                    fator_x * (
                        (h_r  * u_r * v_r) - (h_l  * u_l * v_l)
                    ) - \
                    fator_y * (
                        ((h_up * (v_up ** 2)) + (1/2 * g * (h_up ** 2))) - \
                        ((h_d  * (v_d  ** 2)) + (1/2 * g * (h_d  ** 2)))
                    )

        
            prox_v[i, j] = hv_ij / (h_ij)
    
    return prox_h, prox_u, prox_v

#%% Cálculo

# Resolve para cada passo de tempo

for t in time_interval:
    

    h, u, v = resolve_pdes(h_ant, u_ant, v_ant)  # valores das variaveis no tempo = t + dt
   
    
    # adicionar essas variáveis em listas pra conseguir plotar dps 
    hist_h.append(h); hist_u.append(u); hist_v.append(v)
    
      
    # Coloca as variaveis atuais como anteriores pro proximo calculo
    h_ant = np.copy(h)
    u_ant = np.copy(u)
    v_ant = np.copy(v)

#%% Gráfico animado

# Reorganiza os vetores para plotar

x_2d = X[0]
y_2d = Y[0]

for i in range(1,len(X)):
    
    x_2d = np.append(x_2d, X[i])
    y_2d = np.append(y_2d, Y[i])

def animate(i):
    plt.clf()  # limpa a figura, pra nao ficar sobrepondo figs
    
    # titulos
    plt.suptitle('Evolução da onda', fontsize=14)
    plt.title(f'Tempo: {round(dt*i, 3)}', fontsize=12)
    
    # plot
    plt.scatter(x_2d, y_2d, c=hist_h[0::8][i], marker='.')
    plt.colorbar()
    
    # axis


# fig = plt.figure(figsize=(16,9))
# ax = fig.gca(projection='3d')

fig, ax = plt.subplots()
ani = animation.FuncAnimation(fig, animate, frames = len(hist_h[0::8]), repeat=False, interval=0.1)

#%% Salvar o gif
ani.save('onda.gif', writer='imagemagick', fps=5)