Equação de Águas Rasas

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

Introdução

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

Formação de um Tsunami

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

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

IMAGEM

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

IMAGEM

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

IMAGEM

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

IMAGEM

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

Teoria

Derivação das EQs. de Águas Rasas

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

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

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

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

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

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

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

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

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

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

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

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

Teorema de Leibniz:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Substituindo a pressão em (14):

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Substituindo (30) em (29) obtemos:

Falhou ao verificar gramática (função desconhecida '\math'): {\displaystyle \frac{\tau_x}{\rho} = \frac{gn^{2}}{D^{7/3}} M(M^2 +N^2)^{1/2} \qquad (30) <\math> === 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: * O comprimento da onda é muito maior que as contribuições na direção <math>\vec{z}}

• A aceleração na direção da velocidade na direção ${\vec {z}}$ é zero
• As componentes das velocidades em ${\vec {x}}$ e em ${\vec {y}}$ (${\vec {u}}$ e ${\vec {v}}$) não variam em ${\vec {z}}$

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

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

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

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

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

Desenvolvimento do cálculo

Forma conservativa 2D

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

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

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

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

No desenvolvimento do programa não foi conseguido alcançar os resultados esperados, pois o sistema não converge após um tempo. Mesmo assim, vamos apresentar o código ? escrito na linguagem Python.

#%% Bibliotecas 
import numpy as np

import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation
import matplotlib.patches as mpatches

#%% Parametros

L_xf = 10.0  # m
L_x0 = -L_xf

NX = 100

dx = (L_xf - L_x0) / NX


L_yf = 10.0  # m
L_y0 = -L_yf

NY = 100

dy = (L_yf - L_y0) / NY

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

g = 9.8 # m /s^2

# Tempo
dt = 0.002
Nt = 1000

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

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

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

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

sigma = 0.6
sigma_v = 0.6

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