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 Equações 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}+\mathrm{\nabla }.(\rho \mathbf{𝐮})=0(3)$

$\frac{D\mathbf{𝐮}}{Dt}+\frac{1}{\rho }\mathrm{\nabla }p+\frac{1}{\rho }\mathrm{\nabla }.\mathit{\tau }+\mathbf{𝐠}=0(4)$

$\rho$ é a densidade; p é a pressão; $\mathbf{𝐮}=(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{𝐠}$ é o vetor aceleração da gravidade; $\mathit{\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{𝐯}.\mathrm{\nabla }\eta =w$ , onde $z=\eta (x,y,t)(4)$

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

$\eta$ é o deslocamento vertical da água sobre a superfície em repouso, $\mathbf{𝐯}=(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.

$\mathrm{\nabla }.\mathbf{𝐮}=0(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:

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

Teorema de Leibniz:

$\frac{d}{dx}(\underset{a(x)}{\overset{b(x)}{\int }}f(x,t)dt)=f(x,b(x))\cdot \frac{d}{dx}b(x)-f(x,a(x))\cdot \frac{d}{dx}a(x)+\underset{a(x)}{\overset{b(x)}{\int }}\frac{\partial }{\partial x}f(x,t)dt(8)$

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

$\frac{\partial }{\partial x}\underset{-h}{\overset{\eta }{\int }}udz+\frac{\partial }{\partial y}\underset{-h}{\overset{\eta }{\int }}vdz-w{|}_{eta}-\mathbf{𝐯}.\mathrm{\nabla }\eta =0(9)$

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

$\frac{\partial }{\partial x}\underset{-h}{\overset{\eta }{\int }}udz+\frac{\partial }{\partial y}\underset{-h}{\overset{\eta }{\int }}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(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}\underset{-h}{\overset{\eta }{\int }}udz=uD(11)$

$N=\frac{\partial }{\partial y}\underset{-h}{\overset{\eta }{\int }}vdz=vD(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(10)$

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

$$

$$

$\frac{1}{\rho }\frac{\partial P}{\partial x}+g_z=0(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)(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(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(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(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}\left(\frac{M^2}{D}\right)+\frac{\partial }{\partial y}\left(\frac{MN}{D}\right)+gD\frac{\partial \eta }{\partial x}=0(20)$

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

Iremos escrever as equações das águas rasas considerando o tensor de estresse $\mathit{\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\ne j$, e as perpendiculares por elementos ${\tau }_{ij}$ onde $i=j$

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

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

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

$\frac{1}{\rho }(\frac{\partial }{\partial x}{\tau }_{zx}+\frac{\partial }{\partial y}{\tau }_{zy}+\frac{\partial }{\partial z}{\tau }_{zz})(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 (\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})(25)$

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

${\tau }_{yz}={\tau }_{zy}=\mu (\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z})(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:

$\underset{-h}{\overset{\eta }{\int }}\frac{1}{\rho }(\frac{\partial }{\partial x}{\tau }_{xx}+\frac{\partial }{\partial y}{\tau }_{xy}+\frac{\partial }{\partial z}{\tau }_{xz})=\frac{{\tau }_x}{\rho }-A\left(\frac{{\partial }^2}{\partial x^2}M\frac{{\partial }^2}{\partial x^2}M\right)(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 $$ é:

$\frac{{\tau }_x}{\rho }=\frac{fM}{2D^2}(M^2+N^2{)}^{1/2}(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{2g{n}^2}{D^{1/3}}(30)$

Substituindo (30) em (29) obtemos:

$\frac{{\tau }_x}{\rho }=\frac{g{n}^2}{D^{7/3}}M(M^2+N^2{)}^{1/2}(31)$

Generalizando a expressão (31) para a componente y. $\frac{{\tau }_y}{\rho }=\frac{g{n}^2}{D^{7/3}}N(M^2+N^2{)}^{1/2}(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}\left(\frac{M^2}{D}\right)+\frac{\partial }{\partial y}\left(\frac{MN}{D}\right)+gD\frac{\partial \eta }{\partial x}+\frac{g{n}^2}{D^{7/3}}M(M^2+N^2{)}^{1/2}=0(33)$

$\frac{\partial N}{\partial t}+\frac{\partial }{\partial y}\left(\frac{N^2}{D}\right)+\frac{\partial }{\partial x}\left(\frac{MN}{D}\right)+gD\frac{\partial \eta }{\partial x}+\frac{g{n}^2}{D^{7/3}}N(M^2+N^2{)}^{1/2}=0(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 $\overrightarrow{z}$
• A aceleração na direção da velocidade na direção $\overrightarrow{z}$ é zero
• As componentes das velocidades em $\overrightarrow{x}$ e em $\overrightarrow{y}$ ($\overrightarrow{u}$ e $\overrightarrow{v}$) não variam em $\overrightarrow{z}$

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

$\frac{\partial h}{\partial t}+\frac{\partial hu}{\partial x}+\frac{hv}{\partial y}=0(35)$

$\frac{\partial hu}{\partial t}+\frac{\partial (h{u}^2+\frac{1}{2}g{h}^2)}{\partial x}+\frac{huv}{\partial y}=-gh\frac{\partial b}{\partial x}(35)$

$\frac{\partial hv}{\partial t}+\frac{huv}{\partial x}+\frac{\partial (h{v}^2+\frac{1}{2}g{h}^2)}{\partial y}=-gh\frac{\partial b}{\partial y}(36)$

Onde $h$ é a altura do fluido desde a base, $\overrightarrow{u},\overrightarrow{v}$ são as velocidades médias na direções $\overrightarrow{x},\overrightarrow{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 $\mathrm{\Delta }{𝒪}(x^2)$ na região espacial, enquanto na temporal é de ordem $\mathrm{\Delta }{𝒪}(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(36)$

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

$n_{i,j}^{n+1}=n_{i,j}^n-\mathrm{\Delta }t(\frac{M_{i+1,j}^n-M_{i-1,j}^n}{2\mathrm{\Delta }x}+\frac{M_{i,j+1}^n-M_{i,j-1}^n}{2\mathrm{\Delta }y})(38)$

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

$\Rightarrow \frac{\partial M}{\partial t}=-\left(\frac{\partial }{\partial x}\right(\frac{M^2}{D})+\frac{\partial }{\partial y}(\frac{MN}{D})+gD\frac{\partial \eta }{\partial x}+\frac{g{n}^2}{D^{7/3}}M(M^2+N^2{)}^{1/2})(39)$

Definindo as quantidades:

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

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

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

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-\mathrm{\Delta }t(\frac{M{2}_{i+1,j}^n-M{2}_{i-1,j}^n}{2\mathrm{\Delta }x}+\frac{(MN{)}_{i,j+1}^n-(MN{)}_{i,j-1}^n}{2\mathrm{\Delta }y}+g{D}_{i,j}^n\frac{{\eta }_{i+1,j}^n-{\eta }_{i-1,j}^n}{2\mathrm{\Delta }x}+f_{i,j}^n)(43)$

Generalizando a expressão (43) para o fluxo de descarga N temos:

$N_{i,j}^{n+1}=N_{i,j}^n-\mathrm{\Delta }t(\frac{N{2}_{i,j+1}^n-M{2}_{i,j-1}^n}{2\mathrm{\Delta }y}+\frac{(MN{)}_{i+1,j}^n-(MN{)}_{i-1,j}^n}{2\mathrm{\Delta }x}+g{D}_{i,j}^n\frac{{\eta }_{i+1,j}^n-{\eta }_{i-1,j}^n}{2\mathrm{\Delta }y}+f_{i,j}^n)(44)$

Simulações Computacionais de Tsunamis

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.

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

$\frac{hu{)}_{i,j}^{t+\mathrm{\Delta }t}-(hu{)}_{i,j}^t}{\mathrm{\Delta }t}+\left[\frac{(h{u}^2+\frac{1}{2}g{h}^2{)}_{i+1,j}^t-(h{u}^2+\frac{1}{2}g{h}^2{)}_{i-1,j}^t}{2\mathrm{\Delta }x}\right]+\left[\frac{(huv{)}_{i,j+1}^t-(huv{)}_{i,j-1}^t}{2\mathrm{\Delta }y}\right]=-g{h}_{i,j}^t{b}_{x.i,j}$

$\frac{(hv{)}_{i,j}^{t+\mathrm{\Delta }t}-(hv{)}_{i,j}^t}{\mathrm{\Delta }t}+\left[\frac{(huv{)}_{i+1,j}^t-(huv{)}_{i-1,j}^t}{2\mathrm{\Delta }x}\right]+\left[\frac{(h{v}^2+1/2g{h}^2{)}_{i,j+1}^t-(h{v}^2+1/2g{h}^2{)}_{i,j-1}^t}{2\mathrm{\Delta }y}\right]=-g{h}_{i,j}^t{b}_{y.i,j}$

Para os contornos foi utilizado que:

• $u(x_f,y)=-u(x_f-\mathrm{\Delta }x,y)$

• $u(x_i,y)=-u(x_i+\mathrm{\Delta }x,y)$

• $v(x,y_f)=-v(x,y_f-\mathrm{\Delta }y)$

• $v(x,y_i)=-v(x,y_i+\mathrm{\Delta }y)$

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

• Nos contornos de y: $\frac{\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)

Forma dissipativa 2D

Os exemplos que seguem utilizam as equações de ondas rasas (38),(43) e (44) para calcular os passos de tempo de $\eta$, $M$, $N$, onde as funções em python atualiza_eta, atualiza_M, e atualiza_N implementam computacionalmente isto. Para implementar estás funções e outras ideias do nosso programa, o seguinte código fonte da referência foi usado como base.

def atualiza_eta(eta, M, N, dx, dy, dt, nx, ny):
    
    for j in range(1,nx-1):
        for i in range(1,ny-1):
            
            dMdx = (M[j+1,i] - M[j-1,i]) / (2. * dx)
            dNdy = (N[j,i+1] - N[j,i-1]) / (2. * dy)
            eta[j, i] = eta[j, i] - dt * (dMdx + dNdy)
    
    #Condições de contorno do problema
    eta[0,:] = eta[1,:]
    eta[-1,:] = eta[-2,:]
    eta[:,0] = eta[:,1]
    eta[:,-1] = eta[:,-2]
    
    return eta

def atualiza_M(eta, M, N, D, g, h, n, dx, dy, dt, nx, ny):
    
    M2 = M **2 / D
    MN = M * N / D
    fric = g * n**2 * M * np.sqrt(M**2 + N**2) / D**(7./3.)
    
    for j in range(1,nx-1):
        for i in range(1,ny-1):            
            
            dM2dx = (M2[j+1,i] - M2[j-1,i]) / (2. * dx)
            dMNdy = (MN[j,i+1] - MN[j,i-1]) / (2. * dy)
            dETAdx = (eta[j+1,i] - eta[j-1,i]) / (2. * dx)
            
            M[j, i] = M[j, i] - dt * (dM2dx + dMNdy + g * D[j,i] * dETAdx + fric[j,i])
            
    return M

def atualiza_N(eta, M, N, D, g, h, n, dx, dy, dt, nx, ny):
    
    
    MN = M * N / D
    N2 = N**2 / D
    fric = g * n**2 * N * np.sqrt(M**2 + N**2) / D**(7./3.)
    
    for j in range(1,nx-1):
        for i in range(1,ny-1):            
            
            dMNdx = (MN[j+1,i] - MN[j-1,i]) / (2. * dx)
            dN2dy = (N2[j,i+1] - N2[j,i-1]) / (2. * dy)
            dETAdy = (eta[j,i+1] - eta[j,i-1]) / (2. * dy)
            
            N[j, i] = N[j, i] - dt * (dMNdx + dN2dy + g * D[j,i] * dETAdy + fric[j,i])
                        
    return N

A função shallow water waves recebe os parâmetros iniciais do nosso programa, executa o Loop responsável pela atualização das variáveis da amplitude da onda e do fluxo de descarga com o tempo, através da chamada das funções atualiza M,N e eta. Posteriormente, a cada passagem dentro do loop um plot do sistema é feito. Obs: não colocamos todo código da função shallow_water na imagem a seguir, apenas a que mencionamos neste parágrafo.

def shallow_water(eta0, M0, N0, h, g, n, nt, dx, dy, dt, X, Y):
     
    eta = eta0.copy()
    M = M0.copy()
    N = N0.copy()
    
    D = eta + h

    # ... 
    
    for k in range(1,nt):
        

        eta = atualiza_eta(eta, M, N, dx, dy, dt, nx, ny)
        M = atualiza_M(eta, M, N, D, g, h, n, dx, dy, dt, nx, ny)
        N = atualiza_N(eta, M, N, D, g, h, n, dx, dy, dt, nx, ny)

        D = eta + h


        fig = plt.figure(figsize=(8.,6.))

        fundo = plt.imshow(-h, 'Purples', interpolation = 'nearest', extent = limites)
        amp = plt.imshow(eta, extent = limites, interpolation = 'sinc', cmap = 'seismic', alpha= 0.75, vmin=-0.4, vmax= 0.4)
    
        #plt.title('tempo = %f', dt*n )
        #plt.plot(f'Tempo {round(k*dt,3)} s')
        plt.xlabel('x [m]')
        plt.ylabel('y [m]')
        cbar_amp = plt.colorbar(amp)
        cbar_fundo = plt.colorbar(fundo)
        cbar_fundo.set_label(r'$-h$ [m]')
        cbar_amp.set_label(r'$\eta$ [m]')

        plt.show()