Equação de Águas Rasas: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
 
(191 revisões intermediárias por 2 usuários não estão sendo mostradas)
Linha 1: Linha 1:
'''Em construção'''
'''Grupo: Gabriel Schmökel e Julia Remus'''
'''Grupo: Gabriel Schmökel e Julia Remus'''


O objetivo deste trabalho é buscar a solução das equações de águas rasas, por meio de métodos de integração numérica, para resolução de equações diferenciais parciais (EDP's) e posteriormente apresentar uma breve interpretação física das soluções. Demonstramos, nesta página, a derivação das equações, junto com a explicação de cada quantidade física presente. A discretização das equações de águas rasas conservativas e não conservativas são feitas por FTCS explícito, também foi aplicado FTCS implícito para as equações em uma forma simplificada, representada pela equação da deriva. Os resultados obtidos para equação simplificada são comparados com a solução exata, e exemplos mais complexos são solucionados para as formas conservativa e não conservativa. Ao final, uma comparação é feita entre os resultados das equações conservativas e dissipativas. 


== Introdução ==
== 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.  
As equações de águas rasas têm aplicações físicas na previsão de tsunâmis, em fluxos atmosféricos, ondas de tempestade e fluxos planetários. Na descrição física dos problemas de fluxos de fluído em ondas, as equações de águas rasas em uma dimensão são dadas por:
 
 
==== 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.
 
[[Arquivo:imagem1.jpeg|250px|thumb|center]]
 
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. 


[[Arquivo:imagem2.jpeg|250px|thumb|center]]
<math>\dfrac{\partial \eta}{\partial t} + \dfrac{\partial u D}{\partial x} = 0 \qquad (1) </math>


II. A onda gerada se propaga ao longo de todas as direções do plano da água.
<math>\dfrac{\partial \eta u}{\partial t} + \dfrac{\partial \left (\eta u^2 + \dfrac{1}{2}g \eta^2 \right)}{\partial x} = -g \eta \dfrac{\partial h}{\partial x} \qquad (2) </math>


[[Arquivo:imagem3.jpeg|250px|thumb|center]]
As componentes da equação de águas rasas podem ser melhor interpretadas através da seguinte figura:


III. A medida que a onda se aproxima da superfície ela diminui sua velocidade e aumenta sua amplitude
[[Arquivo:TSUNAMI3-06.png|600px|center|Componentes das equações de águas rasas]]


[[Arquivo:imagem4.jpeg|250px|thumb|center]]
<math> \eta = \eta(x,t) </math> corresponde a amplitude da onda, <math> h=h(x) </math> determina a profundidade do mar em repouso, <math> D</math> é o deslocamento total da água, <math> u</math> é a velocidade do fluído. Resolvendo a EDP da equação de águas rasas, obtemos como a amplitude da onda se comporta ao longo do tempo e do espaço.
 
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 ==
== Teoria ==


===Derivação das Equações de Águas Rasas===
===Derivação das Equações de Águas Rasas===
Iremos demonstrar como chegamos nas equações de águas rasas em duas dimensões, nas formas conservativa e dissipativa, em representações do fluxo de descarga e de velocidade. Posteriormente, tendo as equações em 2D, iremos simplificar elas para a forma unidimensional. Neste processo de demonstração, iremos explicar a interpretação física de cada quantidade presente nas equações.


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:
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:


<math> \frac{d\rho}{dt} +\nabla . (\rho \mathbf{u}) = 0 \qquad (3) </math>
<math> \nabla . (\rho \mathbf{u}) = \frac{d\rho}{dt} \qquad (3) </math>


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


A equação da continuidade em (3) descreve o balanço de massas para os elementos de volume infinitesimais que pertencem ao fluído, onde a quantidade do lado esquerdo da equação informa o fluxo de massa que entra e sai pelo elemento de volume, e a quantidade do lado direito está relacionada com a massa que se acumula ao longo do tempo <ref>Equação da continuidade mássica: balanços de massa diferenciais. Bloom Consultoria.Disponível em: <https://www.youtube.com/watch?v=pEip-GvO0LM&list=PL1yqHjPQz-Lqjri07DqZ3RsSWJfICvdiu&index=3></ref>. Nesta expressão <math> \rho  </math> é a densidade, e <math> \mathbf u=(u,v,w,t) </math> é o campo de velocidades, onde u,v e w são as velocidades das partículas que compõe o fluído nas direções x,y,z.


<math> \rho  </math> é a densidade; p é a pressão; <math> \mathbf u=(u,v,w) </math> é 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; <math> \mathbf{g} </math> é o vetor aceleração da gravidade; <math> \boldsymbol{\tau} </math> é o tensor tensão, onde as componentes deste tensor são as tensões normais e tangenciais de cisalhamento, expressas por <math> \tau_{ij} </math>, no qual <math> i </math> indica a direção e <math> j </math> o plano normal.  
As equações de Navier-Stokes em (4) são balanços diferenciais da quantidade de movimento, obtidas através da aplicação da segunda lei de Newton em cada ponto do escoamento <ref>Equação de Navier-Stokes (Parte 1) - Derivadas materiais.
Bloom Consultoria.Disponível em: https://www.youtube.com/watch?v=FLoZODPpayM</ref> <ref>Equação de Navier-Stokes (Parte 2) - Equação diferencial da quantidade de movimento. Bloom Consultoria.Disponível em: https://www.youtube.com/watch?v=e06ZRzdO4iM</ref> <ref>Equação de Navier-Stokes (Parte 3) - Tensões Normais e Cisalhantes. Bloom Consultoria.Disponível em: https://www.youtube.com/watch?v=na2kGOSYNv8</ref>.


Introduzindo as condições de contorno <ref>SEGUR, Harvey; YAMAMOTO, Hiroki. Lecture 8: The Shallow-Water Equations.Disponível em: <https://docplayer.net/49487265-Lecture-8-the-shallow-water-equations.html></ref> para a superfície <math> z(x,y,t) </math> e para a profundidade do oceano <math> h(x,y) </math>:
* <math> Du/Dt </math> é a aceleração da partícula fluída ao longo do campo de velocidade <math> \mathbf u=(u,v,w,t) </math>.
* <math>  \nabla . \boldsymbol{\tau} / \rho</math> está associado as tensões tangenciais e normais atuando sobre os elementos de volume (<math> \boldsymbol{\tau} </math> é o tensor tensão, as componentes deste tensor são as tensões normais e tangenciais de cisalhamento, expressas por <math> \tau_{ij} </math>, no qual <math> i </math> indica a direção e <math> j </math> o plano normal).
* <math> \nabla p / \rho </math> está associado as pressões que atuam sobre os elementos do fluído.
* <math> \mathbf{g} </math> é o vetor aceleração da gravidade atuando sobre os elementos infinitesimais de volume do fluído.   


<math> \frac{D \eta}{Dt} = \frac{\partial \eta}{\partial t} +\mathbf{v} . \nabla \eta = w </math> , onde <math> z= \eta(x,y,t) \qquad (4) </math>
Introduzindo as condições de contorno para a superfície <math> z(x,y,t) </math> e para a profundidade do oceano <math> h(x,y) </math> <ref>SEGUR, Harvey; YAMAMOTO, Hiroki. Lecture 8: The Shallow-Water Equations.Disponível em: <https://docplayer.net/49487265-Lecture-8-the-shallow-water-equations.html></ref>:


<math> \mathbf{u} . \nabla (z + h(x,y)) = </math> , onde <math> z =-h(x,y) \qquad (5)</math>
<math> \frac{D \eta}{Dt} = \frac{\partial \eta}{\partial t} +\mathbf{v} . \nabla \eta = w </math> , onde <math> z= \eta(x,y,t) \qquad (5) </math>


<math> \mathbf{u} . \nabla  (z + h(x,y)) = 0  </math> , onde <math> z =-h(x,y) \qquad (6)</math>


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


A equação da continuidade em (3) pode ser simplificada, que a densidade do fluído no oceano <math> \rho </math> não varia significativamente com o tempo e a posição.
A equação da continuidade em (3) pode ser simplificada pelo fato do fluído ser incompressível, isto implica que a densidade <math> \rho </math> é constante.


<math> \nabla . \mathbf{u} = 0 \qquad (6) </math>  
<math> \nabla . \mathbf{u} = 0 \qquad (7) </math>  


Integrando a expressão da continuidade em (6), utilizando a regra da integral de Leibniz <ref name=Hopf>https://en.wikipedia.org/wiki/Leibniz_integral_rule</ref>, com os limites indo de <math> -h(x,y) </math> até <math> \eta (x,y,t) </math> chegamos na seguinte expressão:
Integrando a expressão da continuidade em (7), utilizando a regra da integral de Leibniz <ref>Leibniz integral rule. Disponível em: https://en.wikipedia.org/wiki/Leibniz_integral_rule </ref>, com os limites indo de <math> -h(x,y) </math> até <math> \eta (x,y,t) </math> chegamos na seguinte expressão:


<math> \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 =
<math> \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} u dz + \frac{\partial}{\partial y} \int_{-h}^{\eta} v dz +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) </math>  
         \frac{\partial}{\partial x} \int_{-h}^{\eta} u dz + \frac{\partial}{\partial y} \int_{-h}^{\eta} v dz +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 (8) </math>  


Teorema de Leibniz:
Teorema de Leibniz:


<math>\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)</math>
<math>\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 (9)</math>


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


<math> \frac{\partial}{\partial x} \int_{-h}^{\eta} u dz + \frac{\partial}{\partial y} \int_{-h}^{\eta} v dz - w \Big |_{eta} -\mathbf{v} . \nabla \eta = 0 \qquad (9) </math>
<math> \frac{\partial}{\partial x} \int_{-h}^{\eta} u dz + \frac{\partial}{\partial y} \int_{-h}^{\eta} v dz - w \Big |_{eta} -\mathbf{v} . \nabla \eta = 0 \qquad (10) </math>


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


<math>  \frac{\partial}{\partial x} \int_{-h}^{\eta} u dz + \frac{\partial}{\partial y} \int_{-h}^{\eta} v dz + \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} </math>
<math>  \frac{\partial}{\partial x} \int_{-h}^{\eta} u dz + \frac{\partial}{\partial y} \int_{-h}^{\eta} v dz + \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} </math>


<math> \Rightarrow \frac{\partial \eta}{\partial t} + \frac{\partial uD}{\partial x}+ \frac{\partial vD}{\partial y} = 0 \qquad (10) </math>
<math> \Rightarrow \frac{\partial \eta}{\partial t} + \frac{\partial uD}{\partial x}+ \frac{\partial vD}{\partial y} = 0 \qquad (11) </math>
 
(11) é a primeira das equações das águas rasas que obtemos, onde <math> D </math> é o comprimento da água total do fundo do oceano até a amplitude da onda.
Podemos expressar (11) através do fluxo de descarga nas direções x e y, estas quantidades estão relacionadas com as velocidades da seguinte forma <ref name=MANUAL> IMAMURA, Fumihiko.Tsunami Modelling Manual.Disponível em: http://www.tsunami.civil.tohoku.ac.jp/hokusai3/J/projects/manual-ver-3.1.pdf </ref>:


(10) é a primeira das equações das águas rasas que obtemos, onde <math> D </math> é o comprimento da água total do fundo do oceano até a amplitude da onda.
<math> M = \frac{\partial}{\partial x} \int_{-h}^{\eta} u dz = uD \qquad (12) </math>
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 <ref name=Hopf>http://www.tsunami.civil.tohoku.ac.jp/hokusai3/J/projects/manual-ver-3.1.pdf</ref>:


<math> M = \frac{\partial}{\partial x} \int_{-h}^{\eta} u dz = uD \qquad (11) </math>
<math> N = \frac{\partial}{\partial y} \int_{-h}^{\eta} v dz = vD \qquad (13) </math>


<math> N = \frac{\partial}{\partial y} \int_{-h}^{\eta} v dz = vD \qquad (12) </math>
Substituindo (12) e (13) em (11) chegamos na representação do fluxo de descarga para uma das equações de águas rasas.


Substituindo (11) e (12) em (10) chegamos na representação do fluxo de descarga para uma das equações de águas rasas.
<math> \Rightarrow \frac{\partial \eta}{\partial t} + \frac{\partial M}{\partial x}+ \frac{\partial N}{\partial y} = 0 \qquad (14) </math>


<math> \Rightarrow \frac{\partial \eta}{\partial t} + \frac{\partial M}{\partial x}+ \frac{\partial N}{\partial y} = 0 \qquad (10) </math>
Conhecendo as taxas dos fluxos de descarga em relação as regiões espaciais, podemos determinar a taxa da variação da amplitude da onda em relação ao tempo.


Escrevendo as quantidades de movimento de Navier-Stokes nas componentes x,y e z:
Vamos buscar obter as outras duas equações de águas rasas restantes, a partir das quantidades de movimento de Navier-Stokes. Nas componentes x,y e z temos:


<math> \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y}
<math> \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) </math>
         + w\frac{\partial u}{\partial z} +\frac{1}{\rho}\frac{\partial P}{\partial x} +g_x = 0 \qquad (15) </math>


<math> \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y}
<math> \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) </math>
         + w\frac{\partial v}{\partial z} +\frac{1}{\rho}\frac{\partial P}{\partial x} +g_y = 0 \qquad (16) </math>


<math> \frac{1}{\rho}\frac{\partial P}{\partial x} +g_z = 0 \qquad (16) </math>
<math> \frac{1}{\rho}\frac{\partial P}{\partial x} +g_z = 0 \qquad (17) </math>


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 <math> g_x</math> e <math> g_y</math> em (14) e passamos a definir <math> g_z = g </math>.
Na componente z em (17) não consideramos a aceleração das partículas, pois a aceleração da gravidade é muito maior. Também tomamos como nulo as componentes <math> g_x</math> e <math> g_y</math> em (15) e (16), assim passamos a definir <math> g_z = g </math>. Neste momento estamos desconsiderando as forças de fricção, por isso o tensor tensão também é nulo.


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


<math> \partial P = \rho g \partial z \Rightarrow P = \rho g (\eta - z) \qquad (17) </math>
<math> \partial P = \rho g \partial z \Rightarrow P = \rho g (\eta - z) \qquad (18) </math>


Substituindo a pressão em (14):
Substituindo a pressão em (15):


<math> \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y}
<math> \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) </math>
         + w\frac{\partial u}{\partial z} +g \frac{\partial \eta}{\partial x} =0 \qquad (19) </math>
 
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 <math>-h(x,y)</math> até <math>\eta(x,y,t)</math> chegamos em outra das equações de águas rasas:
Integrando a expressão (19), utilizando a regra da integral de Leibniz e as condições de contorno (5) e (6), com os limites indo de <math>-h(x,y)</math> até <math>\eta(x,y,t)</math> chegamos em outra das equações de águas rasas:


<math> \frac{\partial uD}{\partial t} + \frac{\partial u^{2}D}{\partial x} + \frac{\partial uvD}{\partial y}
<math> \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) </math>
         + \frac{g}{2} \frac{\partial D^2}{\partial x} =0 \qquad (20) </math>


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


<math> \frac{\partial vD}{\partial t} + \frac{\partial v^{2}D}{\partial y} + \frac{\partial uvD}{\partial x}
<math> \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) </math>
         + \frac{g}{2} \frac{\partial D^2}{\partial y} =0 \qquad (21) </math>


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


<math> \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) </math>
<math> \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 (22) </math>


<math> \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) </math>
<math> \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 (23) </math>


Iremos escrever as equações das águas rasas considerando o tensor de estresse <math> \boldsymbol{\tau} </math>. 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 <math> \tau_{ij}</math> onde <math> i \ne j </math>, e as perpendiculares por elementos <math> \tau_{ij}</math> onde <math> i = j </math>
Com as equações de águas rasas (20), (21), (22) e (23) podemos calcular as taxas de variação dos fluxos de descarga em relação ao tempo. As equações (11), (20) e (21), na representação de velocidades, e as equações (14), (22) e (23), na representação do fluxo de descargas, são as equações de águas rasas conservativas.


Decompondo nas componentes x,y, e z de <math> \frac{1}{\rho} \nabla .  \boldsymbol{\tau} </math> presente em (4):
Iremos buscar pelas equações de águas rasas não conservativas considerando o tensor de estresse <math> \boldsymbol{\tau} </math> em (4). 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 <math> \tau_{ij}</math> onde <math> i \ne j </math>, e as perpendiculares por elementos <math> \tau_{ij}</math> onde <math> i = j </math>.


<math> \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) </math>
Decompondo nas componentes x,y, e z de <math> \frac{1}{\rho} \nabla .  \boldsymbol{\tau} </math> presente em (4) temos:


<math> \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) </math>
<math> \frac{1}{\rho}\Big(\frac{\partial}{\partial x} \tau_{xx} + \frac{\partial}{\partial y} \tau_{xy} + \frac{\partial}{\partial z}\tau_{xz} \Big) \qquad (24) </math>


<math> \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) </math>
<math> \frac{1}{\rho}\Big(\frac{\partial}{\partial x} \tau_{yx} + \frac{\partial}{\partial y} \tau_{yy} + \frac{\partial}{\partial z}\tau_{yz} \Big) \qquad (25) </math>
 
<math> \frac{1}{\rho}\Big(\frac{\partial}{\partial x} \tau_{zx} + \frac{\partial}{\partial y} \tau_{zy} + \frac{\partial}{\partial z}\tau_{zz} \Big) \qquad (26) </math>


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.
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.


<math> \tau_{xy} = \tau_{yx} = \mu \Big( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \Big) \qquad (25) </math>
<math> \tau_{xy} = \tau_{yx} = \mu \Big( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \Big) \qquad (27) </math>


<math> \tau_{xz} = \tau_{zx} = \mu \Big( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \Big) \qquad (26) </math>
<math> \tau_{xz} = \tau_{zx} = \mu \Big( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \Big) \qquad (28) </math>


<math> \tau_{yz} = \tau_{zy} = \mu \Big( \frac{\partial w}{\partial y} + \frac{\partial v}{\partial z} \Big) \qquad (27) </math>
<math> \tau_{yz} = \tau_{zy} = \mu \Big( \frac{\partial w}{\partial y} + \frac{\partial v}{\partial z} \Big) \qquad (29) </math>


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:
Substituindo (27),(28) em (24), integrando em relação a componente z, utilizando a regra de Leibnz e as condições de contorno (5) e (6), obtemos:


<math> \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) </math>
<math> \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 (30) </math>


Onde <math> A </math> é a constante de viscosidade turbulenta, <math> \tau_x </math> é 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. <ref name=Hopf>http://www.tsunami.civil.tohoku.ac.jp/hokusai3/J/projects/manual-ver-3.1.pdf</ref>.
Onde <math> A </math> é a constante de viscosidade turbulenta, <math> \tau_x </math> é uma força de resistividade relativa ao movimento do fluído com o fundo do oceano na direção x. Podemos desconsiderar a constante de turbulência na situação em que não temos inclinações abrutas no fundo  do mar <ref name=MANUAL> IMAMURA, Fumihiko.Tsunami Modelling Manual.Disponível em: http://www.tsunami.civil.tohoku.ac.jp/hokusai3/J/projects/manual-ver-3.1.pdf </ref>.


Considerando que o fluído é uniforme, então a expressão para <math> \frac{\tau_x}{\rho} é </math> é:
Considerando que o fluído é uniforme, então a expressão para <math> \frac{\tau_x}{\rho} é </math> é:


<math> \frac{\tau_x}{\rho} = \frac{fM}{2D^{2}}(M^{2}+N^{2})^{1/2} \qquad (29) </math>  
<math> \frac{\tau_x}{\rho} = \frac{fM}{2D^{2}}(M^{2}+N^{2})^{1/2} \qquad (31) </math>  


<math> f </math> é o coeficiente de fricção, porém o coeficiente de rugosidade de Manning <math> n </math> é mais usado, alguns valores deste coeficiente são:
Onde <math> f </math> é o coeficiente de fricção, porém o coeficiente de rugosidade de Manning <math> n </math> é mais usado, alguns valores deste coeficiente são
<ref> LINSLEY, Ray K.; FRANZINI, Joseph B. Water Resources Engineering. </ref>:


*Cimento puro e metal liso <math> n = 0,010 </math>
{| border="4" cellpadding="2"
*Terra lisa               <math> n = 0,017 </math>
! Material
*Pedras, ervas daninhas   <math> n = 0,035 </math>
! Coeficiente de Rugosidade de Manning <math> (n) </math>
*Péssimo relevo de canal   <math> n = 0,060 </math>
|-
*Bom relevo de canal       <math> n = 0,025 </math>
| Cimento puro e metal liso || 0,010  
|-
| Terra lisa || 0,017
|-
| Pedras, ervas daninhas || 0,035
|-
| Péssimo relevo de canal || 0,060  
|-
| Bom relevo de canal || 0,025  
|}


O coeficiente de fricção <math> f </math> e o de rugosidade de Meanning <math> n </math> estão relacionados por:
O coeficiente de fricção <math> f </math> e o de rugosidade de Meanning <math> n </math> estão relacionados por:


<math> f = \frac{2gn^{2}}{D^{1/3}} \qquad (30) </math>
<math> f = \frac{2gn^{2}}{D^{1/3}} \qquad (32) </math>


Substituindo (30) em (29) obtemos:  
Substituindo (32) em (31) obtemos:  


<math> \frac{\tau_x}{\rho} = \frac{gn^{2}}{D^{7/3}} M(M^2 +N^2)^{1/2} \qquad (31) </math>
<math> \frac{\tau_x}{\rho} = \frac{gn^{2}}{D^{7/3}} M(M^2 +N^2)^{1/2} \qquad (33) </math>


Generalizando a expressão (31) para a componente y.
Generalizando a expressão (29) para a componente y.
<math> \frac{\tau_y}{\rho} = \frac{gn^{2}}{D^{7/3}} N(M^2 +N^2)^{1/2} \qquad (32) </math>


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.
<math> \frac{\tau_y}{\rho} = \frac{gn^{2}}{D^{7/3}} N(M^2 +N^2)^{1/2} \qquad (34) </math>


<math> \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) </math>
Adicionando, repectivamente, (33) e (34) nas expressões (22) e (24), obtemos as equações de águas rasas considerando as forças de fricção do fundo do oceano.


<math> \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) </math>
<math> \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 (35) </math>


=== Forma Conservativa ===
<math> \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 (36) </math>
Um modelo mais simples - desconsiderando a fricção, viscosidade do líquido e as forças de Coriolis sobre ele - pode ser obtido <ref>GARCÍA-NAVARRO, P; et al. The shallow water equations: An example of hyperbolic system. Espanha: 2008. Disponível em: <https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.571.1364&rep=rep1&type=pdf></ref><ref>KUHBACHER, Christian. Shallow Water: Derivation and Applications. Disponível em: <http://www.mathematik.tu-dortmund.de/lsiii/cms/papers/Kuehbacher2009.pdf></ref>. Para desenvolvê-lo são necessárias algumas premissas:


* O comprimento da onda é muito maior que as contribuições na direção <math>\vec{z}</math>
As equações (14), (35) e (36), na representação do fluxo de descargas, são as equações de águas rasas não conservativas.
* A aceleração na direção da velocidade na direção <math>\vec{z}</math> é zero
* As componentes das velocidades em <math>\vec{x}</math> e em <math>\vec{y}</math> (<math>\vec{u}</math> e <math>\vec{v}</math>) não variam em <math>\vec{z}</math>


Em uma dimensão podemos expressar as equações de águas rasas eliminando a componente <math> y</math> das expressões (14),(35) e tomando o fluxo de descarga <math> N </math> como nulo.


O sistema então pode ser descrito pelas seguintes equações:
<math>\frac{\partial \eta}{\partial t} + \frac{\partial M}{\partial x} = 0 \qquad (37) </math>


<math>\dfrac{\partial h}{\partial t} + \dfrac{\partial hu}{\partial x} + \dfrac{hv}{\partial y} = 0 \qquad (35) </math>
<math> \frac{\partial M}{\partial t} + \frac{\partial }{\partial x}\Big(\frac{M^{2}}{D}\Big) + gD \frac{\partial \eta}{\partial x} +\frac{gn^{2}}{D^{7/3}} M^2= 0 \qquad (38) </math>


<math>\dfrac{\partial hu}{\partial t} + \dfrac{\partial \left ( hu^2 + \dfrac{1}{2}g h^2 \right)}{\partial x} + \dfrac{huv}{\partial y} = -gh\dfrac{\partial b}{\partial x} \qquad (35) </math>
===Simplificação das Equações de Águas Rasas===


<math>\dfrac{\partial hv}{\partial t} + \dfrac{huv}{\partial x} + \dfrac{\partial \left ( hv^2 + \dfrac{1}{2}g h^2 \right)}{\partial y}= -gh\dfrac{\partial b}{\partial y}\qquad (36) </math>
As equações de águas rasas podem ser simplificadas para equação de advecção através das seguintes considerações:


Onde <math>h</math> é a altura do fluido desde a base, <math>\vec{u}, \vec{v}</math> são as velocidades médias na direções <math>\vec{x}, \vec{y}</math>, <math>g</math> é a constante gravitacional e <math>b(x, y)</math> é função que descreve a superfície onde acontece o movimento <ref name=Hopf>http://www.tsunami.civil.tohoku.ac.jp/hokusai3/J/projects/manual-ver-3.1.pdf</ref>.
* A velocidade do fluído é constante.
* A profundidade do fundo do oceano é constante.  


=== Forma dissipativa ===
De (39) e das considerações acima temos:


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 <math> \Delta \mathcal{O}(x^2) </math> na região espacial, enquanto na temporal é de ordem <math> \Delta \mathcal{O}(t^1) </math>. O método é conhecido como ''leap-frog method'' devido a discretização central na região espacial.
<math>\frac{\partial \eta}{\partial t} + \frac{\partial M}{\partial x} = \frac{\partial \eta}{\partial t} + \frac{\partial uD}{\partial x} = \frac{\partial \eta}{\partial t} +u\frac{\partial (\eta + h)}{\partial x} = \frac{\partial \eta}{\partial t} + u\frac{\partial \eta}{\partial x} = 0 </math>


Discretizando a expressão (10) pelo ''leap-frog method'':
<math> \Rightarrow \frac{\partial \eta}{\partial t} + u\frac{\partial \eta}{\partial x} = 0 \qquad (39) </math>


<math> \frac{\partial \eta}{\partial t} + \frac{\partial uD}{\partial x}+ \frac{\partial vD}{\partial y} = 0 \qquad (36) </math>
=== Discretização na 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 <ref>GARCÍA-NAVARRO, P; et al. The shallow water equations: An example of hyperbolic system. Espanha: 2008. Disponível em: <https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.571.1364&rep=rep1&type=pdf></ref><ref>KUHBACHER, Christian. Shallow Water: Derivation and Applications. Disponível em: <http://www.mathematik.tu-dortmund.de/lsiii/cms/papers/Kuehbacher2009.pdf></ref>. Para desenvolvê-lo são necessárias algumas premissas:


<math> \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) </math>
* O comprimento da onda é muito maior que as contribuições na direção <math>\vec{z}</math>
* A aceleração na direção da velocidade na direção <math>\vec{z}</math> é zero
* As componentes das velocidades em <math>\vec{x}</math> e em <math>\vec{y}</math> (<math>\vec{u}</math> e <math>\vec{v}</math>) não variam em <math>\vec{z}</math>  


<math> 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) </math>


Discretizando a expressão (33) pelo ''leap-frog method'':
O sistema então pode ser descrito pelas seguintes equações:  


<math> \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) </math>
<math>\dfrac{\partial D}{\partial t} + \dfrac{\partial Du}{\partial x} + \dfrac{\partial Dv}{\partial y} = 0 \qquad (40) </math>


Definindo as quantidades:
<math>\dfrac{\partial Du}{\partial t} + \dfrac{\partial \left (\partial Du^2 + \dfrac{1}{2}g D^2 \right)}{\partial x} + \dfrac{Duv}{\partial y} = -gD\dfrac{\partial b}{\partial x} \qquad (41) </math>


<math> M2 \equiv \frac{M^2(x,y,t)}{D(x,y,t)}\qquad (40) </math>
<math>\dfrac{\partial Dv}{\partial t} + \dfrac{Duv}{\partial x} + \dfrac{\partial \left ( Dv^2 + \dfrac{1}{2}g D^2 \right)}{\partial y}= -gD\dfrac{\partial b}{\partial y}\qquad (42) </math>


<math> MN \equiv \frac{M(x,y,t)N(x,y,t)}{D(x,y,t)} \qquad (41) </math>
Onde <math>D</math> é a altura do fluido desde a base, <math>\vec{u}, \vec{v}</math> são as velocidades médias na direções <math>\vec{x}, \vec{y}</math>, <math>g</math> é a constante gravitacional e <math>b(x, y)</math> é função que descreve a superfície onde acontece o movimento.


<math> f(x,y,t) \equiv \frac{gn^{2}}{D^{7/3}} M(M^2 +N^2)^{1/2} \qquad (42) </math>
Para descrever numericamente as equações de águas rasas na forma conservativa 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.


Das quantidades definidades e da derivada parcial do fluxo de descarga em relação ao tempo temos a respectiva avanço temporal para M:
<math>\dfrac{D^{t + \Delta t}_{i, j} - D^{t}_{i, j}}{\Delta t} + \left [ \dfrac{(Du)^t_ {i+1,j} - (Du)^t_{i-1, j}}{2 \Delta x} \right ] +  \left [ \dfrac{(Dv)^t_ {i,j+1} - (Dv)^t_{i, j-1}}{2 \Delta y} \right ] = 0 \qquad (43)</math>


<math> 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) </math>
<math>\dfrac{(Du)^{t + \Delta t}_{i, j} - (Du)^{t}_{i, j}}{\Delta t} + \left [ \dfrac{(Du^2 + \cfrac{1}{2}gh^2)^t_{i+1, j} - (Du^2 +  \cfrac{1}{2}gh^2)^t_{i-1, j}}{2 \Delta x} \right ] + \left [ \dfrac{(Duv)^t_{i, j+1} - (Duv)^t_{i, j-1}}{2 \Delta y} \right ] = -g D^{t}_{i, j} b_{x. i, j} \qquad (44) </math>


Generalizando a expressão (43) para o fluxo de descarga N temos:
<math>\dfrac{(Dv)^{t + \Delta t}_{i, j} - (Dv)^{t}_{i, j}}{\Delta t} + \left [ \dfrac{(Duv)^t_{i+1, j} - (Duv)^t_{i-1, j}}{2 \Delta x} \right ] + \left [ \dfrac{(Dv^2 + 1/2 gD^2)^t_{i, j+1} - (Dv^2 + 1/2 gD^2)^t_{i, j-1}}{2 \Delta y} \right ] = -g D^{t}_{i, j} b_{y. i, j} \qquad (45) </math>


<math> 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) </math>
Na situação em que temos as equações de águas rasas simplificadas por (39), aplicando FTCS explícito, a evolução temporal da amplitude da onda é:


== Simulações Computacionais de Tsunamis ==
<math>\eta^{t + \Delta t}_{i, j}  = \eta^{t}_{i, j} - \frac{u\Delta t}{2dx} (\eta^t_ {i+1,j} - \eta^t_{i-1, j})  \qquad (46)</math>


=== Solução em 1D ===
Aplicando FTCS implícito temos:


Para a solução em uma direção foi utilizado os mesmos códigos abaixo descritos, desconsiderando a contribuição da direção <math>\vec{y}</math>. Após a descrição dos programas desenvolvidos, será apresentada uma comparação entre ele e o 2D.
<math>\eta^{t + \Delta t}_{i, j} + \frac{u\Delta t}{2dx} (\eta^{t+ \Delta t}_{i+1,j} - \eta^{t- \Delta t}_{i-1, j}) =  \eta^{t}_{i, j}  \qquad (47)</math>


=== Forma conservativa 2D ===
=== Discretização na Forma não Conservativa ===


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.
As equações de águas rasas na forma não conservativa são dadas por (14), (35) e (36), para descrever numericamente estas equações foi utilizada discretização por diferenças finitas, onde realizamos derivadas centradas na região espacial, e para frente na região temporal (FTCS).


<math>\dfrac{h^{t + \Delta t}_{i, j} - h^{t}_{i, j}}{\Delta t} + \left [ \dfrac{(hu)^t_ {i+1,j} - (hu)^t_{i-1, j}}{2 \Delta x} \right ] +  \left [ \dfrac{(hv)^t_ {i,j+1} - (hv)^t_{i, j-1}}{2 \Delta y} \right ] = 0</math>
Discretizando a expressão (14):


<math>\dfrac{hu)^{t + \Delta t}_{i, j} - (hu)^{t}_{i, j}}{\Delta t} + \left [ \dfrac{(hu^2 + \cfrac{1}{2}gh^2)^t_{i+1, j} - (hu^2 +  \cfrac{1}{2}gh^2)^t_{i-1, j}}{2 \Delta x} \right ] + \left [ \dfrac{(huv)^t_{i, j+1} - (huv)^t_{i, j-1}}{2 \Delta y} \right ] = -g h^{t}_{i, j} b_{x. i, j}</math>
<math> \frac{\partial \eta}{\partial t} + \frac{\partial uD}{\partial x}+ \frac{\partial vD}{\partial y} = 0 \qquad (47) </math>


<math>\dfrac{(hv)^{t + \Delta t}_{i, j} - (hv)^{t}_{i, j}}{\Delta t} + \left [ \dfrac{(huv)^t_{i+1, j} - (huv)^t_{i-1, j}}{2 \Delta x} \right ] + \left [ \dfrac{(hv^2 + 1/2 gh^2)^t_{i, j+1} - (hv^2 + 1/2 gh^2)^t_{i, j-1}}{2 \Delta y} \right ] = -g h^{t}_{i, j} b_{y. i, j} </math>
<math> \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 (48) </math>


<math> \eta_{i,j}^{n+1} = \eta_{i,j}^{n} -\Delta t\Bigg(\frac{M_{i+1,j}^{n}-M_{i-1,j}^{n}}{2\Delta x} + \frac{N_{i,j+1}^{n}-N_{i,j-1}^{n}}{2\Delta y}  \Bigg) \qquad (49) </math>


Para os contornos foi utilizado que:
Discretizando a expressão (35):


* <math>u(x_f, y) = - u(x_f - \Delta x, y)</math>
<math> \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 (50) </math>


* <math>u(x_i, y) = - u(x_i + \Delta x, y)</math>
Definindo as quantidades:


* <math>v(x, y_f) = - v(x, y_f - \Delta y)</math>
<math> M2 \equiv \frac{M^2(x,y,t)}{D(x,y,t)}\qquad (51) </math>


* <math>v(x, y_i) = - v(x, y_i + \Delta y)</math>
<math> MN \equiv \frac{M(x,y,t)N(x,y,t)}{D(x,y,t)} \qquad (52) </math>


* Nos contornos de x: <math> \tfrac{\partial h}{\partial x} = 0</math>, discretizando essa derivada temos que: <math>h_{i(+/-)1, j} = h_{i, j}</math>
<math> f(x,y,t) \equiv \frac{gn^{2}}{D^{7/3}} M(M^2 +N^2)^{1/2} \qquad (53) </math>


* Nos contornos de y: <math> \tfrac{\partial h}{\partial y} = 0</math>, discretizando essa derivada temos que: <math>h_{i, j(+/-)1} = h_{i, j}</math>
Das quantidades definidas e da derivada parcial do fluxo de descarga em relação ao tempo temos o respectivo avanço temporal para M:


==== Código ====
<math> 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 (54) </math>
O código foi escrito na linguagem [https://www.python.org/ ''Python''].


<source lang="python">
Generalizando a expressão (54) para o fluxo de descarga N temos:


#%% Bibliotecas
<math> 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 (55) </math>
import numpy as np


import numpy as np
== Simulações Computacionais ==
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
=== Forma conservativa 1D e 2D ===


#%% Parametros
As equações (43),(44) e (45) foram implementadas em python para descrever a evolução temporal das variáveis <math> h</math>, <math> u </math> e <math> v </math> em duas dimensões.


L_xf = 10.0  # m
As condições de contorno dos exemplos obedecem as expressões:
L_x0 = -L_xf


NX = 100
* <math>u(x_f, y) = - u(x_f - \Delta x, y)</math>


dx = (L_xf - L_x0) / NX
* <math>u(x_i, y) = - u(x_i + \Delta x, y)</math>


* <math>v(x, y_f) = - v(x, y_f - \Delta y)</math>


L_yf = 10.0  # m
* <math>v(x, y_i) = - v(x, y_i + \Delta y)</math>
L_y0 = -L_yf


NY = 100
* Nos contornos de x: <math> \tfrac{\partial h}{\partial x} = 0</math>, discretizando essa derivada temos que: <math>h_{i(+/-)1, j} = h_{i, j}</math>


dy = (L_yf - L_y0) / NY
* Nos contornos de y: <math> \tfrac{\partial h}{\partial y} = 0</math>, discretizando essa derivada temos que: <math>h_{i, j(+/-)1} = h_{i, j}</math>


N_INNER = (NX - 2) * (NY - 2)
Para a solução do exemplo em 1D abaixo foi utilizado o mesmo código do conservativo 2D, porém desconsiderando a contribuição da direção <math>\vec{y}</math>.


g = 9.8 # m /s^2
==== Códigos das Equações Conservativas ====
O código 1 foi escrito na linguagem [https://www.python.org/ ''Python''] e resolve os exemplos 1 e 2 que serão apresentados abaixo.


# Tempo
[[Código 1 - Forma Conservativa das Equações de Águas Rasas]]
dt = 0.002
Nt = 1500


time_interval = np.arange(0, Nt*dt, dt)
O código 2 foi escrito na linguagem [https://www.python.org/ ''Python''] e serve para comparar a solução numérica de (39) com a solução exata.


[[Código 2 - Método FTCS explícito]]


#%% Discretização do espaço x-y
O código 3 foi escrito na linguagem [https://www.python.org/ ''Python''] e serve para comparar a solução numérica de (39) com a solução exata.


x = np.linspace(L_x0, L_xf, NX-1) 
[[Código 3 - Método FTCS implícito]]
y = np.linspace(L_y0, L_yf, NY-1)
X, Y = np.meshgrid(x, y)


#%% Condições iniciais, em distribuição gaussiana
==== Resultados ====
 
sigma = 1.0
sigma_v = 1.0


h_2d = 5 * np.ones(shape=np.shape(X)) + np.exp(-(((X)**2 / 2*(sigma**2)) + ((Y)**2 / 2*(sigma**2))))
'''Comparação entre a solução exata e a numérica''': para verificar se o cálculo das equações de águas rasas utilizando FTCS explícito e implícito estão funcionando, partimos da situação simplificada das equações dadas por (46) e (47). Neste problema consideramos a velocidade e a profundidade constantes com  os valores respectivos de 4.3 m/s e 2 m. O deslocamento de água inicial é representado por uma Gaussiana centrada em 80m, de amplitude 2 m e largura 70 m, nesta situação a solução exata de (39) é dada por:
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
<math> \eta = 2 + 3 \cdot exp \Big(-\frac{((x-80) -ut))^2}{70} \Big)  </math>


h_ant = np.copy(h_2d)
Do código 2 obtemos as seguintes curvas ao longo do tempo. Vemos que a solução numérica por FTCS explícito se aproxima do resultado exato para tempos curtos e se afasta da solução exata a medida que o tempo avança.
v_ant = np.copy(v_2d)
u_ant = np.copy(u_2d)


# Inicilização das listas para armazenar os valores
[[Arquivo:Aguas rasas FTCS explicito.gif|thumb|400px|center|Comparação entre a solução exata e a numérica por FTCS explícito. Forma conservativa.]]


hist_h, hist_u, hist_v = [], [], []
Do código 3 obtemos as seguintes curvas ao longo do tempo. Vemos que a solução numérica por FTCS implícito se aproxima do resultado analítico para tempos curtos e se afasta da solução exata a medida que o tempo avança, porém de forma mais lenta que o método explícito.


</source>
[[Arquivo:Aguas rasas FTCS implicito.gif|thumb|400px|center|Comparação entre a solução exata e a numérica por FTCS implícito. Forma conservativa.]]


<math> u \Delta t </math> deve ser significativamente menor que <math> \Delta x </math> para que os métodos FTCS funcionem. 


Função para resolução das equações diferenciais com FTCS.
'''Exemplo 1. Onda confinada em uma caixa com profundidade constante''': obtemos a evolução da amplitude da onda apresentada no GIF abaixo, através do código 1, utilizando que a profundidade é constante, e que tanto a velocidade quanto altura da onda são inicialmente funções gaussianas centradas no espaço.


<source lang="python">
[[Arquivo:conservativa-gausiana.gif|thumb|400px|center|Evolução da amplitude da onda em uma caixa. Forma conservativa]]
# Fator de multiplicacao
fator_x = (dt / (2*(dx)))
fator_y = (dt / (2*(dy)))


def resolve_pdes(h, vx, vy):
Podemos observar que a onda diminui sua amplitude a medida que ela se propaga no espaço. Posteriormente, são observados fenômenos de interferência construtiva e destrutiva devido a reflexão da onda com as paredes da caixa. 


    """
'''Exemplo 2. Onda confinada em uma caixa com profundidade variável''': utilizando a função <math>s = tanh(x)</math> foi possível simular uma onda chegando em uma praia, esse exemplo foi feito em 1D conservativo.
    Função que retorna os valores de profunidade, velocidade em x e em y
    no próximo tempo
   
    Parametros
    -----------
    h : float
                  profundidade no tempo t
    vx : float
        velocidade em x no tempo t
     
    vy : float
        velocidade em y no tempo t
   
    Retorna
    -----------
    prox_h : float
                  profundidade no tempo t + dt
    prox_u : float
        velocidade em x no tempo t + dt
     
    prox_v : float
        velocidade em y no tempo t + dt
   
    """


    # Inicializa os vetores para armazenarem os resultados calculados
Para fazer esse exemplo usamos o código 1, porém desconsideramos as derivadas em <math>x</math> e inicializamos os vetores da altura e velocidade com a discretização em apenas uma dimensão (abaixo descrevemos está parte do código). Deve-se notar que o problema passa a ter um índice, pois a discretização não forma mais uma malha, então pode ser retirado um laço ''for'' do código.
    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]
                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]
<source lang='python'>
                u_r = -vx[i, j]
#%% Discretização do espaço x
                v_r = vy[i, j]
x = np.linspace(L_x0, L_xf, NX)
           
 
            else:
#%% Condições iniciais da superficie e da agua
                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]
           
            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 = - vy[i, j]
           
            else:
               
                h_up = h[i, j + 1]
                v_up = vy[i, j + 1]
                u_up = vx[i, j + 1]
            # --------------------------------------


# Superficie
b = 5 * np.tanh((x- L_xf*0.3) / 2)  # funcao tangente
b += np.full(shape=(np.shape(b)), fill_value=-np.min(b))  # aqui é somado só para a funcao comecar em zero


            ## 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)
                    )


            # # ## Terceira Equação
# Onda e velocidade
       
sigma = 1.0  # distribuicao da onda
            hv_ij = (h[i, j] * vy[i, j]) - \
sigma_v = 1.0 # distribuicao da onda
                    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)
h_2d = 1.1 * max(b) * np.ones(shape=np.shape(x)) + np.exp(-((x)**2/(2*(sigma**2))))  
       
h_2d -= b
    return prox_h, prox_u, prox_v
u_2d =  0.08 * np.exp(-((x)**2/(2*(sigma_v**2))))
</source>
</source>


<source lang="python">
#%% Cálculo


# Resolve para cada passo de tempo
[[Arquivo:gif_conservativa_tanh.gif|thumb|400px|center|Evolução da amplitude da onda em uma caixa com profundidade variável. Forma conservativa]]


for t in time_interval:
=== Forma dissipativa 2D ===
   
 
    h, u, v = resolve_pdes(h_ant, u_ant, v_ant)  # valores das variaveis no tempo = t + dt
   
    # faz isso pq tava dando um erro estranho nesses pontos?
    h[0, :] = h[1, :]
    h[:, 0] = h[:, 1]
   
    v[0, :] = v[1, :]
    v[:, 0] = v[:, 1]
   
    u[0, :] = u[1, :]
    u[:, 0] = u[:, 1]
   
    # 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)
</source>
 
<source lang="python">
#%% Gráfico animado


# Reorganiza os vetores para plotar
Os exemplos que seguem utilizam o código abaixo para calcular as equações de águas rasas não conservativa (38),(43) e (44) nos passos de tempo de <math> \eta </math>, <math> M </math>, <math> N </math>, 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 <ref> KOEHN, Daniel. 2D Shallow Water Equations. Disponível em: <https://github.com/daniel-koehn/Differential-equations-earth-system/blob/master/10_Shallow_Water_Equation_2D/01_2D_Shallow_Water_Equations.ipynb> </ref> foi usado como base.


x_2d = X[0]
==== Códigos das Equações Não Conservativas ====
y_2d = Y[0]


for i in range(1,len(X)):
O código 4 foi escrito na linguagem [https://www.python.org/ ''Python''] e resolve o exemplo 3 que será apresentado abaixo.
   
    x_2d = np.append(x_2d, X[i])
    y_2d = np.append(y_2d, Y[i])


def animate(i):
[[Código 4 - Forma Não Conservativa das Equações de Águas Rasas]]
    plt.clf()  # limpa a figura, pra nao ficar sobrepondo figs
    # ax.cla()
    # titulos
    plt.suptitle('Evolução da onda', fontsize=14)
    plt.title(f'Tempo: {round(dt*i, 3)}', fontsize=12)
   
   
    plt.ylabel('y', fontsize=8)
    plt.xlabel('x', fontsize=8)
   
    # plot
    graph = plt.scatter(x_2d, y_2d, c= hist_h[0::8][i], marker='.')
   
    plt.colorbar()
    # plt.colorbar()
   
    # ax.scatter(x_2d, y_2d, hist_h[0::8][i], marker='.')


    # axis
O código 5 foi escrito na linguagem [https://www.python.org/ ''Python''] e resolve o exemplo 4 que será apresentado abaixo.
    # ax.set_zlim(4.0, 9.0)


    return graph
[[Código 5 - Forma Não Conservativa das Equações de Águas Rasas]]
   


# fig = plt.figure(figsize=(16,9))
==== Exemplo 3 - Onda Confinada em uma Caixa ====
# ax = fig.gca(projection='3d')


fig, ax = plt.subplots()
Através do código 4 fizemos a simulação da propagação de uma onda em uma caixa de <math> L_x =100 </math> m e <math> L_y = 100 </math> m em uma profundidade <math>h(x,y)=50</math> m. O fluxo de descarga e o deslocamento do volume de água inicial, ambas com suas devidas unidades físicas, são Gaussianas centradas em 80, com amplitude igual a 1 e largura igual a 70.
ani = animation.FuncAnimation(fig, animate, frames = len(hist_h[0::8]), repeat=False, interval=0.1)


#%% Salvar o gif
Exemplo 3 com coeficiente de Manning igual a 0.025.
ani.save('onda.gif', writer='imagemagick', fps=5)
</source>


==== Resultados ====
[[Arquivo:Aguas rasas FTCS n=0.025.gif|thumb|400px|center|Evolução da amplitude da onda em uma caixa com profundidade de 50 metros, coeficiente de Manning n = 0.025. Forma dissipativa.]]


1. Utilizando que a superfície é constante e que tanto a velocidade quanto altura da onda são inicialmente funções gaussianas centradas no espaço, obtemos a evolução no espaço que pode ser vista no GIF abaixo.
Exemplo 3 com coeficiente de Manning igual a 20.0.


[[Arquivo:conservativa-gausiana.gif|center|400px|caption:Evolução da altura da onda em uma caixa com superfície constante]]
[[Arquivo:Aguas rasas não conservativo n=20.gif|thumb|400px|center|Evolução da amplitude da onda em uma caixa com profundidade de 50 metros, coeficiente de Manning n = 20. Forma dissipativa.]]


=== Forma dissipativa 2D ===
Dos gráficos acima, podemos observar que quanto maior for o coeficiente de Manning mais rápida a onda irá atenuar a sua amplitude. Está observação faz sentido, já que as forças de fricção sobre o fluído são maiores na situação em que o coeficiente de Manning tem maior valor.


Os exemplos que seguem utilizam as equações de ondas rasas (38),(43) e (44) para calcular os passos de tempo de <math> \eta </math>, <math> M </math>, <math> N </math>, 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 <ref> KOEHN, Daniel. 2D Shallow Water Equations. Disponível em: <https://github.com/daniel-koehn/Differential-equations-earth-system/blob/master/10_Shallow_Water_Equation_2D/01_2D_Shallow_Water_Equations.ipynb> </ref> foi usado como base.
==== Exemplo 4 - Tsunami ====


Simulação da propagação de uma onda em uma caixa de <math> L_x =500</math>  <math> m </math> e <math> L_y =500 </math> <math> m </math>, em uma profundidade variável, que respeita a função tangente hiperbólica. O fluxo de descarga e o deslocamento do volume de água inicial, ambas com suas devidas unidades físicas, são Gaussianas centradas em 300, com amplitude igual a 2 e largura igual a 500.


<source lang="python">
[[Arquivo:Tsunami.gif|thumb|400px|center|Evolução da amplitude da onda em uma caixa com profundidade variáel, coeficiente de Manning n = 0.025. Forma dissipativa.]]
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
</source>


<source lang="python">
A medida que o fundo do oceano diminui sua profundidade a velocidade da onda diminui e a amplitude aumenta. Está simulação mostra o comportamento físico da propagação de um Tsunami em direção ao litoral.
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 
</source>


<source lang="python">
== Comparação entre Métodos ==
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   
</source>
 
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.
 
<source lang="python">
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)
=== Onda confinada em uma caixa ===
        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
Abaixo é mostrado o gráfico de evolução temporal da altura da onda em três pontos distintos do sistema (utilizando os mesmos parâmetros que foi aplicado ao método conservativo com a onda em uma caixa de profundidade constante em todos os outros métodos). Pode-se perceber que com o passar do tempo o movimento das duas equações começa a divergir, mesmo com o fator de fricção baixo.




        fig = plt.figure(figsize=(8.,6.))
[[Arquivo:unificacao_comparacoes.png|thumb|center|700px|Comparação da evolução temporal da altura da onda nos métodos conservativo e dissipativo]]


        fundo = plt.imshow(-h, 'Purples', interpolation = 'nearest', extent = limites)
Para o método 1D é de se esperar que exista menos interferências entre as ondas, pois só existe parede nos finais do sistema, além disso, a onda só precisa se espalhar em uma direção, então a sua amplitude é maior.
        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()
[[Arquivo:unificacao_comparacoes_1D.png|thumb|center|700px|Comparação da evolução temporal da altura da onda nos métodos conservativo e dissipativo 2D e conservativo 1D]]
               
</source>


==== Exemplo 1 - Tsunami Confinada em uma Caixa ====
Também é possível perceber a diferença entre métodos observando apenas em uma direção (<math>\vec{x}</math>), por causa dos termos de fricção e viscosidade.


==== Exemplo 2.1 - Tsunami Propagando-se em Direção a Praia ====
[[Arquivo:unificacao_comparacoes_1D_mesmo.png|thumb|center|700px|Comparação da evolução temporal da altura da onda nos métodos conservativo e dissipativo 1D]]
 
==== Exemplo 2.2 - Tsunami Propagando-se em Direção a Praia ====
 
== Comparação entre Métodos ==


== Referências ==  
== Referências ==  
<references/>
<references/>

Edição atual tal como às 22h03min de 21 de outubro de 2021

Grupo: Gabriel Schmökel e Julia Remus

O objetivo deste trabalho é buscar a solução das equações de águas rasas, por meio de métodos de integração numérica, para resolução de equações diferenciais parciais (EDP's) e posteriormente apresentar uma breve interpretação física das soluções. Demonstramos, nesta página, a derivação das equações, junto com a explicação de cada quantidade física presente. A discretização das equações de águas rasas conservativas e não conservativas são feitas por FTCS explícito, também foi aplicado FTCS implícito para as equações em uma forma simplificada, representada pela equação da deriva. Os resultados obtidos para equação simplificada são comparados com a solução exata, e exemplos mais complexos são solucionados para as formas conservativa e não conservativa. Ao final, uma comparação é feita entre os resultados das equações conservativas e dissipativas.

Introdução

As equações de águas rasas têm aplicações físicas na previsão de tsunâmis, em fluxos atmosféricos, ondas de tempestade e fluxos planetários. Na descrição física dos problemas de fluxos de fluído em ondas, as equações de águas rasas em uma dimensão são dadas por:

As componentes da equação de águas rasas podem ser melhor interpretadas através da seguinte figura:

Componentes das equações de águas rasas

corresponde a amplitude da onda, determina a profundidade do mar em repouso, é o deslocamento total da água, é a velocidade do fluído. Resolvendo a EDP da equação de águas rasas, obtemos como a amplitude da onda se comporta ao longo do tempo e do espaço.

Teoria

Derivação das Equações de Águas Rasas

Iremos demonstrar como chegamos nas equações de águas rasas em duas dimensões, nas formas conservativa e dissipativa, em representações do fluxo de descarga e de velocidade. Posteriormente, tendo as equações em 2D, iremos simplificar elas para a forma unidimensional. Neste processo de demonstração, iremos explicar a interpretação física de cada quantidade presente nas equações.

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:

A equação da continuidade em (3) descreve o balanço de massas para os elementos de volume infinitesimais que pertencem ao fluído, onde a quantidade do lado esquerdo da equação informa o fluxo de massa que entra e sai pelo elemento de volume, e a quantidade do lado direito está relacionada com a massa que se acumula ao longo do tempo [1]. Nesta expressão é a densidade, e é o campo de velocidades, onde u,v e w são as velocidades das partículas que compõe o fluído nas direções x,y,z.

As equações de Navier-Stokes em (4) são balanços diferenciais da quantidade de movimento, obtidas através da aplicação da segunda lei de Newton em cada ponto do escoamento [2] [3] [4].

  • é a aceleração da partícula fluída ao longo do campo de velocidade .
  • está associado as tensões tangenciais e normais atuando sobre os elementos de volume ( é o tensor tensão, as componentes deste tensor são as tensões normais e tangenciais de cisalhamento, expressas por , no qual indica a direção e o plano normal).
  • está associado as pressões que atuam sobre os elementos do fluído.
  • é o vetor aceleração da gravidade atuando sobre os elementos infinitesimais de volume do fluído.

Introduzindo as condições de contorno para a superfície e para a profundidade do oceano [5]:

, onde

, onde

é o deslocamento vertical da água sobre a superfície em repouso, é o vetor velocidade do fluído nas direções horizontais x e y.

A equação da continuidade em (3) pode ser simplificada pelo fato do fluído ser incompressível, isto implica que a densidade é constante.

Integrando a expressão da continuidade em (7), utilizando a regra da integral de Leibniz [6], com os limites indo de até chegamos na seguinte expressão:

Teorema de Leibniz:

Substituindo as condições de contorno da profundidade (6) em (8) obtemos:

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

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

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

Conhecendo as taxas dos fluxos de descarga em relação as regiões espaciais, podemos determinar a taxa da variação da amplitude da onda em relação ao tempo.

Vamos buscar obter as outras duas equações de águas rasas restantes, a partir das quantidades de movimento de Navier-Stokes. Nas componentes x,y e z temos:

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 (15) }

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 (16) }

Na componente z em (17) não consideramos a aceleração das partículas, pois a aceleração da gravidade é muito maior. Também tomamos como nulo as componentes e em (15) e (16), assim passamos a definir . Neste momento estamos desconsiderando as forças de fricção, por isso o tensor tensão também é nulo.

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

Substituindo a pressão em (15):

Integrando a expressão (19), utilizando a regra da integral de Leibniz e as condições de contorno (5) e (6), com os limites indo de até chegamos em outra das equações de águas rasas:

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

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

Com as equações de águas rasas (20), (21), (22) e (23) podemos calcular as taxas de variação dos fluxos de descarga em relação ao tempo. As equações (11), (20) e (21), na representação de velocidades, e as equações (14), (22) e (23), na representação do fluxo de descargas, são as equações de águas rasas conservativas.

Iremos buscar pelas equações de águas rasas não conservativas considerando o tensor de estresse em (4). 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 onde , e as perpendiculares por elementos onde .

Decompondo nas componentes x,y, e z de presente em (4) temos:

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.

Substituindo (27),(28) em (24), integrando em relação a componente z, utilizando a regra de Leibnz e as condições de contorno (5) e (6), obtemos:

Onde é a constante de viscosidade turbulenta, é uma força de resistividade relativa ao movimento do fluído com o fundo do oceano na direção x. Podemos desconsiderar a constante de turbulência na situação em que não temos inclinações abrutas no fundo do mar [7].

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} é } é:

Onde é o coeficiente de fricção, porém o coeficiente de rugosidade de Manning é mais usado, alguns valores deste coeficiente são [8]:

Material Coeficiente de Rugosidade de Manning
Cimento puro e metal liso 0,010
Terra lisa 0,017
Pedras, ervas daninhas 0,035
Péssimo relevo de canal 0,060
Bom relevo de canal 0,025

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

Substituindo (32) em (31) obtemos:

Generalizando a expressão (29) para a componente y.

Adicionando, repectivamente, (33) e (34) nas expressões (22) e (24), obtemos as equações de águas rasas considerando as forças de fricção do fundo do oceano.

As equações (14), (35) e (36), na representação do fluxo de descargas, são as equações de águas rasas não conservativas.

Em uma dimensão podemos expressar as equações de águas rasas eliminando a componente das expressões (14),(35) e tomando o fluxo de descarga como nulo.

Simplificação das Equações de Águas Rasas

As equações de águas rasas podem ser simplificadas para equação de advecção através das seguintes considerações:

  • A velocidade do fluído é constante.
  • A profundidade do fundo do oceano é constante.

De (39) e das considerações acima temos:

Discretização na 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 [9][10]. Para desenvolvê-lo são necessárias algumas premissas:

  • O comprimento da onda é muito maior que as contribuições na direção
  • A aceleração na direção da velocidade na direção é zero
  • As componentes das velocidades em e em ( e ) não variam em


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

Onde é a altura do fluido desde a base, são as velocidades médias na direções , é a constante gravitacional e é função que descreve a superfície onde acontece o movimento.

Para descrever numericamente as equações de águas rasas na forma conservativa 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.

Na situação em que temos as equações de águas rasas simplificadas por (39), aplicando FTCS explícito, a evolução temporal da amplitude da onda é:

Aplicando FTCS implícito temos:

Discretização na Forma não Conservativa

As equações de águas rasas na forma não conservativa são dadas por (14), (35) e (36), para descrever numericamente estas equações foi utilizada discretização por diferenças finitas, onde realizamos derivadas centradas na região espacial, e para frente na região temporal (FTCS).

Discretizando a expressão (14):

Discretizando a expressão (35):

Definindo as quantidades:

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

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

Simulações Computacionais

Forma conservativa 1D e 2D

As equações (43),(44) e (45) foram implementadas em python para descrever a evolução temporal das variáveis , e em duas dimensões.

As condições de contorno dos exemplos obedecem as expressões:

  • Nos contornos de x: , discretizando essa derivada temos que:
  • Nos contornos de y: , discretizando essa derivada temos que:

Para a solução do exemplo em 1D abaixo foi utilizado o mesmo código do conservativo 2D, porém desconsiderando a contribuição da direção .

Códigos das Equações Conservativas

O código 1 foi escrito na linguagem Python e resolve os exemplos 1 e 2 que serão apresentados abaixo.

Código 1 - Forma Conservativa das Equações de Águas Rasas

O código 2 foi escrito na linguagem Python e serve para comparar a solução numérica de (39) com a solução exata.

Código 2 - Método FTCS explícito

O código 3 foi escrito na linguagem Python e serve para comparar a solução numérica de (39) com a solução exata.

Código 3 - Método FTCS implícito

Resultados

Comparação entre a solução exata e a numérica: para verificar se o cálculo das equações de águas rasas utilizando FTCS explícito e implícito estão funcionando, partimos da situação simplificada das equações dadas por (46) e (47). Neste problema consideramos a velocidade e a profundidade constantes com os valores respectivos de 4.3 m/s e 2 m. O deslocamento de água inicial é representado por uma Gaussiana centrada em 80m, de amplitude 2 m e largura 70 m, nesta situação a solução exata de (39) é dada por:

Do código 2 obtemos as seguintes curvas ao longo do tempo. Vemos que a solução numérica por FTCS explícito se aproxima do resultado exato para tempos curtos e se afasta da solução exata a medida que o tempo avança.

Comparação entre a solução exata e a numérica por FTCS explícito. Forma conservativa.

Do código 3 obtemos as seguintes curvas ao longo do tempo. Vemos que a solução numérica por FTCS implícito se aproxima do resultado analítico para tempos curtos e se afasta da solução exata a medida que o tempo avança, porém de forma mais lenta que o método explícito.

Comparação entre a solução exata e a numérica por FTCS implícito. Forma conservativa.

deve ser significativamente menor que para que os métodos FTCS funcionem.

Exemplo 1. Onda confinada em uma caixa com profundidade constante: obtemos a evolução da amplitude da onda apresentada no GIF abaixo, através do código 1, utilizando que a profundidade é constante, e que tanto a velocidade quanto altura da onda são inicialmente funções gaussianas centradas no espaço.

Evolução da amplitude da onda em uma caixa. Forma conservativa

Podemos observar que a onda diminui sua amplitude a medida que ela se propaga no espaço. Posteriormente, são observados fenômenos de interferência construtiva e destrutiva devido a reflexão da onda com as paredes da caixa.

Exemplo 2. Onda confinada em uma caixa com profundidade variável: utilizando a função foi possível simular uma onda chegando em uma praia, esse exemplo foi feito em 1D conservativo.

Para fazer esse exemplo usamos o código 1, porém desconsideramos as derivadas em e inicializamos os vetores da altura e velocidade com a discretização em apenas uma dimensão (abaixo descrevemos está parte do código). Deve-se notar que o problema passa a ter um índice, pois a discretização não forma mais uma malha, então pode ser retirado um laço for do código.

#%% Discretização do espaço x
x = np.linspace(L_x0, L_xf, NX)
  
#%% Condições iniciais da superficie e da agua

# Superficie
b = 5 * np.tanh((x- L_xf*0.3) / 2)  # funcao tangente
b += np.full(shape=(np.shape(b)), fill_value=-np.min(b))  # aqui é somado só para a funcao comecar em zero


# Onda e velocidade
sigma = 1.0   # distribuicao da onda
sigma_v = 1.0  # distribuicao da onda

h_2d = 1.1 * max(b) * np.ones(shape=np.shape(x)) + np.exp(-((x)**2/(2*(sigma**2)))) 
h_2d -= b
 
u_2d =  0.08 * np.exp(-((x)**2/(2*(sigma_v**2))))


Evolução da amplitude da onda em uma caixa com profundidade variável. Forma conservativa

Forma dissipativa 2D

Os exemplos que seguem utilizam o código abaixo para calcular as equações de águas rasas não conservativa (38),(43) e (44) nos passos de tempo de , , , 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 [11] foi usado como base.

Códigos das Equações Não Conservativas

O código 4 foi escrito na linguagem Python e resolve o exemplo 3 que será apresentado abaixo.

Código 4 - Forma Não Conservativa das Equações de Águas Rasas

O código 5 foi escrito na linguagem Python e resolve o exemplo 4 que será apresentado abaixo.

Código 5 - Forma Não Conservativa das Equações de Águas Rasas

Exemplo 3 - Onda Confinada em uma Caixa

Através do código 4 fizemos a simulação da propagação de uma onda em uma caixa de m e m em uma profundidade m. O fluxo de descarga e o deslocamento do volume de água inicial, ambas com suas devidas unidades físicas, são Gaussianas centradas em 80, com amplitude igual a 1 e largura igual a 70.

Exemplo 3 com coeficiente de Manning igual a 0.025.

Evolução da amplitude da onda em uma caixa com profundidade de 50 metros, coeficiente de Manning n = 0.025. Forma dissipativa.

Exemplo 3 com coeficiente de Manning igual a 20.0.

Evolução da amplitude da onda em uma caixa com profundidade de 50 metros, coeficiente de Manning n = 20. Forma dissipativa.

Dos gráficos acima, podemos observar que quanto maior for o coeficiente de Manning mais rápida a onda irá atenuar a sua amplitude. Está observação faz sentido, já que as forças de fricção sobre o fluído são maiores na situação em que o coeficiente de Manning tem maior valor.

Exemplo 4 - Tsunami

Simulação da propagação de uma onda em uma caixa de e , em uma profundidade variável, que respeita a função tangente hiperbólica. O fluxo de descarga e o deslocamento do volume de água inicial, ambas com suas devidas unidades físicas, são Gaussianas centradas em 300, com amplitude igual a 2 e largura igual a 500.

Evolução da amplitude da onda em uma caixa com profundidade variáel, coeficiente de Manning n = 0.025. Forma dissipativa.

A medida que o fundo do oceano diminui sua profundidade a velocidade da onda diminui e a amplitude aumenta. Está simulação mostra o comportamento físico da propagação de um Tsunami em direção ao litoral.

Comparação entre Métodos

Onda confinada em uma caixa

Abaixo é mostrado o gráfico de evolução temporal da altura da onda em três pontos distintos do sistema (utilizando os mesmos parâmetros que foi aplicado ao método conservativo com a onda em uma caixa de profundidade constante em todos os outros métodos). Pode-se perceber que com o passar do tempo o movimento das duas equações começa a divergir, mesmo com o fator de fricção baixo.


Comparação da evolução temporal da altura da onda nos métodos conservativo e dissipativo

Para o método 1D é de se esperar que exista menos interferências entre as ondas, pois só existe parede nos finais do sistema, além disso, a onda só precisa se espalhar em uma direção, então a sua amplitude é maior.

Comparação da evolução temporal da altura da onda nos métodos conservativo e dissipativo 2D e conservativo 1D

Também é possível perceber a diferença entre métodos observando apenas em uma direção (), por causa dos termos de fricção e viscosidade.

Comparação da evolução temporal da altura da onda nos métodos conservativo e dissipativo 1D

Referências

  1. Equação da continuidade mássica: balanços de massa diferenciais. Bloom Consultoria.Disponível em: <https://www.youtube.com/watch?v=pEip-GvO0LM&list=PL1yqHjPQz-Lqjri07DqZ3RsSWJfICvdiu&index=3>
  2. Equação de Navier-Stokes (Parte 1) - Derivadas materiais. Bloom Consultoria.Disponível em: https://www.youtube.com/watch?v=FLoZODPpayM
  3. Equação de Navier-Stokes (Parte 2) - Equação diferencial da quantidade de movimento. Bloom Consultoria.Disponível em: https://www.youtube.com/watch?v=e06ZRzdO4iM
  4. Equação de Navier-Stokes (Parte 3) - Tensões Normais e Cisalhantes. Bloom Consultoria.Disponível em: https://www.youtube.com/watch?v=na2kGOSYNv8
  5. SEGUR, Harvey; YAMAMOTO, Hiroki. Lecture 8: The Shallow-Water Equations.Disponível em: <https://docplayer.net/49487265-Lecture-8-the-shallow-water-equations.html>
  6. Leibniz integral rule. Disponível em: https://en.wikipedia.org/wiki/Leibniz_integral_rule
  7. 7,0 7,1 IMAMURA, Fumihiko.Tsunami Modelling Manual.Disponível em: http://www.tsunami.civil.tohoku.ac.jp/hokusai3/J/projects/manual-ver-3.1.pdf
  8. LINSLEY, Ray K.; FRANZINI, Joseph B. Water Resources Engineering.
  9. GARCÍA-NAVARRO, P; et al. The shallow water equations: An example of hyperbolic system. Espanha: 2008. Disponível em: <https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.571.1364&rep=rep1&type=pdf>
  10. KUHBACHER, Christian. Shallow Water: Derivation and Applications. Disponível em: <http://www.mathematik.tu-dortmund.de/lsiii/cms/papers/Kuehbacher2009.pdf>
  11. KOEHN, Daniel. 2D Shallow Water Equations. Disponível em: <https://github.com/daniel-koehn/Differential-equations-earth-system/blob/master/10_Shallow_Water_Equation_2D/01_2D_Shallow_Water_Equations.ipynb>