Métodos Computacionais para Estudo de Fenômenos Ondulatórios

Filssen Schereiber, João Roth e Lucas Oliveira

Este trabalho tem como objetivo realizar um estudo introdutório da equação da onda, adotando uma abordagem computacional que proporciona uma base teórica sólida e intuitiva. Inicialmente, é realizada a implementação de uma solução numérica (FTCS) para a equação da advecção, visando o estudo da propagação de uma perturbação e explorando aspectos computacionais fundamentais. Em seguida, a equação da onda é aplicada em uma dimensão, permitindo uma investigação detalhada de fenômenos ondulatórios, tais como reflexão, interferência e refração. A abordagem é estendida para duas dimensões, incluindo a análise dos efeitos de contorno nesse contexto. Por fim, o trabalho aprofunda-se no estudo do fenômeno de difração, contribuindo para uma compreensão mais enriquecedora dos aspectos complexos dos fenômenos ondulatórios em múltiplas dimensões.

Introdução

O que é uma onda exatamente? De acordo com 'A Student’s Guide to Waves', a literatura apresenta diversas definições. Contudo, a característica mais comum que define uma onda é sua natureza como uma perturbação, uma alteração no estado de equilíbrio que, inicialmente, permanece inalterado. Quando uma fonte externa, como um objeto vibrante ou uma força inicial, perturba a condição de equilíbrio de um meio, essa perturbação é transmitida de partícula para partícula ao longo do meio. Esse processo de transmissão de energia ocorre sem que as partículas individuais do meio se desloquem significativamente de suas posições de equilíbrio. Em outras palavras, a energia da perturbação é transferida, mas as partículas do meio não se movem em conjunto com a onda.

O estudo da física das ondas desempenha um papel crucial, permeando inúmeras áreas da ciência e da tecnologia ao proporcionar uma compreensão profunda dos fenômenos ondulatórios e suas aplicações práticas. A física das ondas desempenha um papel essencial em nossa vida cotidiana. Seu impacto é evidente na tecnologia de comunicação, onde as ondas de rádio e micro-ondas são a espinha dorsal de redes globais de telefonia móvel e internet sem fio.

Com base no exposto, propõe-se a realização de um estudo abrangente sobre a modelagem de fenômenos ondulatórios. Inicialmente, serão apresentados pontos importantes acerca da equação da onda. Em seguida, dedicaremos atenção ao estudo da propagação de uma perturbação a partir da equação de advecção. Nesse contexto, abriremos a discussão sobre a abordagem computacional, e, tangencialmente, comentaremos sobre outros métodos numéricos para resolver a advecção. Posteriormente, faremos um estudo sobre a solução numérica da equação da onda, permitindo a modelagem de fenômenos ondulatórios, como reflexão, interferência, refração e difração.

Equação da Onda

A função de uma onda é aquela que especifica o valor da perturbação em cada ponto e instante ao longo de seu percurso. Para iniciar nossa análise, abordaremos o exemplo mais simples, em que a propagação ocorre unicamente em uma direção, como é o caso das ondas transversais em uma corda.

Dedução

O perfil da onda em uma corda em um determinado instante é equivalente à forma que a corda apresenta nesse momento, que é dada pela função $y(x,t)$. A perturbação assume a forma de uma onda progressiva, movendo-se como uma entidade coesa para a direita, mantendo sua configuração inalterada, com velocidade $v$. Dessa forma, ao acompanhar a onda em um referencial inercial diferente, em que $y(x,t)=y(x',t)$, onde $x'$ representa o antigo $x$ deslocado por $vt$, a relação entre os dois referenciais é estabelecida por

$x'=x-vt$. $y'=y$.

de modo que, no referencial original,

$y(x,t)=f(x-vt)$.

descreve uma onda progressiva, que se propaga para a direita, com velocidade $v$.

Para associar uma equação de movimento com a propagação da onda, vamos calcular a aceleração num dado ponto $x$. A velocidade e a aceleração em $x$ se obtêm fixando $x$ e derivando em relação ao tempo, o que corresponde a tomar derivadas parciais.

$V={\frac {\partial }{\partial t}}y(x,t)$. $A={\frac {\partial ^{2}}{\partial t^{2}}}y(x,t)$.

pela regra da cadeia

${\frac {\partial y}{\partial t}}={\frac {\partial f}{\partial x'}}{\frac {\partial x'}{\partial t}}=-v{\frac {\partial f}{\partial x'}}$.

Analogamente,

${\frac {\partial ^{2}}{\partial t^{2}}}=-v{\frac {\partial }{\partial t}}{\frac {\partial f}{\partial x'}}=-v{\frac {\partial }{\partial x'}}{\frac {\partial f}{\partial x'}}{\frac {\partial x'}{\partial t}}$.

ou seja,

${\frac {\partial ^{2}y}{\partial t^{2}}}=v^{2}{\frac {\partial ^{2}f}{\partial {x'}^{2}}}$.

Derivando $y$ em relação a $x$ obtemos que:

${\frac {\partial y}{\partial x}}={\frac {\partial f}{\partial x'}}$. ${\frac {\partial ^{2}y}{\partial x^{2}}}={\frac {\partial ^{2}f}{\partial {x'}^{2}}}$.

Assim,

${\frac {\partial ^{2}y}{\partial t^{2}}}-v^{2}{\frac {\partial ^{2}y}{\partial x^{2}}}=0$.

Denominada equação da onda unidimensional, esta é uma das equações fundamentais da física.

Solução Analítica

Agora, procederemos à análise do movimento de uma corda de comprimento $L$ que vibra com extremidades fixas. Como essas extremidades são fixas, temos, nos pontos $x=0$ e $x=L$, os valores $y(0,t)=y(L,t)=0$.

Além disso, supomos que a função $y(x,0)=f(x)$ e ${\frac {\partial y}{\partial t}}=0$ para $t=0$.

Com base nas observações realizadas, torna-se possível identificar um conjunto fundamental de soluções por meio da aplicação do método de separação de variáveis.

Suponhamos que existam duas soluções, uma dependendo exclusivamente de $x$ e outra exclusivamente de $t$, de modo que possamos expressar a função $y(x,t)$ como o produto dessas duas.

$y(x,t)=X(x)T(t)$

Assim, podemos escrever as derivadas parciais como

${\begin{cases}{\frac {\partial ^{2}y}{\partial x^{2}}}={\frac {X''}{X}}\\{\frac {\partial ^{2}y}{\partial t^{2}}}={\frac {T''}{T}}\end{cases}}$

substituindo na equação da onda,

${\frac {T''}{T}}-v^{2}{\frac {X''}{X}}=0$

Dessa forma, torna-se evidente que a razão entre as funções de <mathT</math> consiste em uma constante, assim como $v^{2}{\frac {X''}{X}}$.

[*Equação (12)*]

$v^{2}{\frac {X''}{X}}=\lambda ^{2}\to X''={\bigg (}{\frac {\lambda }{v}}{\bigg )}^{2}X$

[*Equação (13)*]

${\frac {T''}{T}}=\lambda ^{2}\to T''=\lambda ^{2}T$

Observa-se que as equações (12) e (13) são equações diferenciais ordinárias (EDOs) lineares de segunda ordem, cujos coeficientes podem depender da variável independente e, potencialmente, serem funções não constantes. Assim, é possível usarmos da teoria de Sturm-Liouville para resolver essas equações.

Resolvendo a Equação (12)

Se avaliarmos as condições de contorno para os possíveis valores de lambda, que devem ser maiores ou iguais a zero, obtemos apenas a solução trivial. Contudo, se lambda for negativo, surge uma solução para $X(x)$ na forma

$X(x)=a\cos {\bigg (}{\frac {\lambda x}{v}}{\bigg )}+b\sin {\bigg (}{\frac {\lambda x}{v}}{\bigg )}$

Avaliando as condições de contorno, temos:

$X(0)=a\cos(0)+b\sin(0)=a\therefore a=0$ $X(L)=b\sin {\bigg (}{\frac {\lambda L}{v}}{\bigg )}$

Certamente, observa-se que, para alcançarmos uma solução não trivial, é necessário que $\lambda ={\bigg (}{\frac {nv\pi }{L}}{\bigg )}$, onde $n\in \mathbb {N}$.

$X_{n}(x)=\sum _{n=1}^{\infty }b_{n}\sin {\bigg (}{\frac {n\pi x}{L}}{\bigg )}$

Resolvendo a Equação (13)

Analogamente, temos uma solução não trivial para lambda menor que zero.

$T(t)=c\cos(\lambda t)+d\sin(\lambda t)\to T_{n}(t)=\sum _{n=1}^{\infty }c_{n}\cos {\bigg (}{\frac {nv\pi t}{L}}{\bigg )}+d_{n}\sin {\bigg (}{\frac {nv\pi t}{L}}{\bigg )}$

Voltando para a solução geral:

$y(x,t)=\sum _{n=1}^{\infty }b_{n}\sin {\bigg (}{\frac {n\pi x}{L}}{\bigg )}{\bigg (}c_{n}\cos {\bigg (}{\frac {nv\pi t}{L}}{\bigg )}+d_{n}\sin {\bigg (}{\frac {nv\pi t}{L}}{\bigg )}{\bigg )}$

Juntando as constantes

$y(x,t)=\sum _{n=1}^{\infty }\sin {\bigg (}{\frac {n\pi x}{L}}{\bigg )}{\bigg (}A_{n}\cos {\bigg (}{\frac {nv\pi t}{L}}{\bigg )}+B_{n}\sin {\bigg (}{\frac {nv\pi t}{L}}{\bigg )}{\bigg )}$

Avaliando as condições iniciais:

Para $t=0$, temos

$y(x,0)=\sum _{n=1}^{\infty }A_{n}\sin {\bigg (}{\frac {n\pi x}{L}}{\bigg )}=f(x)$

É possível observar que $f(x)$ representa uma série de Fourier de senos, na qual os coeficientes $A_{n}$ são determinados por:

$A_{n}={\frac {2}{L}}\int _{0}^{L}f(x)\sin {\bigg (}{\frac {n\pi x}{L}}{\bigg )}dx$

Avaliando $y'(x,0)$:

${\frac {\partial }{\partial t}}y(x,0)=\sum _{n=1}^{\infty }B_{n}{\bigg (}{\frac {nv\pi }{L}}{\bigg )}\sin {\bigg (}{\frac {n\pi x}{L}}{\bigg )}=g(x)$

Analogamente a $f(x)$, $g(x)$ é uma série de Fourier de senos. Portanto, podemos determinar $B_{n}$ da seguinte maneira:

$B_{n}={\frac {2}{nv\pi }}\int _{0}^{L}g(x)\sin {\bigg (}{\frac {n\pi x}{L}}{\bigg )}dx$

A equação de advecção unidimensional

Partimos da equação de advecção para estudarmos a propagação de uma perturbação.

A equação de onda está intimamente relacionada à chamada equação de advecção, que, em uma dimensão, assume a seguinte forma:

${\frac {\partial u}{\partial t}}+v{\frac {\partial u}{\partial x}}=0$

(4)

A equação a seguir descreve a advecção passiva de um campo escalar $u(x,t)$ transportado por um fluxo de velocidade constante $v$. Como a equação de advecção é um tanto mais simples do que a equação da onda, discutiremos aquela em primeiro lugar. A equação de advecção possui a solução formal:

$u(x,t)=F(x-vt)$

(5)

Essa solução descreve um pulso de forma arbitrária que é arrastado pelo fluxo a uma velocidade constante $v$, sem alterar sua forma.

Esquema Forward Time-Centered Space (FTCS)

Para a abordagem numérica computacional do nosso estudo, utilizamos o esquema FTCS.

O método FTCS é amplamente empregado na resolução de equações diferenciais parciais. Trata-se de um método de diferenças finitas que adota uma abordagem explícita no tempo e central no espaço. Esse esquema apresenta uma acurácia de primeira ordem no tempo e de segunda ordem no espaço.

A discretização no tempo (Forward Time) pode ser feita utilizando a diferença finita para a derivada temporal:

${\frac {\partial u}{\partial t}}\approx {\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}$

(6)

onde $u_{i}^{n}$ representa a amplitude da perturbação no ponto $i$ do espaço e no instante $n$, e $\Delta t$ é o passo de tempo.

A discretização no espaço (Central Space) pode ser feita utilizando a diferença finita para a derivada espacial:

${\frac {\partial u}{\partial x}}\approx {\frac {u_{i+1}^{n}-u_{i-1}^{n}}{2\Delta x}}$

(7)

onde $\Delta x$ é o espaçamento entre os pontos discretizados no espaço.

As equações acima podem ser reescritas como

$u_{i}^{n+1}=u_{i}^{n}+{\frac {r}{2}}(u_{i+1}^{n}-u_{i-1}^{n})$

(8)

onde $r=v{\frac {\Delta t}{\Delta x}}$.

Análise de estabilidade de von Neumann

A análise de estabilidade de von Neumann é uma abordagem teórica utilizada para avaliar a estabilidade numérica de esquemas de diferenças finitas, incluindo o esquema FTCS. Usando o método de von Neumann, supõe-se a seguinte solução tentativa para a EDP:

$u(x,t)=A(t)e^{jkx}$

(8)

onde $j={\sqrt {-1}}$, $k$ é o número de onda, e $A$ é a amplitude da onda, uma função do tempo.

Assim, a solução no passo de tempo $n$ e em $x=x_{i}=i\Delta x$ é escrita como

$A_{n}e^{jki\Delta x}$

(9)

Substituindo esta solução tentativa na expressão da advecção discretizada temos:

$A_{n+1}e^{jki\Delta x}=A_{n}e^{jki\Delta x}-{\frac {u\Delta t}{2\Delta x}}\left(A_{n}e^{jk(i+1)\Delta x}-A_{n}e^{jk(i-1)\Delta x}\right)$

(10)

Seja $\xi$ a razão da amplitude da onda no passo de tempo $n+1$ em relação ao passo $n$. Portanto,

$\xi ={\frac {A_{n+1}}{A_{n}}}$

(11)

Dividir (10) por $A_{n}$ resulta em

$\xi e^{jki\Delta x}=e^{jki\Delta x}-{\frac {u\Delta t}{2\Delta x}}\left(e^{jk(i+1)\Delta x}-e^{jk(i-1)\Delta x}\right)$

(12)

Se dividirmos esse resultado por $e^{jki\Delta x}$ temos:

$\xi =1-{\frac {u\Delta t}{2\Delta x}}\left(e^{jk\Delta x}-e^{-jk\Delta x}\right)=1-{\frac {u\Delta t}{\Delta xj}}\sin(k\Delta x)$

(13)

Assim,

$|\xi |={\sqrt {1+\left({\frac {u\Delta t}{\Delta x}}\sin(k\Delta x)\right)^{2}}}$

(14)

Isso implica que $|\xi |\geq 1$ independentemente do passo de tempo $\Delta t$ selecionado ou do passo espacial $\Delta x$, portanto, FTCS é incondicionalmente instável.

Apesar de ser incondicionalmente instável, o esquema funcionará bem em nossa modelagem, uma vez que o método apresenta também suas vantagens. O esquema FTCS é relativamente simples de implementar, o que facilita sua aplicação, principalmente para problemas bidimensionais ou tridimensionais, como pretendemos abordar posteriormente.

Em alguns casos específicos, o FTCS pode ser computacionalmente eficiente. Se a equação diferencial parcial for relativamente simples e o domínio de interesse for bem comportado, o FTCS pode oferecer uma solução razoável com custo computacional relativamente baixo.

É feita uma comparação entre outros métodos numéricos de diferenças finitas para solução da equação da advecção e comparado o erro relativo à solução analítica no Apêndice A.

Implementação do esquema FTCS

Vamos agora examinar a implementação do esquema FTCS na equação da advecção para visualizarmos a propagação de uma perturbação na forma de $1-\cos(x)$ em um eixo $x$ infinito. Em outras palavras:

• Condição inicial: $u(x,0)=1-\cos(x)$;

• Condições de contorno periódicas.

Neste caso, podemos observar uma perturbação se propagando no sentido negativo do eixo $x$ sem perda de amplitude ou mudança de forma. Isso acontece porque, com a equação da advecção, podemos analisar apenas a cinemática de uma perturbação, a qual, por sua vez, apenas relaciona os atributos como amplitude, comprimento de onda, frequência e velocidade.

Podemos observar, ainda, que a amplitude está relacionada à condição inicial em que o máximo de $1-\cos(x)$ é atingido quando $\cos(x)=-1$. Nesse contexto, considerando $0<x<2\pi$, percebemos que $\cos(x)=-1$ ocorre exatamente no meio do eixo. Em consequência, a amplitude máxima, que é 2, é inicialmente alcançada no ponto central do eixo $x$.

A equação da onda em uma dimensão

A equação da onda, que em uma dimensão assume a forma

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

(15)

onde $u$ geralmente é algum tipo de deslocamento ou perturbação, enquanto $v$ é a velocidade da onda. A equação de onda possui a solução formal

$u(x,t)=F(x+vt)+G(x+vt),$

(16)

onde $F$ e $G$ são funções arbitrárias. A solução acima representa pulsos de onda de formato arbitrário que se propagam com velocidade $v$ nas direções $+x$ e $-x$, respectivamente, sem alterar a forma.

Esquema FTCS

Iniciamos o estudo da solução numérica da equação da onda a partir da implementação do esquema FTCS, pois, como comentado antes, este é o método de implementação mais prático para generalizarmos para mais dimensões espaciais.

As derivadas de segunda ordem podem ser substituídas por diferenças centrais. A aproximação de diferença mais amplamente utilizada para a derivada de segunda ordem é

${\frac {\partial ^{2}}{\partial t^{2}}}u(x_{i},t_{n})\approx {\frac {u_{i}^{n+1}-2u_{i}^{n}+u_{i}^{n-1}}{\Delta t^{2}}}$

(17)

Uma aproximação semelhante para a derivada de segunda ordem na direção $x$ é dada por

${\frac {\partial ^{2}}{\partial x^{2}}}u(x_{i},t_{n})\approx {\frac {u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{\Delta x^{2}}}$

(18)

Agora podemos substituir as derivadas e obter

${\frac {u_{i}^{n+1}-2u_{i}^{n}+u_{i}^{n-1}}{\Delta t^{2}}}=v^{2}{\frac {u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{\Delta x^{2}}}.$

(19)

Assumimos que $u_{i}^{n}$ e $u_{i}^{n-1}$ já foram calculados para $i=0,\ldots ,N_{x}$. A única quantidade desconhecida é, portanto, $u_{i}^{n+1}$, que podemos resolver da seguinte forma:

$u_{i}^{n+1}=-u_{i}^{n-1}+2u_{i}^{n}+r^{2}\left(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}\right),$

(20)

onde introduzimos o parâmetro

$r=v{\frac {\Delta t}{\Delta x}},$

conhecido como número de Courant.

Análise de estabilidade de von Neumann

(ALCINEY DAS NEVES MORAES)

Pela análise de Von Neumann, substituímos $u_{n}$ por $g_{n}e^{i\theta }$ na equação da onda. Após manipulações, a expressão resultante é:

$g=2-g^{-1}+{\frac {C}{2}}(e^{i\theta }-2+e^{-i\theta })$

(21)

Usando a identidade de Euler $e^{i\theta }+e^{-i\theta }=2\cos \theta$, chegamos a:

$g^{2}-2[1-{\frac {C}{2}}(1-\cos \theta )]g+1=0$

(21)

Reescrevendo $\cos \theta =\cos ^{2}\left({\frac {\theta }{2}}\right)-\sin ^{2}\left({\frac {\theta }{2}}\right)$, a equação torna-se:

$g^{2}+2(2C^{2}\sin ^{2}\left({\frac {\theta }{2}}\right)-1)g+1=0$

(22)

Expressamos isso como $g^{2}+2bg+1=0$, onde $b=2C^{2}\sin ^{2}\left({\frac {\theta }{2}}\right)-1$. As raízes são:

$g={\begin{cases}-b\pm {\sqrt {b^{2}-1}},&{\text{se }}|b|>1\\-b\pm i{\sqrt {1-b^{2}}},&{\text{se }}|b|\leq 1\end{cases}}$

(23)

Se $|b|>1$, implica $C\leq 1$. Se $|b|\leq 1$, $C\leq 1$ pela condição de Von Neumann, tornando o método condicionalmente estável pela condição CFL.

Reflexão de onda em uma corda

Quando uma onda se propaga ao longo de uma corda e incide em uma de suas extremidades, ocorre a sua reflexão. Ao atingir um extremo fixo, que não pode oscilar, é possível observar que o pulso refletido apresenta uma orientação oposta à do pulso incidente. Podemos afirmar que o pulso reflete com inversão de fase, pois isso ocorre devido à reação do ponto fixo sobre a corda com uma força de mesmo módulo. Uma situação distinta ocorre quando a extremidade da corda é móvel, ou seja, não há troca de forças entre ela e a corda, resultando em uma reflexão sem inversão de fase.

Além disso, em um sistema fechado e sem dissipação de energia, podemos observar que as ondas refletidas, independentemente da fase, sofrem interferência construtiva e produzem picos de máximo ou mínimo, dependendo do contexto.

Podemos observar esse fenômeno ao implementar a solução da equação da onda unidimensional com condições de contorno fixas e periódicas, respectivamente, para uma corda com extremidades fixas e extremidades móveis.

A equação da onda em duas dimensões

FTCS

Seja a equação da onda em duas dimensões dada por:

${\frac {\partial ^{2}u}{\partial t^{2}}}=v^{2}\left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)$

A discretização desta equação utilizando o esquema FTCS envolve as seguintes aproximações:

${\frac {\partial ^{2}u}{\partial t^{2}}}\approx {\frac {u_{i,j}^{n+1}-2u_{i,j}^{n}+u_{i,j}^{n-1}}{\Delta t^{2}}}$

${\frac {\partial ^{2}u}{\partial x^{2}}}\approx {\frac {u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}}{\Delta x^{2}}}$

${\frac {\partial ^{2}u}{\partial y^{2}}}\approx {\frac {u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{\Delta y^{2}}}$

Assim, a equação da onda discretizada usando FTCS é dada por:

$u_{i,j}^{n+1}=2u_{i,j}^{n}-u_{i,j}^{n-1}+v^{2}\left({\frac {\Delta t^{2}}{\Delta x^{2}}}\cdot (u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n})+{\frac {\Delta t^{2}}{\Delta y^{2}}}\cdot (u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n})\right)$

Para simplificar a notação, fazemos $R1=v^{2}{\frac {\Delta t^{2}}{\Delta x^{2}}}$ e $R2=v^{2}{\frac {\Delta t^{2}}{\Delta y^{2}}}$.

$u_{i,j}^{n+1}=2u_{i,j}^{n}-u_{i,j}^{n-1}+R1\cdot (u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n})+R2\cdot (u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n})$

Reflexão de ondas em duas dimensões

Analogamente ao que analisamos para uma perturbação em uma corda, podemos agora avaliar como uma perturbação evolui ao se propagar em um plano em direção às suas extremidades, e como interage com os contornos, que podem ser fixos ou periódicos. Neste caso, temos uma perturbação senoidal centrada em $x=0,5$ e $y=0,5$, propagando-se isotropicamente em direção aos limites $0$ e $1$ dos dois eixos.

Apêndice A: Métodos numéricos para solução da equação de advecção

Método Lax-Friedrich

Método Lax-Wendroff

Upwind Differencing

Método Leapfrog

Método Lax-Wendroff de dois passos

Comparação dos métodos numéricos