Pêndulo Simples

Considerando o problema do pêndulo simples mostrado na Figura.

Sendo a matriz de transformação é dada por:

\begin{equation} \textbf{T}_{\theta} = \begin{bmatrix} c\theta & s\theta & 0 \\ -s\theta & c\theta & 0 \\ 0 & 0 & 1 \end{bmatrix} \end{equation}

O vetor posição descrito na base $\textit{S}_{1}$ é dada por:

\begin{equation} _{S_1}^{S_1}\textbf{p}^P = \begin{Bmatrix} 0 \\ L \\ 0 \end{Bmatrix} \end{equation}

Enquanto a aceleração do pêndulo pode ser escrita conforme segue.

\begin{equation} _{S_1}^{I}\textbf{a}^P = \; _{S_1}^{I}\textbf{a}^O + \; _{S_1}\textbf{a}^{P/S_1} + \; _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \big( _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \: _{S_1}\textbf{p}^{P/S_1} \big) + 2 \; _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \; _{S_1}\textbf{v}^{P/S_1} + _{S_1}^{I}\boldsymbol{\alpha}^{S_1} \times \; _{S_1}\textbf{p}^{P/S_1} \end{equation}

A partir de uma análise do movimento do sistema observa-se que $\textbf{a}^{O} = \textbf{a}^{P/S_1} = 2 \; ^{I}\boldsymbol{\omega}^{S_1} \times \textbf{v}^{P/S_1} = 0$ e com isso, tem-se:

\begin{equation} _{S_1}^{I}\boldsymbol{\alpha}^{S_1} \times \; _{S_1}\textbf{p}^{P/S_1} = \begin{bmatrix} \textbf{e}_1^{(1)} & \textbf{e}_2^{(1)} & \textbf{e}_3^{(1)} \\ 0 & 0 & \ddot{\theta} \\ 0 & L & 0 \\ \end{bmatrix} = \begin{Bmatrix} 0 \\ -L\ddot{\theta} \\ 0 \end{Bmatrix} \end{equation}\begin{equation} _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \; \big( _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \, _{S_1}\textbf{p}^{P/S_1} \big) = \; _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \; \begin{bmatrix} \textbf{e}_1^{(1)} & \textbf{e}_2^{(1)} & \textbf{e}_3^{(1)} \\ 0 & 0 & \dot{\theta} \\ 0 & L & 0 \\ \end{bmatrix} = \begin{bmatrix} \textbf{e}_1^{(1)} & \textbf{e}_2^{(1)} & \textbf{e}_3^{(1)} \\ 0 & 0 & \dot{\theta} \\ -L\dot{\theta} & 0 & 0 \\ \end{bmatrix} = \begin{Bmatrix} 0 \\ -L\dot{\theta}^2 \\ 0 \end{Bmatrix} \end{equation}

Portanto, a aceleração é dada por:

\begin{equation} _{S_1}^{I}\textbf{a}^P = \begin{Bmatrix} -L\ddot{\theta} \\ -L\dot{\theta}^2 \\ 0 \end{Bmatrix} \end{equation}

Avaliando as forças no pêndulo, considera-se uma força resultante $F$ composta pelo peso, $P$, uma força dissipativa viscosa linear, $F_D$, uma excitação externa, $F_E$ e uma força na haste, $F_H$, conforme mostrado na Figura anterior.

Com isso, tem-se:

\begin{equation} _{S_1}\textbf{P} = \textbf{T}_{\theta}\; _{I}\textbf{P} = \begin{bmatrix} c\theta & s\theta & 0 \\ -s\theta & c\theta & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{Bmatrix} 0 \\ mg \\ 0 \end{Bmatrix} = \begin{Bmatrix} mgs\theta \\ mgc\theta \\ 0 \end{Bmatrix} \end{equation}\begin{equation} _{S_1}\textbf{F}_D = \begin{Bmatrix} CL\dot{\theta} \\ 0 \\ 0 \end{Bmatrix} \end{equation}\begin{equation} _{S_1}\textbf{F}_E = \begin{Bmatrix} -F(t) \\ 0 \\ 0 \end{Bmatrix} \end{equation}\begin{equation} _{S_1}\textbf{F}_H = \begin{Bmatrix} 0 \\ -T \\ 0 \end{Bmatrix} \end{equation}

Desta forma, aplicando a segunda lei de Newton, $m \;_{S_1}^{I}\textbf{a}^P = \,_{S_1}\textbf{F}$, tem-se que:

\begin{equation} m \; _{S_1}\textbf{a}^P = m \begin{Bmatrix} -L\ddot{\theta} \\ -L\dot{\theta}^2 \\ 0 \end{Bmatrix} = \begin{Bmatrix} mgs\theta + CL\dot{\theta} - F(t) \\ mgc\theta - T \\ 0 \end{Bmatrix} \end{equation}

O que resulta nas seguintes equações:

\begin{equation} mL\ddot{\theta} + CL\dot{\theta} + mgs\theta = F(t) \\ T = mL\dot{\theta}^2 + mgc\theta \end{equation}

Note que a segunda esquação estabelece o valor da tensão na haste enquanto a primeira equação define o movimento do pêndulo que, trata-se de uma equação diferencial não-linear. Considerando pequenos ângulos, $\theta ≪ 1$, pode-se assumir que $s\theta \approx \theta$. Essa simplificação usual transforma a equação de movimento em uma equação diferencial linear, que representa um oscilador linear:

\begin{equation} mL\ddot{\theta} + CL\dot{\theta} + mg\theta = F(t) \end{equation}

A fim de obter a solução numérica da equação de movimento utilizando o método de Runge-Kutta de 4ª ordem. Fazemos a seguinte mudança de variáveis: $x_1 = \theta$, $x_2 = \dot{\theta}$, permitindo que a equação de movimento não-linear e linear sejam escritas na forma canônica, conforme expresso a seguir.

$\textbf{Não-linear}$

\begin{equation} \begin{cases} \dot{x}_1 = x_2 \\ \dot{x}_2 = \frac{F(t)}{mL} - \frac{C}{m}x_2 - \frac{g}{L}s(x_1) \end{cases} \end{equation}

$\textbf{Linear}$

\begin{equation} \begin{cases} \dot{x}_1 = x_2 \\ \dot{x}_2 = \frac{F(t)}{mL} - \frac{C}{m}x_2 - \frac{g}{L}x_1 \end{cases} \end{equation}

onde $\textbf{x}$ é o vetor de variáveis de estado dado por,

\begin{equation} \textbf{x} = \begin{Bmatrix} x_1 \\ x_2 \\ \end{Bmatrix} = \begin{Bmatrix} \theta \\ \dot{\theta} \\ \end{Bmatrix} \end{equation}

Em seguida vamos definir uma função em Python para as equações de movimento não-linear e linear.

In [20]:
%matplotlib inline
from numpy import sin, cos
import numpy as np
import matplotlib.pyplot as plt
from RK4 import integrate
#============================================================================#
# Right-hand side of simple pendulum differential equations                  #
#============================================================================#
def lin_sim_pen(x, t):

    ff    = np.zeros_like(x)
    ff[0] = x[1]
    ff[1] = amp*sin(omega*t) - ksi*x[1] -(G/L)*x[0]

    return ff

def nonlin_sim_pen(x, t):

    ff    = np.zeros_like(x)
    ff[0] = x[1]
    ff[1] = amp*sin(omega*t) - ksi*x[1] -(G/L)*sin(x[0])

    return ff

Vamos criar uma função para plotar os resultados obtidos.

In [49]:
def plots(y, X1, X2, label, ls, ax1, ax2):
        
    t     = y[:,0]                      # time array
    theta = y[:,1]                      # angle array
    w     = y[:,2]                      # angular velocity array
    theta = np.degrees(theta)           # converts from rad to degree
    w     = np.degrees(w)               # converts from rad to degree
    #========================================================================#
    # Plot results                                                           #
    #========================================================================#
    ax1.plot(t, theta, ls, label=label, alpha=1.0)
    ax1.set_xlabel(' Tempo (s) ')
    ax1.set_ylabel('$\\theta [\degree]$')
    ax1.set_title(' Pêndulo Simples')
    ax1.set_xlim(t[0], t[-1])
    ax1.set_ylim(-1.5*np.max(theta), 1.5*np.max(theta))
    ax1.legend()
    #
    # State space
    #
    ax2.plot(theta, w, ls, label=label)
    ax2.set_xlabel('$\\theta [\degree]$')
    ax2.set_ylabel('$\\omega [\degree/s]$')
    ax2.set_title(' Pêndulo Simples')
    ax2.set_xlim(1.5*np.min(theta), 1.5*np.max(theta))
    ax2.legend()


Agora, podemos definir os parâmtros do problema e utlizar as funções de integração numérica pelo método de Runge-Kutta, implementadas anteriormente, para obter a resposta do sistema. Por fim plotamos os resultados.

Caso 1: Vibração livre sem amortecimento

Primeiramente, assumimos como condições iniciais, uma condição de pequeno deslocamento, $\theta_0 = -10°$ e velocidade angular de $\dot{\theta}_0 = 0°$/s.

In [50]:
#========================================================================#
# Inputs                                                                 #
#========================================================================#
G  = 9.8        # acceleration due to gravity, in m/s^2
L  = 1.0        # length of pendulum in m
ksi = 0.0       # viscous damping ratio
amp = 0.0       # excitation amplitude
omega = 0.0     # excitation frequency
t0 = 0.         # initial time
dt = 0.05       # time step
N  = 1000       # number of time steps

th1 = -10.0     # initial angle (degrees)
w1  = 0.0       # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1])  # converts to radians


plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Linear Simple Pendulum Solution using RK4                              #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4                           #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
Your browser does not support this tag.

Conforme pode ser observado, para pequenos ângulos, a resposta do sistema são semelhantes para a equação de movimento linear e não-linear.

Caso 2: Vibração livre sem amortecimento

Por outro lado, para ângulos maiores, digamos $\theta_{0} = -120°$ e mantendo $\dot{\theta}_0 = 0°/s$, o pêndulo simples é não-linear e é possível observar a influência na resposta do sistema, conforme pode ser visualizado nos resultados que seguem.

In [51]:
#========================================================================#
# Inputs                                                                 #
#========================================================================#
G  = 9.8        # acceleration due to gravity, in m/s^2
L  = 1.0        # length of pendulum in m
ksi = 0.0       # viscous damping ratio
amp = 0.0       # excitation amplitude
omega = 0.0     # excitation frequency
t0 = 0.         # initial time
dt = 0.05       # time step
N  = 1000       # number of time steps

th1 = -120.0    # initial angle (degrees)
w1  = 0.0       # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1])  # converts to radians


plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Linear Simple Pendulum Solution using RK4                              #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4                           #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
Your browser does not support this tag.

Caso 3: Vibração livre sem amortecimento

Considerando neste momento as seguintes condições iniciais, $\theta_{0} = -120°$ e $\dot{\theta}_0 = 200 °/s$, onde uma velocidade angular não nula é assumida, tem-se o comparativo mostrado a seguir.

In [61]:
#========================================================================#
# Inputs                                                                 #
#========================================================================#
G  = 9.8        # acceleration due to gravity, in m/s^2
L  = 1.0        # length of pendulum in m
ksi = 0.0       # viscous damping ratio
amp = 0.0       # excitation amplitude
omega = 0.0     # excitation frequency
t0 = 0.         # initial time
dt = 0.05       # time step
N  = 1000       # number of time steps

th1 = -120.0    # initial angle (degrees)
w1  =  200.0    # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1])  # converts to radians


plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4                           #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
#========================================================================#
# Linear Simple Pendulum Solution using RK4                              #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
ax1.set_ylim(-200.0, 400);
ax2.set_xlim(-200.0, 1000);
Your browser does not support this tag.

Caso 4: Vibração livre sem amortecimento

Em seguida, considera-se uma velocidade angular iniciais ainda maior, $\theta_0=500°/𝑠$, mantendo o mesmo valor de deslocamento angular inicial. Observa-se que o pêndulo linear tem um comportamento peculiar, que $\textbf{não}$ tem nenhum sentido físico para um sistema sem nenhum tipo de forçamento externo. Isso ocorre porque é assumido a hipótese de pequenas rotações que não condiz com a situação analisada. Em seguida, note o comportamento esperado para o pêndulo quando considera-se a não-linearidade.

In [69]:
#========================================================================#
# Inputs                                                                 #
#========================================================================#
G  = 9.8        # acceleration due to gravity, in m/s^2
L  = 1.0        # length of pendulum in m
ksi = 0.0       # viscous damping ratio
amp = 0.0       # excitation amplitude
omega = 0.0     # excitation frequency
t0 = 0.         # initial time
dt = 0.05       # time step
N  = 1000       # number of time steps

th1 = -120.0    # initial angle (degrees)
w1  =  500.0    # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1])  # converts to radians


plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4                           #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
#========================================================================#
# Linear Simple Pendulum Solution using RK4                              #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
ax1.set_ylim(-300.0, 600);
ax2.set_xlim(-300.0, 1000);
Your browser does not support this tag.
Your browser does not support this tag.

Caso 5: Vibração livre com amortecimento

Considerando um amortecimento viscoso, $\xi = 0.08$, assumimos como condições iniciais, uma condição de pequeno deslocamento, $\theta_0 = -10°$ e velocidade angular de $\dot{\theta}_0 = 0°$/s. Note que a resposta do sistema são semelhantes e apresentam decaimento exponencial com o decorrer do tempo.

In [72]:
#========================================================================#
# Inputs                                                                 #
#========================================================================#
G  = 9.8        # acceleration due to gravity, in m/s^2
L  = 1.0        # length of pendulum in m
ksi = 0.08      # viscous damping ratio
amp = 0.0       # excitation amplitude
omega = 0.0     # excitation frequency
t0 = 0.         # initial time
dt = 0.05       # time step
N  = 1000       # number of time steps

th1 = -10.0    # initial angle (degrees)
w1  =   0.0    # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1])  # converts to radians


plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4                           #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
#========================================================================#
# Linear Simple Pendulum Solution using RK4                              #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
Your browser does not support this tag.

Caso 6: Vibração livre com amortecimento

Considerando um amortecimento viscoso, $\xi = 0.08$ e assumindo um ângulo maior, $\theta_0 = -120°$ e velocidade angular de $\dot{\theta}_0 = 0°$/s. O pêndulo simples é não-linear e é possível observar a influência na resposta do sistema, conforme pode ser visualizado nos resultados que seguem.

In [73]:
#========================================================================#
# Inputs                                                                 #
#========================================================================#
G  = 9.8        # acceleration due to gravity, in m/s^2
L  = 1.0        # length of pendulum in m
ksi = 0.08      # viscous damping ratio
amp = 0.0       # excitation amplitude
omega = 0.0     # excitation frequency
t0 = 0.         # initial time
dt = 0.05       # time step
N  = 1000       # number of time steps

th1 = -120.0   # initial angle (degrees)
w1  =   0.0    # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1])  # converts to radians


plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4                           #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
#========================================================================#
# Linear Simple Pendulum Solution using RK4                              #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
Your browser does not support this tag.

Caso 7: Vibração livre com amortecimento

Considerando um amortecimento viscoso, $\xi = 0.2$ e assumindo um ângulo maior, $\theta_0 = -120°$ e velocidade angular de $\dot{\theta}_0 = 500°$/s. O pêndulo simples é não-linear e é possível observar a influência na resposta do sistema, conforme pode ser visualizado nos resultados que seguem.

In [83]:
#========================================================================#
# Inputs                                                                 #
#========================================================================#
G  = 9.8        # acceleration due to gravity, in m/s^2
L  = 1.0        # length of pendulum in m
ksi = 0.2       # viscous damping ratio
amp = 0.0       # excitation amplitude
omega = 0.0     # excitation frequency
t0 = 0.         # initial time
dt = 0.05       # time step
N  = 1000       # number of time steps

th1 = -120.0    # initial angle (degrees)
w1  =  500.0    # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1])  # converts to radians


plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4                           #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
#========================================================================#
# Linear Simple Pendulum Solution using RK4                              #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen)  # matrix with results
t     = y[:,0]                        # time
theta = y[:,1]                        # angle
w     = y[:,2]                        # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
ax1.set_ylim(-300.0, 1800);
ax2.set_xlim(-200.0, 1800);
Your browser does not support this tag.