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}}+\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 (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\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\qquad (35)$

${\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}}\qquad (35)$

${\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}}\qquad (36)$

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\qquad (36)$

$\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}}\qquad (37)$

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

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 (39)$

Definindo as quantidades:

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

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

$f(x,y,t)\equiv {\frac {gn^{2}}{D^{7/3}}}M(M^{2}+N^{2})^{1/2}\qquad (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}-\Delta t{\Bigg (}{\frac {M2_{i+1,j}^{n}-M2_{i-1,j}^{n}}{2\Delta x}}+{\frac {(MN)_{i,j+1}^{n}-(MN)_{i,j-1}^{n}}{2\Delta y}}+gD_{i,j}^{n}{\frac {\eta _{i+1,j}^{n}-\eta _{i-1,j}^{n}}{2\Delta x}}+f_{i,j}^{n}{\Bigg )}\qquad (43)$

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

$N_{i,j}^{n+1}=N_{i,j}^{n}-\Delta t{\Bigg (}{\frac {N2_{i,j+1}^{n}-M2_{i,j-1}^{n}}{2\Delta y}}+{\frac {(MN)_{i+1,j}^{n}-(MN)_{i-1,j}^{n}}{2\Delta x}}+gD_{i,j}^{n}{\frac {\eta _{i+1,j}^{n}-\eta _{i-1,j}^{n}}{2\Delta y}}+f_{i,j}^{n}{\Bigg )}\qquad (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.

${\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. O código foi 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.

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 estas 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()