Grupo3 - Ondas2

Introdução

Equações diferenciais parciais (EDP's) hiperbólicas geralmente podem ser formuladas a partir de teoremas de conservação. Um exemplo é a equação do tipo:

$\frac{\partial \overrightarrow{U}}{\partial t}+\mathrm{\nabla }\cdot \overrightarrow{F}(U)=\overrightarrow{S}(U)$,

onde $\overrightarrow{U}$ é o vetor de densidades da quantidade conservada, i.e., $\overrightarrow{U}=(U_1,...,U_n)$, $\overrightarrow{F}$ é o fluxo de densidade e $\overrightarrow{S}$ é um termo genérico representando fontes ou sumidouros.

Uma classe especial de equações hiperbólicas são as chamadas equações de adveção, na qual a derivada temporal da quantidade conservada $\overrightarrow{U}(\overrightarrow{x},t)$ é proporcional à sua derivada espacial. Nesses casos, $\overrightarrow{F}(U)$ é diagonal e dada por:

$\overrightarrow{F}(U)=v\overrightarrow{I}\cdot \overrightarrow{U}$,

onde $\overrightarrow{I}$ é a matriz identidade.

Considerando apenas uma dimensão e com $\overrightarrow{U}\equiv u$, temos a equação de adveção:

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

onde $v$ é a velocidade de propagação do pulso gerado. A equação admite uma solução analítica da forma $u=f(x-vt)$, representando uma onda se movendo na direção $x$.

A equação da onda em uma dimensão é uma EDP hiperbólica de segunda ordem dada por

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

E admite duas soluções, representadas por pulsos, $f(x+vt)$ e $f(x-vt)$.

Assumindo que $v\ne v(x)$ na equação da onda, nos restringimos a problemas lineares. Além disso, se escrevermos

$k=v\frac{\partial u}{\partial x},s=\frac{\partial u}{\partial t},$

então a equação da onda pode ser escrita como um sistema de três equações diferenciais de primeira ordem:

$\{\begin{array}{l}\frac{\partial k}{\partial t}=v\frac{\partial s}{\partial x}\\ \frac{\partial s}{\partial t}=v\frac{\partial k}{\partial x}\\ \frac{\partial u}{\partial t}=s\end{array}$

Em notação vetorial, o sistema acima pode ser reescrito na forma conservativa como: $\frac{\partial U}{\partial t}+\frac{\partial F(U)}{\partial x}=0$,

onde $U=\left(\begin{array}{c}k\\ s\end{array}\right),\text{e}F(U)=\left(\begin{array}{cc}0& -v\\ -v& 0\end{array}\right)$

Condição CFL

Uma condição necessária mas nem sempre suficiente para a convergência de equações diferenciais parciais resolvidas a partir de métodos de diferença, formulada por Courant, Friedrichs e Lewy em 1928 (aprender a referenciar), conhecida como condição CFL, é formulada a partir do termo de “domínio de dependência”.

Considerando, por exemplo, a equação da advecção (citar a equação da advecção) em uma aproximação a partir de métodos de diferença em sua forma explícita $\frac{U_i^{n+1}-U_i^n}{\mathrm{\Delta }t}+v\frac{U_i^n-U_{i-1}^n}{\mathrm{\Delta }x}=0$, diz-se que o valor de $U_i^{n+1}$ depende dos valores anteriores de $U_i^n$ e $U_{i-1}^n$, e esses dois pontos dependem, novamente, de outros dois pontos em tempos anteriores. Todos esses pontos dependentes formam o domínio de dependência do valor $U_i^{n+1}$, representado abaixo (botar figura do domínio de dependência).

A condição CFL enuncia que a condição mínima para que haja estabilidade em métodos de diferenças o domínio de dependência da equação diferencia parcial, dado por sua equação característica, deve estar situado dentro do domínio de dependência do esquema numérico.

A partir dessa condição, define-se o número CFL como

$r\equiv v\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}$

onde $\mathrm{\Delta }t$ é um passo na malha temporal e $\mathrm{\Delta }x$ é um passo na malha espacial. Para casos unidimensionais, como o que será tratado aqui, a condição CFL é satisfeita se r

$\left|r\right|\le 1$

Um uso da condição CFL é determinar o tamanho do passo temporal, sabendo-se $v$ e $\mathrm{\Delta }x$:

$\mathrm{\Delta }t=c_{cfl}\frac{\mathrm{\Delta }x}{\left|v\right|}$

onde o fator $c_{cfl}<1$.

O número CFL pode ser entendido como um fluxo numérico advectivo adimensionalisado pelas malhas espacial e temporal do problema. De um ponto de vista matemático, ele garante que o domínio numérico de dependência é sempre maior que o domínio físico. De um ponto de vista físico, é garantido que a velocidade de propagação de qualquer perturbação, como uma onda, seja menor ou ao menos igual que a velocidade de propagação numérica, fazendo com que a distância propagada pela perturbação não seja maior do que a divisão da malha espacial:

$\left|v\right|\mathrm{\Delta }t\le \mathrm{\Delta }x$

O Problema Físico

O Modelo de Corda Ideal

Para uma primeira abordagem da equação da onda, podemos primeiro dividir o comprimento $L$ da corda em $K$ intervalos de comprimentos iguais, dessa forma $\mathrm{\Delta }x=\frac{L}{K}$. Cada intervalo é discretizado, representado por $x_i$, $i=0,1,...,K$. Também podemos dividir o tempo em intervalos iguais $\mathrm{\Delta }t$ e denotá-los como $t_n$, $n=0,1...$.

Tendo feita a discretização das variáveis, podemos aproximar a equação da onda por diferenciação finita, utilizando derivadas centradas da seguinte forma:

$\frac{{\partial }^{2}u}{\partial t^2}=\frac{U_i^{n+1}-2U_i^n+U_i^{n-1}}{\mathrm{\Delta }{t}^2}$

$\frac{{\partial }^{2}u}{\partial x^2}=\frac{U_{i+1}^n-2U_i^n+U_{i-1}^n}{\mathrm{\Delta }{x}^2}$

onde $U_i^n$ representa o valor discretizado de $u(x_i,t_n)$.

Assim, chegamos em uma equação discretizada:

$\frac{U_i^{n+1}-2U_i^n+U_i^{n-1}}{(\mathrm{\Delta }t{)}^2}=v^2\frac{U_{i+1}^n-2U_i^n+U_{i-1}^n}{(\mathrm{\Delta }x{)}^2}$.

Sabendo que essa discretização da equação da onda pode ser verificada como sendo o método Leapfrog (ver seção do método de Leapfrog), podemos resolver a equação para $U_i^{n+1}$ para sabermos o deslocamento de uma partição da corda no momento de tempo seguinte, assim obtendo

$U_i^{n+1}=2(1-r^2)U_i^n+r^2[U_{i+1}^n+U_{i-1}^n]-U_i^{n-1}$,

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

Um Quadro Mais Realístico - O Modelo de Corda Rígida

Para nos aproximarmos de um modelo mais real, podemos adicionar um termo à equação original da onda que corresponde ao efeito de fricção em uma corda. De acordo com , a equação da onda mais geral com efeito de fricção pode ser escrita como

$\frac{{\partial }^{2}u}{\partial t^2}=v^2(\frac{{\partial }^{2}u}{\partial x^2}-ϵL^2\frac{{\partial }^{4}u}{\partial x^4}),$

onde $v$ é a velocidade transversal de propagação do pulso na corda, dada pela relação $v=\sqrt[]{\frac{T}{\rho }}$ (sendo $T$ a tensão na corda e $\rho$ a densidade linear da mesma), $ϵ$ é um parâmetro adimensional de fricção que representa a rigidez da corda e $L$ o comprimento da corda.

O parâmetro $ϵ$ é dado por

$$,

onde $\kappa$ é o raio da corda, $E$ é o Módulo de Young e $S$ a área da secção da corda.

Discretizamos a equação da seguinte maneira:

$\frac{{\partial }^{2}u}{\partial t^2}=\frac{U_i^{n+1}-2U_i^n+U_i^{n-1}}{\mathrm{\Delta }{t}^2}$

$\frac{{\partial }^{2}u}{\partial x^2}=\frac{U_{i+1}^n-2U_i^n+U_{i-1}^n}{\mathrm{\Delta }{x}^2}$

$\frac{{\partial }^{4}u}{\partial x^4}=\frac{U_{i+2}^n-4U_{i+1}^n+6U_i^n-4Y_{i-1}^n+U_{i-2}^n}{\mathrm{\Delta }{x}^4}$

e resolvemos para $U_i^{n+1}$, obtendo:

$U_i^{n+1}=[2-2r^2-6ϵr^2{K}^2]U_i^n-U_i^{n-1}+r^2[1+4ϵK^2][U_{i+1}^n+U_{i-1}^n]-ϵr^2{K}^2[U_{i+2}^n+U_{i-2}^n].$

O fato de essa discretização depender do deslocamento da corda em posições $i-2$ e $i+2$ implica em precisarmos simular "pontos fantasmas" quando integramos os extremos das cordas. Para fazermos isso, podemos ou utilizar a aproximação $u_{-1}^n=-u_{+1}^n$ ou podemos considerar esses "pontos fantasmas" como pontos presos e, portanto, sempre iguais a zero.

Os Métodos Utilizados

Foi realizada uma abordagem ao problema da corda real a partir de três métodos diferentes de integração numérica. Os três são métodos para fins de resolução de equações diferenciais parciais da forma apresentada anteriormente.

O método mais básico é chamado de FTCS (Forward-Time-Centered-Space) e consiste em duas expansões de Taylor ao redor do ponto $x_j$:

$u(x_i+\mathrm{\Delta }x,t^n)=u(x_i,t^n)+\frac{\partial u}{\partial x}(x_i,t^n)\mathrm{\Delta }x+\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}(x_i,t^n)\mathrm{\Delta }{x}^2+{𝒪}(\mathrm{\Delta }{x}^3),$

$u(x_i-\mathrm{\Delta }x,t^n)=u(x_i,t^n)-\frac{\partial u}{\partial x}(x_i,t^n)\mathrm{\Delta }x+\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}(x_i,t^n)\mathrm{\Delta }{x}^2+{𝒪}(\mathrm{\Delta }{x}^3).$

Subtraindo as duas expressões, encontramos a expressão

${\frac{\partial u}{\partial x}|}_i^n=\frac{u_{i+1}^n-u_{i-1}^n}{2\mathrm{\Delta }x}+{𝒪}(\mathrm{\Delta }{x}^2)$,

A qual podemos substituir na equação da onda, juntamente com a discretização da derivada parcial temporal. Temos então que, para um sistema linear de equações hiperbólicas:

${\text{U}}_i^{n+1}={\text{U}}_i^n-\frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}[{\text{F}}_{i+1}^n-{\text{F}}_{i-1}^n]+{𝒪}(\mathrm{\Delta }{t}^2,\mathrm{\Delta }{x}^2\mathrm{\Delta }t)$

Visto que essa última notação é mais genérica, ela será utilizada para a explicação dos métodos posteriores.

O Método de Lax-Friedrichs

O método de Lax-Friedrichs consiste em substituir o termo ${\text{U}}_i^n$ com sua respectiva média espacial, i.e., $u_i^n=(u_{i+1}^n+u_{i-1}^n)/2$. Logo, temos a seguinte equação de recorrência:

${\text{U}}_i^{n+1}=\frac{1}{2}({\text{U}}_{i+1}^n+{\text{U}}_{i-1}^n)-\frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}[{\text{F}}_{i+1}^n-{\text{F}}_{i-1}^n]+{𝒪}(\mathrm{\Delta }{x}^2)$

O Método de Leapfrog

Tanto o método FTCS quanto o método de Lax-Friedrichs são métodos de primeira ordem para a derivada temporal. Nessas circunstâncias, $v\mathrm{\Delta }t$ deve ser significantemente menor do que $\mathrm{\Delta }x$, muito abaixo do limite imposto pela condição de Courant.

Uma nova expressão para a derivada temporal, com precisão de segunda ordem é dada por

${\frac{\partial u}{\partial t}|}_i^n=\frac{u_i^{n+1}-u_i^{n-1}}{2\mathrm{\Delta }t}+{𝒪}(\mathrm{\Delta }{t}^2).$

Substituindo a nova expressão acima no método de FTCS discutido anteriormente, encontramos o método de Leapfrog:

${\text{U}}_i^{n+1}={\text{U}}_i^{n-1}-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}[{\text{F}}_{i+1}^n-{\text{F}}_{i-1}^n]+{𝒪}(\mathrm{\Delta }{x}^2).$

Como o método de Leapfrog foi o mais aplicado na resolução do problema em questão, é interessante um aprofundamento maior do método. Podemos adaptar o método de Leapfrog para o sistema de equações definido para a equação da onda ao fazermos

$k_{i+\frac{1}{2}}^n\equiv v{\frac{\partial u}{\partial x}|}_{i+\frac{1}{2}}^n=v\frac{u_{i+1}^n-u_i^n}{\mathrm{\Delta }x}+{𝒪}(\mathrm{\Delta }x)(1)$

$s_i^{n+\frac{1}{2}}\equiv {\frac{\partial u}{\partial t}|}_i^{n+\frac{1}{2}}=\frac{u_i^{n+1}-u_i^n}{\mathrm{\Delta }t}+{𝒪}(\mathrm{\Delta }t)(2)$

Com a representação Leapfrog das equações do sistema de três equações, temos:

$k_{i+\frac{1}{2}}^{n+1}=k_{i+\frac{1}{2}}^n+r(s_{i+1}^{n+\frac{1}{2}}-s_i^{n+\frac{1}{2}})+{𝒪}(\mathrm{\Delta }{x}^2)(3)$

$s_i^{n+\frac{1}{2}}=s_i^{n-\frac{1}{2}}+r(k_{i+\frac{1}{2}}^n-k_{i-\frac{1}{2}}^n)+{𝒪}(\mathrm{\Delta }{x}^2)(4)$

Com essas duas equações, podemos fazer uma integração utilizando o método de Euler para obter $u_i^{n+1}$, ou seja, o deslocamento de um determinado ponto no próximo instante de tempo:

$u_i^{n+1}=u_i^n+\frac{\mathrm{\Delta }t}{2}{s}_i^{n+\frac{1}{2}}+{𝒪}(\mathrm{\Delta }{x}^2).$

Contudo, podemos fazer uma simples substituição das equações $(1)$ e $(2)$ nas equações $(3)$ e $(4)$ e, assim, obtemos que a representação de Leapfrog da equação da onda é dada pela discretização de segunda ordem da própria equação da onda, com ${𝒪}(\mathrm{\Delta }{t}^2,\mathrm{\Delta }{x}^2)$. Isso nos dá uma solução de "um passo", onde só precisamos efetuar o cálculo da equação discretizada.

O Método de Lax-Wendroff

O método de Lax-Wendroff é a extensão do método de Lax-Friedrichs de segunda ordem. Calculamos o vetor $U$ a partir de um passo médio de Lax-Friedrichs:

${\text{U}}_{i+\frac{1}{2}}^{n+\frac{1}{2}}=\frac{1}{2}[{\text{U}}_{i+1}^n+{\text{U}}_i^n]-\frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}[{\text{F}}_{i+1}^n-{\text{F}}_i^n]+{𝒪}(\mathrm{\Delta }{x}^2)$,

${\text{U}}_{i-\frac{1}{2}}^{n+\frac{1}{2}}=\frac{1}{2}[{\text{U}}_i^n+{\text{U}}_{i-1}^n]-\frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}[{\text{F}}_i^n-{\text{F}}_{i-1}^n]+{𝒪}(\mathrm{\Delta }{x}^2)$,

E encontramos os fluxos ${\text{F}}_{i\pm \frac{1}{2}}^{n+\frac{1}{2}}$ a partir dos valores de ${\text{U}}_{i\pm \frac{1}{2}}^{n+\frac{1}{2}}.$

Logo, com um meio passo de Leapfrog, temos a expressão final do método:

${\text{U}}_i^{n+1}={\text{U}}_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}[{\text{F}}_{i+\frac{1}{2}}^{n+\frac{1}{2}}-{\text{F}}_{i-\frac{1}{2}}^{n+\frac{1}{2}}]+{𝒪}(\mathrm{\Delta }{x}^2)$

Análise e Discussão dos Resultados

Escolhemos para simular quatro diferentes cordas: as cordas C2, C4 e C7 de um piano cujos dados foram obtidos no artigo sobre cordas do Chaigne [3] e uma corda com dados pré-estabelecidos encontrados no livro de física computacional do Giordano [1]. Os dados das cordas C2, C4 e C7 estão na tabela abaixo.

	C2	C4	C7
Comprimento ($m$)	$1.90$	$0.62$	$0.09$
Massa($g$)	$35$	$3.93$	$0.467$
Tensão ($N$)	$750$	$670$	$750$
Divisões	$100$	$50$	$16$
Amostragem de sinal	$16\times 10^3$	$32\times 10^3$	$96\times 10^3$
Parâmetro de fricção	$7.5\times 10^{-6}$	$3.82\times 10^{-5}$	$8.67\times 10^{-4}$

Desses dados, temos que a densidade linear de massa das cordas é dada por

$\rho =\frac{M}{l},$

onde $M$ é a massa e $l$ o comprimento da corda,

$v=\sqrt{\frac{T}{\rho }},$

onde $T$ é a tensão na corda,

$\mathrm{\Delta }x=\frac{l}{K},$

$\mathrm{\Delta }t=\frac{1}{f_e},$

e

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

Para a corda com dados pré-estabelecidos, utilizamos uma corda com 2 metros de comprimento, com velocidade de propagação da onda sendo 300$m/s$, $\mathrm{\Delta }x$ de 0.01$m$, $\mathrm{\Delta }t=\frac{\mathrm{\Delta }x}{4v}$ e parâmetro de fricção de $7.5\times 10^{-6}$. Supondo a corda inicialmente em repouso, temos que em $t=0$ a corda recebe em seu centro o equivalente à batida do martelo do piano. Supomos que esse estímulo possuía o formato aproximado de uma Gaussiana com amplitude de $\frac{1}{4}$ do comprimento da corda. Então, com o estado inicial sendo um pulso com o formato de um pacote gaussianico e os dados da tabela, simulamos a propagação de ondas em cada uma das cordas.

Utilizando o método de Leapfrog, foi realizada uma primeira simulação para uma onda ideal em uma corda:

Corda C2:

Corda C4:

Corda C7:

Corda com as definições do Giordano:

Análise de Erro e Estabilidade dos Métodos

Estabilidade do método Leapfrog

Pela estabilidade de Von Neumann, podemos escrever que

$r\le 1.$

Para $r=1$, a equação da discretização da onda pode ser reescrita como

$U_i^{n+1}=U_{i+1}^n+U_{i-1}^n-U_i^{n-1}.$

Essa escolha com

$\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}=v$

nos dá a solução exata sem dispersão numérica. Contudo, $r=1$ é válido somente no caso de uma corda ideal. É conveniente escrever a condição acima em termos da amostragem de sinal $f_e=\frac{1}{\mathrm{\Delta }t}$ e a frequência fundamental da corda $f_1=\frac{v}{2L}$, o que nos leva a

$K{f}_1\le \frac{f_e}{2}$

O teorema de Nyquist diz que a frequência superior no espectro deve ser menor do que $\frac{f_e}{2}$ para evitar serrilhamento e para garantir uma reconstrução única e contínua. Logo, no caso ideal quando as autofrequências da corda são igualmente espaçadas ($\mathrm{\Delta }f=f_1$), a condição de Nyquist indica que o número máximo de frequências no espectro é $K$. Isso significa que essa condição pode ser usada para selecionar o número apropriado de pontos espaciais $K$ para a corda. Entretanto, como $K$ é um inteiro, apenas valores discretos da frequência natural podem ser obtidos sem erros de trucamento, ou seja, usando $r=1$. Como series discretas não costumam ser utilizadas, precisamos aceitar pequenos erros de truncamento para ajustar $f_1$, ou seja, utilizando $r<1$. No caso de uma onda com fricção, temos que $r=\frac{1}{4}$ é um valor de boa estabilidade.

---