Grupo2 - Ondas1

Introdução

Integrantes do grupo: Rodrigo Zamin Ferreira (262692), Leonardo Xavier Rodrigues (262696), Maurício Gomes de Queiroz () e Rodrigo Lopes de Sousa Silva (262705)

A modelagem numérica vem se tornando cada vez mais uma ferramenta indispensável para um engenheiro. Tal modelagem pode trazer informações importantes para entender como melhor abordar o desenvolvimento de um projeto, neste caso, um que envolva ondas. Nós, como futuros engenheiros físicos, pensamos em trazer um problema mais "concreto", de engenharia costeira e portuária, que pode ou não surgir em nossas vidas profissionais mas cujo método de solução certamente estará presente. Aqui será apresentado um modelo baseado em uma condição inicial e um perfil topográfico do local estudado que descreve a evolução temporal de uma onda.

Para testarmos os diferentes métodos, utilizaremos a equação da onda em uma dimensão, que é uma equação diferencial parcial de segunda ordem, para modelarmos uma corda:

${\frac {\partial ^{2}u}{\partial t^{2}}}=c^{2}{\frac {\partial ^{2}u}{\partial x^{2}}}$

em que $u(x,t)$ é o deslocamento vertical da corda e $0<x<L$, com $L$ o comprimento da corda.

Admitindo $c=1$,

${\frac {\partial }{\partial t}}{\Big (}{\frac {\partial u}{\partial t}}{\Big )}={\frac {\partial }{\partial x}}{\Big (}{\frac {\partial u}{\partial x}}{\Big )}$ .

Uma vez que os métodos citados abaixo são para equações de primeira ordem, é necessário separarmos a equação em um sistema de equações, fazendo a substituição $v={\frac {\partial u}{\partial t}}$ e $w={\frac {\partial u}{\partial x}}$:

${\begin{cases}{\dfrac {\partial v}{\partial t}}={\dfrac {\partial w}{\partial x}}\\\\{\dfrac {\partial v}{\partial x}}={\dfrac {\partial w}{\partial t}}\\\end{cases}}$

As condições de contorno são $u(0,t)=u(L,t)=0$ (pontas fixas), e as condições iniciais são $u(x,0)=\sin {\frac {\pi x}{L}}$ e ${\frac {\partial u}{\partial t}}(x,0)=0$

Algoritmos

Apresentaremos aqui três abordagens diferentes para a solução da equação diferencial parcial apresentada, e após, seus respectivos erros associados. A respeito das discretizações, $j$ corresponde à posição, e $n$ representa o tempo.

Método de Lax-Friedrichs - $O(\Delta x^{2},\Delta t)$

Esse método consiste em inicialmente discretizar as equações no esquema FTCS (Forward Time Centered Space), ou seja, discretizando a derivada temporal utilizando os tempos n e n+1 e a derivada espacial através das posições j-1 e j+1:

${\frac {1}{\Delta t}}(v_{j}^{n+1}-v_{j}^{n})={\frac {1}{2\Delta x}}(w_{j+1}^{n}-w_{j-1}^{n})$,

${\frac {1}{\Delta t}}(w_{j}^{n+1}-w_{j}^{n})={\frac {1}{2\Delta x}}(v_{j+1}^{n}-v_{j-1}^{n})$.

Resultando em

$v_{j}^{n+1}=v_{j}^{n}+{\frac {\Delta t}{2\Delta x}}(w_{j+1}^{n}-w_{j-1}^{n})$,

$w_{j}^{n+1}=w_{j}^{n}+{\frac {\Delta t}{2\Delta x}}(v_{j+1}^{n}-v_{j-1}^{n})$.

Entretanto, ao se realizar uma análise de estabilidade de Von Neumann, conclui-se que esse método é instável[3] . Para torná-lo estável, é necessário trocarmos os termos $v_{j}^{n}$ e $w_{j}^{n}$ por suas médias espaciais, chegando, assim, na expressão do esquema de Lax-Friedrichs

$v_{j}^{n+1}={\Big (}{\frac {v_{j-1}^{n}+v_{j+1}^{n}}{2}}{\Big )}+{\frac {\Delta t}{2\Delta x}}(w_{j+1}^{n}-w_{j-1}^{n})$,

$w_{j}^{n+1}={\Big (}{\frac {w_{j-1}^{n}+w_{j+1}^{n}}{2}}{\Big )}+{\frac {\Delta t}{2\Delta x}}(v_{j+1}^{n}-v_{j-1}^{n})$.

Para obtermos o valor de $u_{j}^{n+1}$, que é o nosso objetivo, discretizamos a equação ${\frac {\partial u}{\partial t}}=v$

$u_{j}^{n+1}=u_{j}^{n}+v_{j}^{n}\Delta t$,

Embora as médias espaciais sejam necessárias para a estabilidade do método, elas introduzem um problema: surge um efeito chamado de dissipação numérica, ou seja, a amplitude da solução diminui com o tempo. Isso pode ser observado através da análise de Von Neumann ou de uma investigação da equação do esquema Lax-Friedrichs [3]. Por este método, observa-se que ao inserirmos as médias, mudamos a equação original do problema, pois agora há também um termo do tipo difusivo (uma derivada segunda).

Agora vamos unir todas as equações, utilizando, além da equação para $v$ obtida acima, as discretizações de $v={\frac {\partial u}{\partial t}}$ e $w={\frac {\partial u}{\partial x}}$

$v_{j}^{n}={\frac {1}{\Delta t}}(u_{j}^{n+1}-u_{j}^{n})$,

$w_{j}^{n}={\frac {1}{2\Delta x}}(u_{j+1}^{n}-u_{j-1}^{n})$.

Assim, obtemos

$u_{j}^{n+1}=u_{j}^{n}+{\Big (}{\frac {u_{j-1}^{n}+u_{j+1}^{n}}{2}}{\Big )}-{\Big (}{\frac {u_{j-1}^{n-1}+u_{j+1}^{n-1}}{2}}{\Big )}+{\frac {(\Delta t)^{2}}{4(\Delta x)^{2}}}{\Big (}u_{j-2}^{n-1}-2u_{j}^{n-1}+u_{j+2}^{n-1}{\Big )}$.

Método de Leapfrog - $O(\Delta x^{2},\Delta t^{2})$

Neste método utilizamos os pontos intermediários da discretização.

Para $v$ temos

${\frac {1}{\Delta t}}{\bigg (}v_{j}^{n+{\frac {1}{2}}}-v_{j}^{n-{\frac {1}{2}}}{\bigg )}={\frac {1}{\Delta x}}{\Big (}w_{j+{\frac {1}{2}}}^{n}-w_{j-{\frac {1}{2}}}^{n}{\bigg )}\Rightarrow v_{j}^{n+{\frac {1}{2}}}=v_{j}^{n-{\frac {1}{2}}}+{\frac {\Delta t}{\Delta x}}{\bigg (}w_{j+{\frac {1}{2}}}^{n}-w_{j-{\frac {1}{2}}}^{n}{\bigg )}$,

Para $w$ temos

${\frac {1}{\Delta x}}{\bigg (}v_{j+1}^{n+{\frac {1}{2}}}-v_{j}^{n+{\frac {1}{2}}}{\bigg )}={\frac {1}{\Delta t}}{\bigg (}w_{j+{\frac {1}{2}}}^{n+1}-w_{j+{\frac {1}{2}}}^{n}{\bigg )}\Rightarrow w_{j+{\frac {1}{2}}}^{n+1}=w_{j+{\frac {1}{2}}}^{n}+{\frac {\Delta t}{\Delta x}}{\bigg (}v_{j+1}^{n+{\frac {1}{2}}}-v_{j}^{n+{\frac {1}{2}}}{\bigg )}$,

Para $u$ temos

$v_{j}^{n+{\frac {1}{2}}}={\frac {\partial u}{\partial t}}{\bigg |}_{j}^{n+{\frac {1}{2}}}={\frac {1}{\Delta t}}(u_{j}^{n+1}-u_{j}^{n})\Rightarrow u_{j}^{n+1}=u_{j}^{n}+v_{j}^{n+{\frac {1}{2}}}\Delta t$,

Utilizando o fato de que

$v_{j}^{n+{\frac {1}{2}}}={\frac {1}{\Delta t}}(u_{j}^{n+1}-u_{j}^{n})$,

$w_{j+{\frac {1}{2}}}^{n}={\frac {1}{\Delta x}}(u_{j+1}^{n}-u_{j}^{n})$,

chegamos na equação para $u_{j}^{n+1}$

$u_{j}^{n+1}=2u_{j}^{n}-u_{j}^{n-1}+{\frac {(\Delta t)^{2}}{(\Delta x)^{2}}}(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n})$,

o que é equivalente a discretizarmos a equação da onda diretamente, utilizando que, para uma função $f(x)$,

${\frac {d^{2}f}{dx^{2}}}{\bigg |}_{j}={\frac {1}{(\Delta x)^{2}}}(f_{j-1}-2f_{j}+f_{j+1})$,

sendo $j$ a discretização em $x$.

Método de Lax-Wendroff de Dois Passos - $O(\Delta x^{2},\Delta t^{2})$

O primeiro passo consiste em calcular o valor de $v^{n+{\frac {1}{2}}}$ e $w^{n+{\frac {1}{2}}}$ utilizando o método de Lax-Friedrichs, para posterior cálculo de $v^{n+1}$ e $w^{n+1}$:

$v_{j+{\frac {1}{2}}}^{n+{\frac {1}{2}}}={\frac {1}{2}}(v_{j+1}^{n}+v_{j}^{n})+{\frac {\Delta t}{2\Delta x}}(w_{j+1}^{n}-w_{j}^{n})$,

$v_{j-{\frac {1}{2}}}^{n+{\frac {1}{2}}}={\frac {1}{2}}(v_{j}^{n}+v_{j-1}^{n})+{\frac {\Delta t}{2\Delta x}}(w_{j}^{n}-w_{j-1}^{n})$,

$w_{j+{\frac {1}{2}}}^{n+{\frac {1}{2}}}={\frac {1}{2}}(w_{j+1}^{n}+w_{j}^{n})+{\frac {\Delta t}{2\Delta x}}(v_{j+1}^{n}-v_{j}^{n})$,

$w_{j-{\frac {1}{2}}}^{n+{\frac {1}{2}}}={\frac {1}{2}}(w_{j}^{n}+w_{j-1}^{n})+{\frac {\Delta t}{2\Delta x}}(v_{j}^{n}-v_{j-1}^{n})$,

Agora, no tempo $n+1$:

$v_{j}^{n+1}=v_{j}^{n}+{\frac {\Delta t}{\Delta x}}{\bigg (}w_{j+{\frac {1}{2}}}^{n+{\frac {1}{2}}}-w_{j-{\frac {1}{2}}}^{n+{\frac {1}{2}}}{\bigg )}$,

$w_{j}^{n+1}=w_{j}^{n}+{\frac {\Delta t}{\Delta x}}{\bigg (}v_{j+{\frac {1}{2}}}^{n+{\frac {1}{2}}}-v_{j-{\frac {1}{2}}}^{n+{\frac {1}{2}}}{\bigg )}$,

Agrupando as equações,

$w_{j}^{n+1}=w_{j}^{n}+{\frac {\Delta t}{\Delta x}}{\Bigg [}{\frac {1}{2}}(v_{j+1}^{n}-v_{j-1}^{n})+{\frac {\Delta t}{2\Delta x}}(w_{j+1}^{n}-2w_{j}^{n}+w_{j-1}^{n}){\Bigg ]}$,

$v_{j}^{n+1}=v_{j}^{n}+{\frac {\Delta t}{\Delta x}}{\Bigg [}{\frac {1}{2}}(w_{j+1}^{n}-w_{j-1}^{n})+{\frac {\Delta t}{2\Delta x}}(v_{j+1}^{n}-2v_{j}^{n}+v_{j-1}^{n}){\Bigg ]}$,

E finalmente temos a equação unificada em u, utilizando a expressão para $v_{j}^{n+1}$ e as discretizações de $v={\frac {\partial u}{\partial t}}$ e $w={\frac {\partial u}{\partial x}}$, como obtidas na seção sobre o Método de Lax-Friedrichs:

$u_{j}^{n+1}=u_{j}^{n}+v_{j}^{n}\Delta t$

$u_{j}^{n+1}=2u_{j}^{n}-u_{j}^{n-1}+{\frac {(\Delta t)^{2}}{2(\Delta x)^{2}}}{\Bigg [}{\frac {1}{2}}u_{j+2}^{n-1}-{\frac {1}{2}}u_{j-2}^{n-1}+u_{j+1}^{n}-u_{j+1}^{n-1}-2u_{j}^{n}+u_{j}^{n-1}+u_{j-1}^{n}-u_{j-1}^{n-1}{\Bigg ]}$,

Programas

Ao implementarmos o método, surgem dois problemas: o problema não é auto-inicializável, pois para calcularmos o valor de $u_{j}^{1}$, necessitamos de $u_{j}^{-1}$ (além de $u_{j}^{0}$). Entretanto, isto é rapidamente solucionado quando discretizamos a condição inicial de que ${\frac {\partial u}{\partial t}}(x,0)=0$:

${\frac {1}{2\Delta t}}(u_{j}^{1}-u_{j}^{-1})=0\Rightarrow u_{j}^{1}=u_{j}^{-1}$,

ou seja, para o cálculo de $u_{j}^{1}$, utilizamos que $u_{j}^{1}=u_{j}^{-1}$. Para o método de Leapfrog, dessa forma conseguimos isolar $u_{j}^{1}$, não havendo nenhum outro termo no tempo $n=1$:

$u_{j}^{-1}=2u_{j}^{0}-u_{j}^{-1}+{\frac {(\Delta t)^{2}}{(\Delta x)^{2}}}(u_{j+1}^{0}-2u_{j}^{0}+u_{j-1}^{0})$,

$u_{j}^{-1}=u_{j}^{0}+{\frac {(\Delta t)^{2}}{2(\Delta x)^{2}}}(u_{j+1}^{0}-2u_{j}^{0}+u_{j-1}^{0})$.

Porém, isso não ocorre com os outros dois métodos, pois surgem termos em diferentes posições para o tempo $n=0$ (de $u_{-2}^{0}$, $u_{-1}^{0}$, até $u_{2}^{0}$), sendo necessário resolvermos o sistema como um todo simultaneamente, ou seja, teríamos que inverter uma matriz. Por isso, foi utilizado o método de Leapfrog para o cálculo de $u_{j}^{-1}$ em todos os métodos, devido a sua simplicidade.

Além disso, são necessários valores de $u_{-1}^{n}$ e de $u_{j_{max}+1}^{n}$, com $j_{max}$ correspondendo a $x=L$, para calcularmos $u_{1}^{n}$ e $u_{j_{max}}^{n}$, para qualquer tempo, utilizando os métodos de Lax-Wendroff de dois passos e Lax-Friedrichs. A solução a este problema foi utilizarmos

$u_{-1}^{n}=-u_{1}^{n}$

$u_{j_{max}+1}^{n}=-u_{j_{max}-1}^{n}$.

Pensando na condição inicial $u(x,0)=\sin {\frac {\pi x}{L}}$, e estendendo para além da corda (pensando no seno de $-\infty <x<\infty$), observamos que ela respeita as equações acima.

Análise de erros

A análise de erros se torna mais evidente durante a escolha do parâmetro $k$, onde $k={\frac {\Delta t}{\Delta x}}$. Valores grandes trazem pouca acurácia, e valores pequenos necessitam de muito poder de computação (tempo e dinheiro). Trazemos problemas mais simplificados como um "guia" de escolha do parâmetro.

A partir do cálulo da solução analítica da equação da onda, podemos calcular quanto o valor obtido pelos métodos difere da solução real, o que leva a uma visualização do erro corrente em cada método de integração. Nesse caso, a solução é $u(x,t)=\cos {{\Big (}{\frac {\pi t}{L}}{\Big )}}\sin {{\Big (}{\frac {\pi x}{L}}{\Big )}}$.

Podemos observar a ordem com que os erros crescem à medida que o parâmetro k se torna maior. Lembrando que os valores da constante são determinados pela discretização do espaço e do tempo.

Aqui variamos o valor de $\Delta t$, fixando $\Delta x=1$

Comparação do erro global entre os três métodos estudados

Observamos que se determinarmos a reta que melhor se ajusta às curvas dos métodos de Leapfrog e Lax-Wendroff, ela tem inclinação aproximadamente 1, já que a inclinação dá a ordem do erro global (já que fizemos o cálculo do erro após muitos passos transcorridos), que é sempre uma ordem menor do que o erro local (erro de um passo).

Simulação de Propagação de Onda 2D Dependente de Topografia

O modelo mais simples parte da equação da onda [1], acrescentando o termo $H(x,y,t)$.

${\frac {\partial ^{2}u}{\partial t^{2}}}={\Big (}{\frac {\partial }{\partial x}}H(x,y,t){\frac {\partial u}{\partial x}}{\Big )}+{\Big (}{\frac {\partial }{\partial y}}H(x,y,t){\frac {\partial u}{\partial y}}{\Big )}-{\frac {\partial ^{2}H}{\partial t^{2}}}$ ,

Sendo $H(x,y,t)$ uma representação da profundidade em águas calmas. Em uma situação real, pode-se obtê-la por mapeamento eletrônico do terreno por sistema de sonar.

Exemplo de mapeamento de terreno sub - calota polar feito por AUV (autonomous underwater vehicle)[2]

Como primeira abordagem visando uma análise em 2D, a integração da equação em 1D (mesmo sendo uma situação muito idealizada) já traz resultados interessantes. Podemos observar, por exemplo, que a amplitude da onda cresce perto da costa. Esta informação por si só ajuda na construção de proteção contra quebra de ondas, pois é obtido o tamanho que as mesmas atingem.

Simulação em 1D de ondas perto da margem

É importante notar o quão poderosa é a integração de equações parciais na vida de um engenheiro.

A dependência em $t$ de $H(x,y,t)$ permite um modelo no qual o terreno se modifica com o tempo. Isto é, pode-se observar o efeito que o deslocamento de placas tectônicas, deslizamentos, e até explosões provocam no comportamento das ondas na costa de um país e o reconhecimento de áreas críticas.

Estendendo o algoritmo do Leap-Frog à situação 2D, obtemos, para uma dada condição inicial e $H(x,y,t)=C$, onde $C$ é uma constante:

Simulação em 2D de ondas em águas com profundidade constante, visão de cima

Simulação em 2D de onda em águas com profundidade constante, visão em ângulo

Podemos então, analisar como a mesma condição inicial se porta quando $H(x,y,t)$ descreve uma gaussiana na origem:

Simulação em 2D de ondas em águas com gaussiana na origem, visão de cima

Simulação em 2D de onda em águas com gaussiana na origem, visão em ângulo

Bibliografia

¹"The Wave Equation in 1D and 2D," por Knut–Andreas Lie, Dept. of Informatics, University of Oslo; disponível em: [1]; Último acesso em 23/10/2017.

²"Digital terrain mapping of the underside of sea ice from a small AUV," por Wadhams, M. J. Doble; disponível em: DOI: 10.1029/2007GL031921 ; Último acesso em 23/10/2017.

³ Press, William H.; Teukolsky, Saul A.; Vetterling, William T.; Flannery, Brian P. (2007). Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8.

⁴ "Métodos Computacionais da Física C - Aula 2 - 2016/1" por Heitor C. M. Fernandes, Instituto de Física, UFRGS; Último acesso em 23/10/2017.