Um Método Integral de Contorno para a Modelagem da Propagação de Ondas Internas em um Sistema de Dois Fluidos

RESUMO

Este trabalho tem como objetivo utilizar a formulação integral de contorno na construção de um método numérico para modelar a propagação de ondas internas na interface entre dois fluidos. Apresentamos vários exemplos numéricos para ilustrar a acurácia do método proposto e também mostrar sua utilidade na simulação das interações de ondas não lineares.

Palavras-chave:
ondas aquáticas internas; método da integral de contorno; discretização do operador Dirichlet-Neumann

ABSTRACT

In this work, we use boundary integral formulation to develop a numerical method to study the propagation of internal waves at the interface of two fluids. We present several numerical examples in order to illustrate the accuracy of the proposed method and also to show that its usefulness in the simulation of nonlinear waves interactions.

Keywords:
internal water waves; boundary integral method; discretization of the Dirichlet-Neumann operator.

1 INTRODUÇÃO

As diferenças de temperatura e salinidade da água provocam a formação de camadas de diferentes densidades em águas oceânicas, lagos e reservatórios. Sob os efeitos da força gravitacional, as partículas localizadas na interfaces entre duas camadas são deslocadas produzindo movimentos ondulatórios. De fato, quando uma parcela do fluido acima da interface é deslocada para baixo a força de empuxo produz o deslocamento para cima de uma parte do fluido abaixo da interface, consequentemente se este fluido for mais pesado do que o seu entorno então a força restauradora da gravidade irá desloca-lo para baixo novamente movimentando também a própria interface. Esse movimento de cada parcela do fluido, análogo ao sistema massa-mola, acontece ao longo de toda a interface que como consequência oscila no espaço e no tempo.

De forma mais geral, chamamos de ondas internas ao movimento das interfaces entre as camadas de fluidos com diferentes densidades¹⁰. Nas águas oceânicas, o período desse tipo de onda varia de alguns minutos até meia hora e o comprimento de onda pode variar entre algumas centenas de metros e vários quilômetros.

Um caso muito especial de ondas internas, são as ondas que se formam na interface entre a água e o ar. Nesta situação o efeito da camada de ar pode ser desprezado e os movimentos ondulatórios da superfície da água são chamados de ondas aquáticas de superfície.

No estudo da propagação de ondas internas os efeitos da viscosidade e da fricção podem ser desprezados pois o comprimento dessas ondas é relativamente grande^{10), (11}, assim os modelos mais simples consideram o caso quando as camadas são formadas por fluidos não viscosos, incompressíveis e homogêneos que não se misturam, e estão separadas por interfaces bem definidas. O escoamento em cada camada é governado pelas equações de Euler, nas interfaces a pressão e a velocidade normal são contínuas e nas fronteiras rígidas a velocidade normal é conhecida^{10), (12), (14}.

As simplificações descritas acima podem ser aplicadas sob certas restrições: a estratificação em camadas homogêneas é aplicável quando as variações espaciais da densidade acontecem em uma escala espacial bem menor que o comprimento das ondas em estudo; a incompressibilidade representa uma forma simples de filtrar (homogeneizar) os efeitos das ondas sonoras; a continuidade da velocidade normal e seu comportamento nas fronteiras rígidas é indicativo de que os efeitos da fricção e da rugosidade são desprezíveis; a continuidade da pressão é válida quando desconsideramos os efeitos da tensão superficial. Mais detalhes sobre as diferentes hipóteses consideradas e outros modelos de ondas internas podem ser encontrados em^{10), (11}.

Destacamos que modelos matemáticos baseados nas equações de Euler também são muito usados no estudo das ondas aquáticas de superfície (veja^{12), (14} e as referências contidas nele).

Mesmo considerando apenas duas camadas, esse tipo de modelo é muito útil, por exemplo pode ser usado na descrição de regiões costeiras onde a água proveniente de um rio se encontra com a água do mar, e também em regiões oceânicas onde a diferença de salinidade ocorre de forma abrupta¹⁰. Na aplicação desse modelo, o efeito da interface água-ar é desprezado já que as ondas de superfície correspondentes apresentam amplitudes bem mais baixas que as ondas que aparecem na interface entre as duas camadas de água com diferentes densidades, considera-se assim que a camada superior da água está limitada por uma tampa rígida.

Além disso, o estudo desse modelo é bastante desafiador do ponto de vista matemático e computacional. Como consequência, modelos simplificados ou reduzidos que exploram diferentes regimes assintóticos são amplamente estudados na literatura e devido a sua simplicidade também são usados frequentemente em simulações numéricas¹⁴.

As equações de Euler tem sido usadas com sucesso no desenvolvimento de métodos numéricos para a simulação de ondas aquáticas de superfície (veja por exemplo^{5), (17} e as referências contidas lá). Em geral, nesses trabalhos é usado o operador Dirichlet-Neumann para reduzir a quantidade de incógnitas e a dimensão do espaço subjacente^{14), (19}. Dessa forma o problema original é transformado de um problema onde temos incógnitas em cada ponto da região ocupada pelo fluido, em um onde as incógnitas estão associadas apenas com a fronteira dessa região. A aproximação eficiente e precisa desse operador tem representado um dos maiores desafios nesses trabalhos. No trabalho recente¹⁶ é apresentada uma comparação de vários métodos desenvolvidos com esse objetivo.

Os métodos baseados na formulação do operador Dirichlet-Neumann através da integral de contorno possuem boa precisão e estabilidade numérica¹⁶. Uma aplicação recente desse método foi apresentada em⁹, para fazer um estudo numérico das órbitas de partículas associadas com ondas de Stokes progressivas sobre a influência de uma correnteza uniforme. Na aproximação do operador Dirichlet-Neumann os autores usam várias ideias apresentadas em⁶. Destacamos que no estudo desse problema não foi necessário modelar a evolução temporal da onda, o que limita a aplicabilidade do método numérico desenvolvido.

Para o problema de ondas internas, em¹³ foi apresentado um método numérico baseado em uma formulação integral das equações de Euler para estudar o comportamento de ondas solitárias (levando em conta o efeito da tensão superficial). Na formulação do problema não foi usado o operador Dirichlet-Neumann e também não foi necessário introduzir explicitamente a evolução no tempo. O cálculo de ondas solitárias (desconsiderando o efeito da tensão superficial) foi tratado em⁷ aplicando uma formulação integral. Nesse trabalho também foi apresentado um método numérico baseado na formulação integral para o caso de ondas internas em 3 dimensões, mas usando apenas os primeiros termos da expansão do operador Dirichlet-Neumann.

Assim, o principal objetivo deste trabalho é apresentar um método numérico para o cálculo da evolução temporal de ondas internas em um sistema de dois fluidos usando o operador Dirichlet-Neumann sem precisar introduzir qualquer tipo de truncamento e também mostrar exemplos numéricos da sua aplicação. O método proposto é baseado na representação desse operador através de uma integral de contorno usando resultados apresentados em⁶.

Destacamos que as simplificações introduzidas no modelo de ondas internas tem como objetivo produzir um método numérico simples cuja implementação computacional não demande o uso de técnicas mais elaboradas, mas ainda possa ser usado para estudar problemas interessantes. Por isso, neste trabalho é considerado o caso de duas dimensões espaciais. Também devemos enfatizar que devido à forte não linearidade, esse modelo apesar de simples ainda não tem sido compreendido totalmente do ponto de vista teórico³.

O restante deste trabalho está organizado da seguinte forma, na seção 2 apresentamos as equações de Euler que modelam a propagação de ondas internas na interface entre dois fluidos e fazemos a redução do número de variáveis e da dimensão espacial. Na seção 3 é apresentado o método integral de contorno. Na seção 4 são apresentados e discutidos os resultados numéricos das simulações. As conclusões foram apresentadas na seção 5.

2 MODELAGEM DE ONDAS INTERNAS EM UM SISTEMA DE DOIS FLUIDOS

2.1 Equações de Euler

Nesta seção, vamos apresentar um modelo simplificado para a propagação de ondas internas a partir do qual será desenvolvido o método numérico. Neste modelo estão sendo consideradas todas as simplificações mencionadas na introdução.

Consideraremos a evolução temporal da interface entre duas camadas de fluidos de diferentes densidades e que não se misturam, o fluido mais denso ocupa a camada inferior. Os fluidos são inviscidos, incompressíveis, homogêneos e estão sob a influência de um campo gravitacional constante. Consideraremos que os fluidos estão limitados por uma tampa e um fundo planos sem rugosidades, e que o escoamento é irrotacional e bidimensional. Adotaremos o índice 1 para referenciar os elementos da camada superior e 2 para os elementos da camada inferior conforme ilustra a Figura 1. A posição da interface está definida pela equação y = η(x, t), onde x, y e t representam as coordenadas horizontal e vertical, e o tempo, respectivamente. A posição de repouso da interface corresponde à reta y = 0, as posições da tampa e do fundo são definidas pelas retas y = h ₁ e y = -h ₂, respectivamente, em que h ₁, h ₂ > 0. A densidade dos fluidos são representadas por ρ ₁ e ρ ₂.

Figura 1
Configuração física do sistema com duas camadas de fluidos.

Vamos considerar que as camadas estão definidas pelas equações

$\begin{array}{ll}\hfill & {\Omega }_1=\left\{0\le x\le 2L\text{e}\eta \left(x,t\right)\le y\le h_1\right\}\hfill \\ \hfill & \partial {\Omega }_{1r}=\left\{y=h_1\right\}\hfill \\ \hfill & {\Omega }_2=\left\{0\le x\le 2L\text{e}-h_2\le y\le \eta \left(x,t\right)\right\}\hfill \\ \hfill & \partial {\Omega }_{2r}=\left\{y=-h_2\right\}\hfill \end{array}$

onde 2L > 0 representa o comprimento do domínio de interesse e ∂Ω_ir corresponde às fronteiras horizontais desse domínio.

Assim as equações de Euler para o sistema de dois fluidos¹² tem a forma

$\{\begin{array}{ll}\hfill & {\varphi }_{i,xx}+{\varphi }_{i,yy}=0em{\Omega }_i\hfill \\ \hfill & {\varphi }_{i,y}=0em\partial {\Omega }_{ir}\hfill \\ \hfill & {\eta }_t+{\varphi }_{i,x}{\eta }_x-{\varphi }_{i,y}=0emy=\eta \left(x,t\right)\hfill \\ \hfill & {\varphi }_{i,t}+\frac{1}{2}\left({\varphi }_{i,x}^2+{\varphi }_{i,y}^2\right)+\frac{p_i}{{\rho }_i}+g\eta =0emy=\eta \left(x,t\right)\hfill \\ \hfill & p_1=p_{2}emy=\eta \left(x,t\right)\hfill \end{array}$ (2.1)

em que _^ϕi , _^ρi e _^ρi (i = 1, 2) representam o potencial de velocidade, a densidade e a pressão correspondente à camada i, respectivamente, e os índices t, x, y estão associados com derivadas parciais. Por simplicidade nas fronteiras laterais consideramos condições de contorno periódicas, ou seja todas as funções envolvidas são periódicas na variável x com período 2L.

A condição de continuidade da pressão ao atravessar a interface, do ponto de vista físico significa que a tensão superficial é desprezível, a introdução dos seus efeitos não apresenta nenhuma dificuldade.

As simplificações associadas com a tampa e o fundo plano, o escoamento bidimensional e as condições de contorno periódicas tem como objetivo desenvolver um método numérico usando aproximações de alta precisão e também obter um modelo discretizado com poucos graus de liberdade.

2.2 Formulação integral de contorno e os operadores Neumann-Dirichlet e Dirichlet-Neumann

Nesta seção introduzimos os operadores de Neumann-Dirichlet e Dirichlet-Neumann a partir da formulação integral de contorno².

Para uma função harmônica na região Ω de ℝ², a formulação integral de contorno permite relacionar os valores desta função na fronteira ∂Ω com os valores da sua derivada normal na fronteira². Dessa forma podemos definir os operadores Dirichlet-Neumann e Neumann-Dirichlet, que relacionam os dados de Dirichlet e Neumann correspondentes a uma função harmônica¹⁴.

Seja Ω uma região de ℝ² com fronteira ∂Ω suave por partes e u uma função harmônica em Ω. A partir da terceira Identidade de Green⁸, temos a formulação integral de contorno a seguir²

$\pi u\left(P_0\right)={\sout{\int }}_{\partial \Omega }\left(G\left(P;P_0\right)\frac{\partial u}{\partial n_P}\left(P\right)-u\left(P\right)\frac{\partial G}{\partial n_P}\left(P;P_0\right)\right)d{s}_P,P_0\in \partial \Omega$ (2.2)

em que P ∊ ∂Ω, _^nP representa o vetor normal à fronteira externo a Ω em P e G(P; P ₀) é uma função harmônica em relação a P com singularidade logarítmica - 1n| P - P ₀| em P = P ₀. Observe que, devido à singularidade de G, a integral (2.2) é definida no sentido do Valor Principal de Cauchy. A partir de (2.2) definem-se os seguintes operadores integrais que atuam sobre as funções w definidas na fronteira ∂Ω

$Aw\left(P_0\right):={\sout{\int }}_{\partial \Omega }w\left(P\right)G\left(P;P_0\right)d{s}_P,P_0\in \partial \Omega$ (2.3)

$\mathcal{B}w\left(P_0\right):=\underset{\partial \Omega }{{\displaystyle \int }}w\left(P\right)\frac{\partial G}{\partial n_P}\left(P;P_0\right)d{s}_p,P_0\in \partial \Omega$ (2.4)

Levando em conta as relações (2.2)-(2.4), os operadores Dirichlet-Neumann e Neumann-Dirichlet denotados como 𝒟 e 𝒩, respectivamente, são definidos por

$\mathcal{D}w\left(P_0\right):={\mathcal{A}}^{-1}\left[\pi \mathcal{I}+\mathcal{B}\right]w\left(P_0\right),P_0\in \partial \Omega$ (2.5)

$\mathcal{N}w\left(P_0\right):={[\pi \mathcal{I}+\mathcal{B}]}^{-1}\mathcal{A}w\left(P_0\right),P_0\in \partial \Omega$ (2.6)

em que ℐ representa o operador identidade. Esses operadores são responsáveis por relacionar as condições de fronteira de Dirichlet e Neumann correspondentes a uma função harmônica. Observamos que a função 𝒩w está definida a menos de uma constante, portanto para defini-la de maneira única é necessário impor alguma condição adicional, por exemplo, que tenha média zero.

2.3 Sistema de equações reduzidas

Sistemas reduzidos equivalentes ao sistema de equações de Euler podem ser encontrados em^{14), (19}, mas a forma que apresentamos aqui tem como objetivo simplificar o desenvolvimento dos métodos numéricos.

Para reduzir o número de incógnitas, vamos considerar os traços dos potenciais da velocidade na interface e reescrever as equações de Euler incluindo apenas esses traços e a própria interface. Nesse processo os operadores de Neumann-Dirichlet e Dirichlet-Neumann desempenham um papel muito importante.

Começamos desenvolvendo uma descrição dos operadores 𝒟 e 𝒩 mais apropriada e especifica para o sistema de duas camadas. Para isso, as nossas considerações serão feitas com respeito a camada inferior (i = 2) e as relações referentes a camada superior serão apresentadas posteriormente com as devidas adaptações. Com isso, afim de simplificar a notação, será desconsiderado o índice i. Como ϕ é uma função harmônica na região Ω, os operadores (2.5) e (2.6), nos dão relações entre ϕ e _^ϕn para P ₀ ∊ ∂Ω, onde _^ϕn representa a derivada normal. Sendo assim, temos

${\varphi }_n\left(x,\eta \left(x,t\right),t\right)=\mathcal{D}\varphi \left(x,\eta \left(x,t\right),t\right),\varphi \left(x,\eta \left(x,t\right),t\right)=\mathcal{N}{\varphi }_n\left(x,\eta \left(x,t\right),t\right).$

Utilizando o Método das Imagens para um ponto P ₀ = (x ₀, y ₀), consideraremos a função harmônica G(P; P ₀) = K(x - x ₀, y - y ₀) + K(x - x ₀, -2h - y - y ₀) em que K é a função de Green dada por ⁶

$K\left(x,y\right)=-\frac{1}{2}\ln \left[1+e^{2\pi \frac{y}{L}}-2e^{\pi \frac{y}{L}}\cos \left(\frac{\pi x}{L}\right)\right].$

Tal escolha para G é feita pois, além de ser periódica e possuir apenas singularidades logarítmicas, também satisfaz a condição de Neumann do problema em ∂Ω_r . Considerando ∂Ω = Υ_d ⋃ Γ ⋃ Υ_e ⋃ ∂Ω_r , em que Γ representa a interface e Υ_d , Υ_e são as partes laterais da fronteira ∂Ω. Note que (2.3) e (2.4) restritos à Υ_d e Υ_e se anulam devido à periodicidade das funções envolvidas e em ∂Ω_r é nulo devido ao Método das Imagens. Sendo assim, os operadores 𝒜 e ℬ estão restritos apenas a Γ como pode ser visto a seguir,

$A{\varphi }_n\left(P_0\right)=\underset{0}{\overset{2L}{\sout{\int }}}\left[{\varphi }_n\left(P_s\right)G\left(P_s;P_0\right)\right)]\parallel \left(\mathrm{1,}{\eta }_s\left(s\right)\right)\parallel ds$

$\mathcal{B}\varphi \left(P_0\right)={\int }_0^{2L}\left[\varphi \left(P_s\right)\frac{\partial G}{\partial n}\left(P_s;P_0\right)\right]\parallel \left(\mathrm{1,}{\eta }_s\left(t,s\right)\right)\parallel ds$

em que _{^{Ps =}} (x(s), y(s)) representa uma parametrização de Γ.

Considerando o traço do potencial na interface Φ(x, t) = ϕ(x, η(x, t), t), suas derivadas normais, tangenciais e as derivadas com respeito à x e y em y = η (x, t) estão relacionadas de modo que:

$\widetilde{{\Phi }_n}\left(x,t\right)=\parallel \left(\mathrm{1,}{\eta }_x\left(x,t\right)\right)\parallel {\Phi }_n\left(x,t\right)=\left[-{\eta }_x{\varphi }_x+{\varphi }_y\right]\left(x,\eta \left(x,t\right),t\right)$ (2.7)

${\Phi }_x\left(x,t\right)=\parallel \left(\mathrm{1,}{\eta }_x\left(x,t\right)\right)\parallel {\Phi }_{\tau }\left(x,t\right)=\left[{\varphi }_x+{\varphi }_y{\eta }_x\right]\left(x,\eta \left(x,t\right),t\right)$ (2.8)

${\Phi }_t\left(x,t\right)=\left[{\varphi }_t+{\varphi }_y{\eta }_t\right]\left(x,\eta \left(x,t\right),t\right)$

onde τ representa o vetor tangente à Γ e Φ_τ a derivada na direção de τ.

A partir de (2.7) e (2.8) encontramos

${\varphi }_x\left(x,\eta \left(x,t\right),t\right)=\frac{{\Phi }_x-{\tilde{\Phi }}_n{\eta }_x}{1+{\eta }_x^2\left(x,t\right)}\text{e}{\varphi }_y\left(x,\eta \left(x,t\right),t\right)=\frac{{\eta }_x{\Phi }_x+{\tilde{\Phi }}_n}{1+{\eta }_x^2\left(x,t\right)}.$

Retornando com a notação que indica a camada abordada, η₁ e η₂ representam, respectivamente, os vetores normais externos à ∂Ω₁ e ∂Ω₂ e na interface Γ das duas camadas adotamos η = η₂ = ?η₁.

Substituindo o resultado anterior e (2.7) na terceira e quarta equação de (2.1), temos que

$\{\begin{array}{l}{\eta }_t={\tilde{\Phi }}_{\mathrm{2,}n}\hfill \\ {\Phi }_{\mathrm{2,}t}=-\frac{1}{2}\left[\frac{{\Phi }_{\mathrm{2,}x}^2-{\tilde{\Phi }}_{\mathrm{2,}n}^2}{1+{\eta }_x^2}\right]+\frac{\eta {\Phi }_{\mathrm{2,}x}{\tilde{\Phi }}_{\mathrm{2,}n}}{1+{\eta }_x^2}-g\eta -\frac{p_2}{{\rho }_2}.\hfill \end{array}$

Note que, ${\Phi }_{\mathrm{1,}{n}_1}=-{\Phi }_{\mathrm{2,}{n}_2}$, logo Φ_{1, n} .= Φ_{2, n} . Portanto, os termos de (2.2) que dependem do vetor normal, quando referentes a camada superior, sofrerão alterações. Feitas as devidas adaptações, os operadores envolvidos serão tais que

$\begin{array}{ll}{\mathcal{D}}_1\left(\eta \right){\Phi }_1\hfill & ={[-{\mathcal{A}}_1\left(\eta \right)]}^{-1}\left[\pi \mathcal{I}-{\mathcal{B}}_1\left(\eta \right)\right]{\Phi }_1\hfill \\ {\mathcal{N}}_1\left(\eta \right){\Phi }_{\mathrm{1,}n}\hfill & ={[\pi \mathcal{I}-{\mathcal{B}}_1\left(\eta \right)]}^{-1}\left[-{\mathcal{A}}_1\left(\eta \right)\right]{\Phi }_{\mathrm{1,}n}\hfill \\ {\mathcal{D}}_2\left(\eta \right){\Phi }_2\hfill & ={[{\mathcal{A}}_2\left(\eta \right)]}^{-1}\left[\pi \mathcal{I}+{\mathcal{B}}_2\left(\eta \right)\right]{\Phi }_2\hfill \\ {\mathcal{N}}_2\left(\eta \right){\Phi }_{\mathrm{2,}n}\hfill & ={[\pi \mathcal{I}+{\mathcal{B}}_2\left(\eta \right)]}^{-1}\left[{\mathcal{A}}_2\left(\eta \right)\right]{\Phi }_{\mathrm{2,}n}\hfill \end{array}$

com

$\begin{array}{ll}{\mathcal{A}}_i\left(\eta \right){\Phi }_{i,n}\left(P_0\right)\hfill & ={\sout{\int }}_{\Gamma }{\Phi }_{i,n}\left(P_s\right)G_i\left(P_s;P_0\right)ds\hfill \\ {\mathcal{B}}_i\left(\eta \right){\Phi }_i\left(P_0\right)\hfill & ={\int }_{\Gamma }{\Phi }_i\left(P_s\right)\frac{\partial G_i}{\partial n}\left(P_s;P_0\right)ds\hfill \end{array}$

onde _^Gi representa a função de Green referente a camada i com respectiva profundidade. Com esses operadores, é possível relacionar os Φ₁ e Φ₂ pois

$\begin{array}{ll}{\Phi }_1\hfill & ={\mathcal{N}}_1\left(\eta \right){\Phi }_{\mathrm{1,}n}\hfill \\ \hfill & ={\mathcal{N}}_1\left(\eta \right){\Phi }_{\mathrm{2,}n}\hfill \\ \hfill & ={\mathcal{N}}_1\left(\eta \right)\left[{\mathcal{D}}_2\left(\eta \right){\Phi }_2\right].\hfill \end{array}$ (2.9)

Além disso, como p ₁(x, y, t) = p ₂(x, y, t) em y = η(x, t), considerando Ψ = ρ ₁Φ₁ - ρ ₂Φ₂ temos

$\begin{array}{l}{\Psi }_t=-\frac{1}{2}\left[\frac{{\rho }_1{\Phi }^2{}_{\mathrm{1,}x}-{\rho }_2{\Phi }^2{}_{\mathrm{2,}x}-\left({\rho }_1-{\rho }_2\right){\tilde{\Phi }}^2{}_{\mathrm{1,}n}}{1+{\eta }^2{}_x}\right]\text{}\\ +\frac{{\eta }_x{\Psi }_x{\tilde{\Phi }}_{\mathrm{1,}n}}{1+{\eta }^2{}_x}-\left({\rho }_1-{\rho }_2\right)g\eta .\end{array}$ (2.10)

Ao utilizar o operador Neumann-Dirichlet, para que função Φ_i fique definida de maneira única, consideraremos que ${\widehat{\Phi }}_i$ = Φ_i - ??Φ_i ?? em que ??Φ_i ?? representa a média da função Φ_i em Γ. Assim, a partir de (2.9), concluímos que

$\widehat{\Psi }=\left[{\rho }_1{\mathcal{N}}_1\left(\eta \right){\mathcal{D}}_2\left(\eta \right)-{\rho }_2\mathcal{I}\right]{\widehat{\Phi }}_2.$

Definimos os operadores ℛ₁ e ℛ₂ que são responsáveis pela recuperação de ${\widehat{\Phi }}_1$ e ${\widehat{\Phi }}_2$ a partir de $\widehat{\Psi }$ segundo as equações

${\widehat{\Phi }}_2={\mathcal{R}}_2\left(\eta \right)\Psi ={\left[{\rho }_1{\mathcal{N}}_1\left(\eta \right){\mathcal{D}}_2\left(\eta \right)-{\rho }_2\mathcal{I}\right]}^{-1}\widehat{\Psi }$ (2.11)

${\widehat{\Phi }}_1={\mathcal{R}}_1\left(\eta \right)\Psi ={\rho }_1^{-1}\left(\widehat{\Psi }-{\rho }_2{\mathcal{R}}_2\left(\eta \right)\Psi \right).$ (2.12)

Portanto, para recuperar os potenciais a partir de Ψ, consideraremos que Φ₁ possui média zero, portanto

$${\Phi }_1={\mathcal{R}}_1\left(\eta \right)\Psi {\Phi }_2={\mathcal{R}}_2\left(\eta \right)\Psi -{\rho }_2^{-1}〈\Psi 〉.$$

Unindo (2.10) com o resultado anterior, temos o seguinte sistema

$\{\begin{array}{l}{\eta }_t=\mathcal{F}\left(\eta ,\Psi \right)\hfill \\ {\Psi }_t=\mathcal{G}\left(\eta ,\Psi \right)\hfill \end{array}$ (2.13)

onde ℱ e 𝒢 são da seguinte forma

$$\{\begin{array}{l}\mathcal{F}\left(\eta ,\Psi \right)={\mathfrak{D}}_1{\Re }_1\Psi ,\hfill \\ \begin{array}{l}\mathcal{G}\left(\eta ,\Psi \right)=-\frac{1}{2}\left[\frac{{\rho }_1{\left({\Re }_1\Psi \right)}_x^2-{\rho }_2{\left({\Re }_2\Psi \right)}_x^2-\left({\rho }_1-{\rho }_2\right){\left({\mathfrak{D}}_1{\Re }_1\Psi \right)}^2}{1+{\eta }_x^2}\right]\\ \begin{array}{ccc}& & \end{array}+\frac{{\eta }_x{\Psi }_x\left({\mathfrak{D}}_1{\Re }_1\Psi \right)}{1+{\eta }_x^2}-\left({\rho }_1-{\rho }_2\right)g\eta .\end{array}\hfill \end{array}$$

Note que para simplificar a notação na definição de ℱ e 𝒢 foi omitida a dependência de η nos diferentes operadores envolvidos.

O sistema (2.13) deve ser tratado como um problema de valor inicial, para isso incluímos condições iniciais apropriadas η(x, 0) e Ψ(x, 0) (no tempo t = 0). Observamos que se a posição da interface η(x, 0) e pelo menos um dos potenciais Φ_i (x, 0) forem conhecidos, então será possível determinar uma condição inicial Ψ(x, 0), compatível com esse problema.x

Na próxima seção serão apresentados métodos numéricos para a solução aproximada deste sistema.

3 MÉTODOS NUMÉRICOS E ALGORITMO

3.1 Tratamento das integrais singulares

Uma característica intrínseca dos operadores A e 𝔅 definidos na seção anterior é a presença de termos singulares, pelo fato de K(x - x ₀, y - y ₀) ter uma singularidade em P = P ₀.

A fim de estudar o comportamento de A, para cada tempo t fixo, introduzimos a desingularização da função G proposta em¹⁸ (veja também²)

$$\begin{array}{l}G\left(x,y;x_0,y_0\right)=K\left(x-x_0,y-y_0\right)-\frac{1}{2}\ln \left\{{\left[2\sin \left(\frac{\pi \left(x-x_0\right)}{2L}\right)\right]}^2\right\}\text{}\\ \begin{array}{cccc}& & & \end{array}+\frac{1}{2}\ln \left\{{\left[2\sin \left(\frac{\pi \left(x-x_0\right)}{2L}\right)\right]}^2\right\}+K\left(x-x_0,2h-y-y_0\right).\end{array}$$

Assim, rescrevemos G como G (x, y; x ₀, y ₀) = K ₁(x, y; x ₀, y ₀) + K ₂(x - x ₀) onde

$$\begin{array}{l}{K}_1\left(x,y;x_0,y_0\right)=\frac{\pi }{2L}\left(y-y_0\right)-\frac{1}{2}\ln \left\{1+\frac{\sin h^2\left(\frac{\pi \left(y-y_0\right)}{2L}\right)}{{\sin }^2\left(\frac{\pi \left(x-x_0\right)}{2L}\right)}\right\}\text{}\\ \begin{array}{ccccc}& & & & \end{array}+K\left(x-x_0,2h-y-y_0\right),\text{}\\ \begin{array}{cc}& \end{array}{K}_2\left(x-x_0\right)=\frac{1}{2}\ln \left[4{\sin }^2\left(\frac{\pi \left(x-x_0\right)}{2L}\right)\right].\end{array}$$

Por simplicidade, omitiremos a dependência temporal de η. Observe que quando x → x ₀, sinh²(η(x) - η(x ₀)) e sin²(x - x ₀)) são proporcionais a (η(x) - η(x ₀))² e (x - x ₀)², respectivamente. Usando a expansão de Taylor obtemos

$\underset{x\to x_0}{\lim }{K}_1\left(x-x_0,\eta \left(x\right)-\eta \left(x_0\right)\right)=-\ln \sqrt{1+{\eta }_x^2\left(x_0\right)}+K\left(\mathrm{0,2}h-2\eta \left(x_0\right)\right).$

Finalmente, obtemos a decomposição 𝒜(η) = 𝒜₀ + $\tilde{\mathcal{A}}$(η), em que

$\begin{array}{l}\tilde{\mathcal{A}}\left(\eta \right){\varphi }_n\left(x_0,y_0\right)={\sout{\int }}_{\Gamma }\left[{\varphi }_n\left(s,\eta \left(s\right)\right)K_1\left(s,\eta \left(s\right);x_0,y_0\right))\right]ds\text{}\\ {\mathcal{A}}_0{\varphi }_n\left(x_0,y_0\right)={\sout{\int }}_{\Gamma }\left[{\varphi }_n\left(s,\eta \left(s\right)\right)K_2\left(s-x_0\right)\right]ds.\end{array}$

Para o operador ℬ, quando $\frac{\partial G}{\partial n_P}$ é expandido em Série de Taylor na vizinhança de x ₀, obtemos que

$\begin{array}{l}\underset{x\to x_0}{\lim }\frac{\partial G}{\partial n}\left(P_x,P_0\right)\parallel \left(\mathrm{1,}{\eta }_x\left(x\right)\right)\parallel \\ =\frac{{\eta }_{xx}\left(x_0\right)}{2+2{\eta }_x^2\left(x_0\right)}+\frac{\partial K}{\partial y}\left(\mathrm{0,}2\left(h-\eta \left(x_0\right)\right)\right)-\frac{\pi }{2L}.\end{array}$

Portanto, nesse caso o efeito da singularidade de K será eliminado se redefinimos apropriadamente $\frac{\partial G}{\partial n}\left(P_x,P_0\right)|\left|\left(\mathrm{1,}{\eta }_x\left(x\right)\right)\right|$ para x = x ₀.

3.2 Discretização e esquema numérico

Para a discretização das equações utilizamos vários métodos numéricos disponíveis na literatura^{2), (15}. A escolha dos métodos foi feita com o intuito de se obter um esquema numérico com boa acurácia para isto foram selecionados métodos que apresentam precisão espectral.

Consideramos uma discretização de Γ associada à discretização do intervalo [0, 2L] dada por _{^{xj = j}} Δx para j = 0, 1, 2, ..., _^Nx - 1 com $\Delta x=\frac{2L}{N_x}$ em que _^Nx é um inteiro par. Dessa forma temos que determinar as funções vetoriais $\eta \left(t\right),\Psi \left(t\right)\in {ℝ}^{N_x}$ cujas componentes aproximam η(_{^{xj, t}} ) e Ψ(_{^{xj, t}} ), respectivamente.

Os operadores $\tilde{\mathcal{A}}$ e ℬ foram discretizados aproximando as integrais envolvidas pela Regra dos Trapézios. Já para 𝒜₀ foi usada sua representação através de uma série de Fourier e a mesma foi aproximada utilizando a Transformada Discreta de Fourier¹⁵. Dessa forma, os operadores 𝒜 e ℬ foram discretizados na forma matricial e a partir deles foi possível aproximar os operadores Dirichlet-Neumann, Neumann-Dirichlet assim como os outros operadores envolvidos.

Para aproximar as derivadas de primeira e segunda ordem de η e Φ_i com relação à variável x foram utilizadas as matrizes de diferenciação espectral $D_N^{\left(1\right)}$ e $D_N^{\left(2\right)}$ obtidas a partir da Transformada Discreta de Fourier, conforme apresentado em¹⁵.

Finalmente, todas essas aproximações foram usadas para aproximar a parte direita do sistema (2.13), obtendo dessa forma um sistema de equações diferenciais ordinárias (EDOs) para determinar η (t), Ψ(t). É importante destacar que devido à periodicidade das funções envolvidas todas as aproximações introduzidas possuem uma alta ordem de precisão O((Δx)^m ) em que m ≥ 1 representa a ordem de diferenciabilidade da solução. Em particular, quando a solução é uma função analítica temos precisão espectral $O\left(e^{-\frac{c}{\Delta x}}\right)$ onde c > 0 é uma constante que depende da solução^{2), (15}.

Para a discretização do tempo utilizamos o esquema Leap-frog¹⁵. Assim, a partir do sistema (2.13), escolhendo um passo de tempo Δt as aproximações de η (nΔt), Ψ(nΔt) representadas por _^ηn , Ψ ⁿ , respectivamente, são calculadas pelo seguinte esquema totalmente discreto

$\{\begin{array}{l}{\eta }^{n+1}={\eta }^{n-1}+2\Delta tF\left({\eta }^n,{\Psi }^n\right)\hfill \\ {\Psi }^{n+1}={\Psi }^{n-1}+2\Delta tG\left({\eta }^n,{\Psi }^n\right)\hfill \end{array}paran\ge \mathrm{1,}$ (3.1)

onde F(_^ηn , Ψ ⁿ ) e G(_^ηn , Ψ ⁿ ) são obtidos a partir dos métodos aproximados discutidos acima.

Este esquema possui ordem de precisão quadrática e tem um baixo custo computacional no sentido que é preciso fazer uma única avaliação das funções F e G em cada passo de tempo.

Como o esquema Leap-frog é um esquema de dois passos, é necessário conhecer os valores nos passos n - 1 e n para obter uma aproximação para o passo n + 1. Porém, por se tratar de um problema de valor inicial, conhecemos apenas η ⁰ e Ψ ⁰ . Uma maneira de aproximar a solução para n = 1 é utilizando um esquema explícito de passo um, por exemplo, o esquema de Euler explícito ¹⁵. Em conclusão, no esquema numérico proposto a partir das aproximações iniciais η ⁰ e Ψ ⁰:

• Calculamos η ¹ e Ψ ¹ usando o esquema de Euler explícito

${\eta }^1={\eta }^0+\Delta tF\left({\eta }^0,{\Psi }^0\right),{\Psi }^1={\Psi }^0+\Delta tG\left({\eta }^0,{\Psi }^0\right).$

• Para n ≥ 2 calculamos _^ηn e Ψ ⁿ aplicando o esquema Leap-frog (3.1).

É importante notar que de acordo com os métodos aproximados escolhidos para a discretização das derivadas espaciais, dos operadores integrais envolvidos e do sistema de EDOs correspondente, temos que a ordem de precisão do esquema numérico proposto é O((Δx)^m ) + O((Δt)²) em que m ≥ 1 depende da ordem de diferenciabilidade da solução^{2), (15}. Além disso, para soluções analíticas em x temos precisão espectral em relação a essa variável. Em relação à estabilidade numérica, para o sistema linearizado com respeito ao estado de repouso temos uma condição de Courant-Friedrichs-Lewy (CFL) da forma _^NC = _^co Δt/Δx ≤ 1 onde _^co é a velocidade de propagação da onda e _^NC é o número de Courant¹⁵. Espera-se que para ondas de pequena amplitude seja suficiente restringir mais um pouco a condição CFL, usando por exemplo a condição _^NC ≤ C ₁ ≤ 1, para garantir a estabilidade e nesse caso obter uma ordem de convergência compatível com a ordem de precisão. Destacamos que uma análise rigorosa da estabilidade numérica e a convergência do esquema proposto será desenvolvida em trabalhos futuros.

Na próxima seção, apresentaremos alguns resultados numéricos obtidos a partir de uma imple mentação desse algoritmo realizada no Matlab.

4 RESULTADOS NUMÉRICOS

Para validar o algoritmo e sua implementação, vamos comparar os resultados numéricos com soluções aproximadas disponíveis na literatura. Posteriormente faremos simulações onde a so lu ção é desconhecida.

4.1 Estudo de convergência

Considerando os parâmetros da discretização Δt e Δx determinamos o erro na aproximação de η como

$E\left(\eta \right)=\max \left\{\left|\eta \left(j\Delta x,n\Delta t\right)-{\eta }_j^n\right|:j=\mathrm{0,}\dots ,N_x-\mathrm{1,}n=\mathrm{0,}\dots \right\}$

e de forma semelhante define-se o erro E(Ψ). Considerando que os parâmetros da discretização Δt e Δx satisfazem a condição Δt = O(Δx) espera-se que os erros sejam de ordem quadrática em Δt, ou seja E = O(Δt ²).

Para diferentes discretizações com tamanho de passo Δ_^ti obtemos os erros _^Ei (i = 0, 1, ...), assim a taxa de convergência empírica _^pi , observada nessas simulações será calculada de acordo com a equação

$$p_i=\frac{\log \left(E_i/E_{i-1}\right)}{\log \left(\Delta t_i/\Delta t_{i-1}\right)},i=\mathrm{1,}\dots .$$

De acordo com as observações apresentadas acima espera-se obter valores de _^pi ≈ 2, ou seja _^Ei ≈ C(Δ_^ti )² onde C > 0 é uma constante que não depende da discretização.

Nos exemplos a seguir adotamos g = 9.81 m/s², as densidades dos fluidos ρ₁ = 1000 kg/m³ e ρ₂ = 1023 kg/m³, as profundidades das camadas h ₁ = 10 m, h ₂ = 200 m e um comprimento de onda λ = 30 m.

Como a velocidade de propagação da onda _^co é conhecida (aproximadamente) a discretização foi feita de forma tal que o número de Courant $N_C=\frac{c_o\Delta t}{\Delta x}=\frac{2}{3}$.

4.1.1 Ondas estacionárias de pequena amplitude

A solução do sistema linearizado (aproximação de primeira ordem) para o regime de ondas estacionárias próximas do repouso apresentado em¹² é dado por:

$\begin{array}{l}\eta \left(x,t\right)=a\sin \left(\omega t\right)\cos \left(kx\right)\text{}\\ {\Phi }_1\left(x,t\right)=-\frac{\omega }{k}a\coth \left(k{h}_1\right)\cos \left(\omega t\right)\cos \left(kx\right)\text{}\\ {\Phi }_2\left(x,t\right)=-\frac{\omega }{k}a\coth \left(k{h}_2\right)\cos \left(\omega t\right)\cos \left(kx\right)\end{array}$ (4.1)

em que k = 2π/λ é o número de onda, ω a frequência e a representa a amplitude. Substituindo (4.1) em (2.1) temos a relação de dispersão a seguir.

${\omega }^2=\frac{gk\left({\rho }_2-{\rho }_1\right)}{{\rho }_1\coth \left(k{h}_1\right)+{\rho }_2\coth \left(k{h}_2\right)}$ (4.2)

Os erros obtidos nas simulações e as taxas de convergência empírica, considerando diferentes valores de a, são apresentados nas Tabelas 1-3. Os valores dos erros para η e Ψ apresentados nas tabelas foram adimensionalizados de acordo com

$\overline{E}\left(\eta \right)=\frac{E\left(\eta \right)}{a},\overline{E}\left(\Psi \right)=\frac{E\left(\Psi \right)}{{\rho }_{1}a\omega /k}$

Tabela 1
Erro e taxa de convergência empírica de η e Ψ para ondas estacionárias com a = 10^{- 1}.

		Erro (adimensional)		Taxa de convergência
Δx	Δt	η	Ψ	η	Ψ
0.9375	0.8453	0.018	27.3772	-	-
0.4688	0.4270	5.3321 × 10-3	6.7831	2.0603	2.0135
0.2344	0.2146	9.5037 × 10-4	1.5271	2.1885	2.1505
0.1172	0.1076	2.1425 × 10-4	0.2541	2.1492	2.5892
0.0586	0.0539	1.5195 × 10-4	0.2803	0.4957	-0.1424

Tabela 2
Erro e taxa de convergência empírica de η e Ψ para ondas estacionárias com a = 10^{- 2}.

		Erro (adimensional)		Taxa de convergência
Δx	Δt	η	Ψ	η	Ψ
0.9375	0.8453	1.8238 × 10-2	27.6423	-	-
0.4688	0.4270	4.5025 × 10-3	7.0003	2.0181	1.9803
0.2344	0.2146	1.1209 × 10-3	1.7535	2.0061	1.9996
0.1172	0.1076	2.7872 × 10-4	0.4418	2.0078	1.9924
0.0586	0.0539	6.8375 × 10-5	0.1131	2.0273	1.9606

Tabela 3
Erro e taxa de convergência empírica de η e Ψ para ondas estacionárias com a = 10^{- 3}.

		Erro (adimensional)		Taxa de convergência
Δx	Δt	η	Ψ	η	Ψ
0.9375	0.8453	1.8239 × 10-2	27.6147	-	-
0.4688	0.4270	4.5042 × 10-3	7.0417	2.0177	1.9809
0.2344	0.2146	1.1226 × 10-3	1.7950	2.0044	2.0012
0.1172	0.1076	2.8042 × 10-4	0.4384	2.0012	1.9965
0.0586	0.0539	7.0079 × 10-5	0. 1099	2.0005	1.9955

Observamos que quanto menor o valor de a menores são os erros, esse resultado é esperado já que estamos calculando os erros em relação a uma solução aproximada. Também temos que quando a aproximação linear é suficientemente boa (a = 10^-3) o método apresenta taxas de convergência empírica _^pi próximas de 2 que é a ordem de convergência do método Leap-frog usado na discretização no tempo.

4.1.2 Ondas de Stokes de 3ª Ordem

No seguinte exemplo vamos um pouco além da teoria linear, considerando ondas progressivas de Stokes de terceira ordem¹⁰. As solução para este tipo de ondas é dada por

$\eta \left(x,t\right)=a_1\cos \left(kx-\omega t\right)+a_2\cos 2\left(kx-\omega t\right)+a_3\cos 3\left(kx-\omega t\right)$ (4.3)

em que os coeficientes de a ₂, a ₃ e ω dependem de a ₁, k, h ₁ e h ₂ (as expressões podem ser achadas em¹⁰, p. 33])

Os resultados das simulações apresentados nas Tabelas 4 e 5 mostram que houve uma melhoria nos resultados numéricos se comparado com os resultados do caso linearizado. A taxa de convergência empírica também ficou próxima do valor esperado p = 2.

Tabela 4
Erro e taxa de convergência empírica de η para ondas de Stokes de terceira ordem com a₁ = 10^{- 1}.

		Erro (adimensionalizado)	Precisão
Δx	Δt	η	η
1.8750	1.7083	0.0775	-
0.9375	0.8453	0.0196	2.0249
0.4688	0.4270	5.3031 × 10-3	1.8435
0.2344	0.2146	1.9393 × 10-3	1.4513
0.1172	0.1076	1.1018 × 10-3	0.8157

Tabela 5
Erro e taxa de convergência empírica de η para ondas de Stokes de terceira ordem com a₁ = 10^{- 3}.

		Erro	Precisão
Δx	Δt	η(x,t)	η(x,t)
1.8750	1.7083	6.2499 × 10-2	-
0.9375	0.8453	1.6252 × 10-2	2.0016
0.4688	0.4270	4.2349 × 10-3	1.9693
0.2344	0.2146	1.0951 × 10-3	1.9659
0.1172	0.1076	2.8332 × 10-4	1.9579
0.0586	0.0539	7.6833 × 10-4	1.8862

4.2 Simulação de ondas não lineares

4.2.1 Evolução de um perfil harmônico

Neste exemplo consideramos a evolução de uma perturbação inicial correspondente a um perfil harmônica de amplitude grande.

Consideramos g=9.81 m/s², as densidades dos fluidos ρ ₁ = 1000 kg/m³ e ρ ₂ = 1023 kg/m³ kg/m³ e as profundidades das camadas h ₁ = 17.2 m, h ₂ = 57.8 m.

Como condição inicial adotamos a perturbação harmônica

$\begin{array}{l}\eta \left(x,0\right)=a\cos \left(kx\right),\text{}\\ {\Phi }_1\left(x,t\right)=\frac{\omega }{k}a\coth \left(k{h}_1\right)\sin \left(kx\right),\text{}\\ {\Phi }_2\left(x,t\right)=\frac{\omega }{k}a\coth \left(k{h}_2\right)\sin \left(kx\right)\end{array}$

em que ω é dado por (4.2). Fazemos a = 7 m e o número de onda $k=\frac{\pi }{200}{m}^{-1}$ (correspondente a λ = 400 m), notamos que a solução de Stokes de 3^a ordem (4.3) considerando a ₁ = 7 m não possui sentido físico pois a altura dessa onda resulta em um valor maior que h ₁.

O domínio físico nas simulações corresponde ao intervalo [0, 2λ] e consideramos Δx = 12.5 m e o Δt corresponde a um número de Courant _^NC = 0.2 (considerando a velocidade de propagação da onda linearizada).

Na Figura 2 mostramos o perfil da interface no momento inicial e após 26118.0 s. Observamos que a altura da onda final é próxima de 10.5 m e que a largura da crista é maior que a do vale, isto mostra o caráter não linear da onda. Os gráficos da Figura 3 mostram a evolução temporal da interface em dois locais diferentes com coordenadas x = 0 e x = λ/2. Eles sugerem que no movimento da interface acontece uma superposição de ondas com períodos curtos e longos.

Figura 2
Evolução de uma perturbação inicial harmônica.

Figura 3
Gráficos da evolução temporal do deslocamento da interface nos locais com coordenadas x = 0 e x = λ/2.

Determinando o espectro de Fourier de η(x, t), na Figura 4 mostramos sua magnitude normalizada. Observamos que o espectro se encontra concentrado apenas em alguns poucos pontos. De fato, temos 14 pontos da forma (k’, ω’) em que os k’ ficam próximos de ±k, ±2k e ±3k com k = π/λ ≈ 0.157 m^-1, e k ₂ ≈ 2k ₁ e k ₃ ≈ 3k ₁ , entanto que ω’ serão próximos de ±ῶ, ±2ῶ e ±3ῶ com ῶ ≈ 0.0269 s^-1. Notamos que ῶ está próximo do valor da frequência de uma onda linear correspondente a k dado por (4.2). Assim a evolução da interface é dada pela superposição de ondas harmônicas em ressonância com o perfil inicial.

Figura 4
Magnitude do espectro de Fourier do sinal η(x, t).

4.2.2 Decomposição em ondas solitárias

No seguinte exemplo, consideramos soluções do tipo ondas solitárias, mais especificamente mostraremos a decomposição de uma perturbação inicial da interface em dois grupos de ondas solitárias ¹.

Sabemos que quando h ₁ ≪ h ₂ = O(λ) existem soluções aproximadas (assintóticas) na forma de ondas solitárias que dependem do parâmetro 0 < β < π/2 e satisfazem a equação de ondas longas intermediárias (ILW)¹⁴. Essa solução é dada por

${\eta }_{ILW}\left(x,t\right)=\frac{a{\cos }^2\beta }{{\cos }^2\beta +{\sinh }^2\left(\frac{x-ct}{{\lambda }_{ILW}}\right)}$ (4.4)

em que a = 4c ₂ β tan β/(c ₁ h ₂), c = c ₀ - 2c ₂ β tan β/h ₂ e λ_ILW = h ₂/β representam a amplitude, velocidade de propagação e comprimento de onda efetivo, respectivamente, observe que seu calculo envolve os seguintes parâmetros: $c_0^2=g{h}_1^2$(ρ ₂ - ρ ₁)/ρ ₁, c ₁ = -1.5c ₀/h ₁ e c ₂ = 0.5c ₀ h ₁ ρ ₂/ρ ₁ (que dependem apenas da configuração física).

Neste exemplo consideramos como condições iniciais uma interface com perfil η ₀(x) = 2^_ηILW (0, x) (para o caso β = 1) e potenciais nulos. Consideramos g = 9.81 m/s², as densidades dos fluidos ρ ₁ = 1000 kg/m³ e ρ ₂ = 1023 kg/m³ e as profundidades das camadas h ₁ = 10 m, h ₂ = 200 m. Observamos que neste caso o perfil inicial representa uma onda de depressão. O domínio físico nas simulações foi o intervalo [0, 2000] usamos um Δx = 15.6250 m e um Δt = 1.4128 s que corresponde a um número de Courant ^_NC = 0.125 (com relação a velocidade de propagação da onda solitária).

Na Figura 5 os resultados mostram a decomposição do perfil inicial em dois pulsos viajando para a direita e a esquerda, respectivamente. Observamos que esses dois pulsos se deslocam sem praticamente mudar sua forma, devido à periodicidade em x elas voltam a interagir e se reencontrar no centro.

Figura 5
Decomposição de um perfil inicial em dois pulsos viajando para a direita e a esquerda, respectivamente. Mostramos os gráficos de -η(x, t).

Na Figura 6 é apresentado em detalhe do pulso que viaja para a direita quando t = 867.0664 em conjunto com o perfil de uma onda solitária correspondente à equação ILW (4.4). Podemos notar que a onda solitária viaja um pouco mais devagar, mas comparando o perfil do pulso com a onda solitária deslocado observamos que eles são muito semelhantes.

Figura 6
Detalhe do perfil da onda viajando para esquerda. Comparação com o perfil de uma onda solitária da equação ILW.

4.2.3 Colisão de ondas solitárias

Neste exemplo apresentamos a colisão entre duas ondas solitárias correspondentes à equação ILW. Consideramos os valores de ρ ₁, ρ ₂, h ₁ e h ₂ são os mesmos que no exemplo anterior. O domínio físico nas simulações foi o intervalo [0, 6000] usamos um Δx = 15 m e um Δt que corresponde a um número de Courant ^_NC = 0.2 (com relação a velocidade da onda solitária mais rápida). A primeira onda é dada por (4.4) com β = 0.8 e a segunda corresponde a β = 1.

A Figura 7 (a) mostra o perfil das duas ondas separadas inicialmente uma distância de 3000 m no tempo inicial e sua evolução após 2170.1 s. Uma imagem das trajetórias das duas ondas é mostrada na Figura 7 (b), onde pode ser observado que após transcorridos aproximadamente 1000 s acontece a colisão das duas ondas solitárias. Notamos que as trajetórias são muito próximas de linhas retas.

Figura 7
Colisão de duas ondas solitárias.

Na Figura 8 (a), mostramos os perfis das duas ondas pouco tempo antes da colisão e no momento da colisão. Observamos que antes da colisão a forma das ondas era praticamente a mesma que no momento inicial e no momento da colisão ambas se unem para formar uma única onda. Na Figura 8 (b), apresentamos os perfis das ondas no momento da colisão e após a mesma. Notamos que após a interação direta entre eles o pulso maior vai deixando para atrás um pequeno rastro (consequência do caráter dispersivo do problema).

Figura 8
Perfil das ondas instantes antes, durante e após a colisão.

5 CONCLUSÕES

Neste trabalho apresentamos um método integral de contorno para a solução numérica do problema da propagação de ondas internas na interface entre duas camadas de fluidos. Nos exemplos numéricos mostramos que o comportamento dos erros e a taxa de convergência do método dependem principalmente da discretização temporal que foi usando o esquema Leap-frog, visto que os métodos associados com a discretização no espaço possuem uma alta ordem de precisão. Além disso, foram mostrados outros exemplos onde os efeitos não lineares jogam um papel importante. Os resultados indicam que o método proposto pode ser usado na modelagem computacional das interações não lineares que acontecem na propagação de ondas internas.

REFERÊNCIAS

• ¹
  M. Ablowitz & H. Segur. Solitons and the Inverse Scattering Transform. SIAM, (1981).
• ²
  K. Atkinson & W. Han. Theoretical Numerical Analysis. Springer Science, (2009).
• ³
  C. Bardos & A. Fursikov (ed.). Instability in Models Connected with Fluid Flows I, II. Springer, (2008).
• ⁴
  L. C. Evans. Partial Differential Equations. 2nd. ed., AMS, (2010).
• ⁵
  D. Fructus, D. Clamond, J. Grue & Ø. Kristiansen. An efficient model for three-dimensional surface wave simulations Part I: Free space problems. J. of Computational Physics, 205 (2005), 665-685.
• ⁶
  P.A. Guidotti. A first kind boundary integral formulation for the Dirichlet-to-Neumann map in 2D. J. of Computational Physics, 190 (2008), 325-345.
• ⁷
  J. Grue. Interfacial Wave Motion of Very Large Amplitude: Formulation in Three Dimensions and Numerical Experiments. Procedia IUTAM, 8 (2013), 129-143.
• ⁸
  F. John. Partial Differential Equations. 4th ed., Springer-Verlag, (1982).
• ⁹
  A. Nachbin & R. Ribeiro-Junior. A boundary integral formulation for particle trajectories in Stokes waves. Discrete and Continuous Dynamical Systems, 34 (2014), 3135-3153.
• ¹⁰
  R. Massel. Internal Gravity Waves in the Shallow Seas. Springer, (2015).
• ¹¹
  Yu.Z. Miropolsky & O.D. Shiskina. Dynamics of Internal Gravity Waves in the Ocean. Springer, (2001).
• ¹²
  H. Lamb. Hydrodynamics. Cambridge Univ. Press, (1895).
• ¹³
  O. Laget & F. Dias. Numerical computation of capillary? gravity interfacial solitary waves. J. Fluid Mech., 349 (1997), 221-251.
• ¹⁴
  J-C. Saut. Asymptotic models for surface and internal waves. Publicações matemáticas, Sociedade Brasileira de Matemática, (2013).
• ¹⁵
  L. N. Trefethen. Spectral Methods in MATLAB. SIAM, (2001).
• ¹⁶
  J. Wilkening & V. Vasan. Comparison of five methods of computing the Dirichlet-Neumann operator for the water wave problem. Contemp. Math., 635 (2015), 175-210.
• ¹⁷
  L. Xu & P. Guyenne. Numerical simulation of three-dimensional nonlinear water waves. J. of Computational Physics, 228 (2009), 8446-8466.
• ¹⁸
  Y. Yan & I.H. Sloan. On integral equations of the first kind with logarithmic kernels. J. of Integral Equations and Applications, 1 (1988), 549-579.
• ¹⁹
  V.E. Zakharov & E. Kutznetzov. Hamiltonian formalism for nonlinear waves. Uspekhi Fiz. Nauk, 11 (1997), 1137-1167.

Datas de Publicação

• Publicação nesta coleção
  May-Aug 2017

Histórico

• Recebido
  30 Set 2015
• Aceito
  21 Abr 2017