Física Computacional - Pêndulo Amortecido

Dando prosseguimento ao estudo do movimento oscilatório iniciado neste post, agora estarei considerando um fator de amortecimento na descrição do movimento do pêndulo.

A equação diferencial do pêndulo amortecido é a seguinte:

$$\frac{d^2\theta}{dt^2} + q\frac{d\theta}{dt} + \frac{g}{l}\theta = 0$$

Onde $q$ é o fator de amortecimento.

Esta equação é uma EDO homogênea de segunda ordem e possui solução analítica do tipo $\theta(t)=Ce^{\lambda t}$. Substituindo essa solução na EDO, temos:

$\lambda^2 + q\lambda + \frac{g}{l} = 0$

$\lambda = \frac{-q \pm \sqrt{q^2 - 4\frac{g}{l}}}{2}$

$\lambda = \frac{-q}{2} \pm \sqrt{\frac{q^2}{4} - \frac{g}{l}}$

Dessa forma, temos três casos possíveis de amortecimento:

Primeiro: $\frac{q^2}{4} - \frac{g}{l} > 0$

Este caso é conhecido como amortecimento supercrítico. Nele, o termo dentro da raiz quadrada é positivo, assim a solução da EDO fica sendo:

$\theta(t) = C_1e^{\lambda_{1}t} + C_{2}e^{\lambda_{2}t}$

Onde,

$\lambda_{1} = \frac{-q}{2} + \sqrt{\frac{q^2}{4}-\frac{g}{l}}$
$\lambda_{2} = \frac{-q}{2} - \sqrt{\frac{q^2}{4}-\frac{g}{l}}$
$C_{1}$ e $C_{2}$, são constantes á serem determinadas usando as condições iniciais.


Segundo: $\frac{q^2}{4} - \frac{g}{l} = 0$

Este caso é conhecido como amortecimento crítico. Nele, o termo dentro da raiz quadrada é nulo, assim a solução da EDO fica sendo:

$\theta(t)=C_{1}e^{\lambda t}+C_{2}te^{\lambda t}$

Onde,

$\lambda = -\frac{q}{2}$
$C_{1}$ e $C_{2}$, são constantes á serem determinadas usando as condições iniciais.

Terceiro: $\frac{q^2}{4} - \frac{g}{l} < 0$

Este caso é conhecido como amortecimento subcrítico. Nele, o termo dentro da raiz quadrada é negativo, assim a solução da EDO fica sendo:

$\theta(t)=e^{\alpha t}(C_{1}cos(\beta t) + C_{2}sin(\beta t))$

Onde,

$\lambda=\alpha + \beta i$
$\alpha = -\frac{q}{2}$
$\beta = \sqrt{\frac{g}{l} - \frac{q^2}{4}}$

Solução Numérica: Método de Euler-Cromer

Conforme foi discutido do post anterior, foi escolhido o método de Euler-Cromer para resolver a EDO numericamente.

Para aplicar o método foi escrito um programa em linguagem C. Segue abaixo o código do programa:

/*
Programa: euler_cromer.c
Objetivo: Usar o método de Euler-Cromer para resolver a equação
                  diferencial do Pêndulo Simples, amortecido, mas
                  não-forçado e para pequenas oscilações
Autor: Alex Souza Marques
*/

#include <stdio.h>
#include <math.h>

#define g   9.8     // Aceleração da gravidade
#define l   1         // Comprimento do pêndulo
#define dt  0.04  // Tamanho do passo
#define m   1       // Massa do pêndulo
#define q   0.       // Constante de amortecimento

int main(){
    FILE *out;
    float t=0;
    float w=0, o=0.2;   // Omega e Theta
    float K, U;         // Energia Cinética e Energia Potencial

    out=fopen("euler_cromer.dat","w");   // Arquivo para saída dos dados
    K=0.5*m*l*l*w*w;
    U=m*g*l*(1-cos(o));
    fprintf(out,"%f\t%f\t%f\t%f\t%f\t%f\n",t,o,w,K,U,K+U);  // Imprime as condicoes iniciais

    /*Aplicando o método de Euler-Cromer*/
    for(t=dt;t<=5*2*M_PI*sqrt(l/g);t+=dt){  // Itera por 5 períodos do pêndulo
        w=w-(g/l)*o*dt -q*w*dt;
        o=o+w*dt;
        K=0.5*m*l*l*w*w;
        U=m*g*l*(1-cos(o));
        fprintf(out,"%f\t%f\t%f\t%f\t%f\t%f\n",t,o,w,K,U,K+U);  // Imprime o resultado
    }

    fclose(out);

    return 0;
}

Esse programa gera o arquivo "euler_cromer.dat" com a saída dos resultados. Foi usado o software gnuplot para ilustrar os resultados. Abaixo segue um exemplo de script, que foi criado para gerar os gráficos no gnuplot.

# Script: pendulo.gpl
# Objetivo: Imprimir os dados do programa euler_cromer.c
# Autor: Alex Souza Marques

reset;
set grid;
set xlabel "Tempo";
set ylabel "Energia";
set title "Amortecimento Super-Crítico - Energia";
#set key left;
#set yrange [0:0.27];
plot "euler_cromer_12_3.dat" u 1:4 title "Energia Cinética" w l;
rep "euler_cromer_12_3.dat" u 1:5 title "Energia Potêncial" w l;
rep "euler_cromer_12_3.dat" u 1:6 title "Energia Total" w l;
#rep "euler_cromer_9_3.dat" u 2:3 title "q=9.3" w l;
#rep "euler_cromer_12_3.dat" u 2:3 title "q=12.3" w l;

Usando a solução analítica e os valores de $l$ e $g$ definidos anteriormente, foi determinado o valor de $q$ no regime de transição (crítico) como sendo aproximadamente $6,3$. Abaixo, seguem os resultados para vários valores possíveis de $q$.








Foi usado o Gnuplot para ilustrar também a perda de energia do pêndulo em relação ao tempo nos três regimes de amortecimento possíveis. 









Conclusão

Os resultados foram compatíveis com o que era esperado pelo resultado analítico, ilustrando o quanto o método de Euler-Cromer é preciso e eficaz para esse tipo de problema.

No próximo post estarei falando do caso não-linear. Considerando uma força externa e fora do limite de pequenas oscilações, será falado sobre o limiar entre um movimento regular e predizível e o caos. Até lá.



Comentários