28.40 A dinâmica de um sistema massa-mola amortecido forçado pode ser representada pela seguinte EDO de segunda ordem:
onde
. Use um método numérico para determinar o deslocamento
e a velocidade (
) como uma função do tempo com as condições iniciais
. Expresse seus resultados graficamente como gráficos de séries temporais (
e
em função de
) e com um gráfico no plano de fase (v versus
). Faça simulações para ambas as molas (a) linear (
) e (b) não-linear (
).
Passo 1
O enunciado trata de um oscilador massa–mola amortecido forçado descrito pela EDO de segunda ordem
com
,
,
,
e condições iniciais
,
onde
. Vamos resolver numericamente esse sistema (método de Runge–Kutta de
ª ordem), obtendo
,
e traçando:
) Séries temporais
e
.
) Plano de fase
versus
.
Faremos duas simulações: (a) mola linear ( ,
) e (b) mola não linear (
,
).
Faremos duas simulações: (a) mola linear (
Passo 2
Transformamos em sistema de primeira ordem definindo
Então
Substituímos
,
,
,
, fixamos o passo de tempo
pequeno (por ex.
) e intervalo total
.
Passo 3
Implementamos o método de Runge–Kutta de
ª ordem para o vetor
. Em pseudocódigo Python:
import numpy as np
# parâmetros
m\(=1.0\); c\(=0.4\); P\(=0.5\); omega\(=0.5\)
def deriv(x, t, k\(1\), k\(3\)):
x\(1\), x\(2=\) x
dx\(1=\) x\(2\)
dx\(2=\) (P*np.cos(omega*t) - c*x\(2\) - k\(1\)*x\(1\) - k\(3\)*x\(1\)**\(3\))/m
return np.array([dx\(1\), dx\(2\)])
def rk\(4\)(f, x\(0\), t, k\(1\), k\(3\)):
n \(=\) len(t)
x \(=\) np.zeros((n,\(2\)))
x[\(0\)] \(=\) x\(0\)
for i in range(n-\(1\)):
dt \(=\) t[i\(+1\)]-t[i]
\(k_1 =\) f(x[i], t[i], k\(1\), k\(3\))
\(k_2 =\) f(x[i]\(+0.5\)*dt*\(k_1\), t[i]\(+0.5\)*dt, k\(1\), k\(3\))
\(k_3 =\) f(x[i]\(+0.5\)*dt*\(k_2\), t[i]\(+0.5\)*dt, k\(1\), k\(3\))
\(k_4 =\) f(x[i]\(+\)dt*\(k_3\), t[i]\(+\)dt, k\(1\), k\(3\))
x[i\(+1\)] \(=\) x[i] \(+\) dt*(\(k_1 +2\)*\(k_2 +2\)*\(k_3 + k_4\))/\(6\)
return x
t \(=\) np.arange(\(0,200\),\(0.01\))
x\(0=\) [\(0.0\), \(0.0\)]
# simulação linear
\(sol_linear =\) rk\(4\)(deriv, x\(0\), t, k\(1=1.0\), k\(3=0.0\))
# simulação não linear
\(sol_nl =\) rk\(4\)(deriv, x\(0\), t, k\(1=1.0\), k\(3=0.5\))
# para plotar, use matplotlib:
# plt.figure(); plt.plot(t, \(sol_linear\)[:,\(0\)]); plt.title('x(t) linear')
# plt.figure(); plt.plot(t, \(sol_linear\)[:,\(1\)]); plt.title('v(t) linear')
# plt.figure(); plt.plot(\(sol_linear\)[:,\(0\)], \(sol_linear\)[:,\(1\)]); plt.title('fase linear')
# e o mesmo para \(sol_nl\).
Passo 4
Para o caso (a) linear (
), os gráficos mostram:
)
e
começam em zero, apresentam um regime transitório amortecido e convergem para uma oscilação quase senoidal de frequência
.
) No plano de fase, a curva tende a uma elipse de orientação inclinada (regime estacionário linear amortecido-forçado).
Passo 5
Para o caso (b) não linear (
), os gráficos mostram:
)
e
também partem de zero, passam por amortecimento transitório e atingem regime permanente, mas a forma de onda é levemente distorcida (harmônicos gerados pela não linearidade cúbica).
) O retrato de fase não é uma elipse perfeita: surgem achatamentos e distorções características da não linearidade, indicando troca de energia entre modos.
Resposta
Os resultados gráficos obtidos com o código acima reproduzem:
(a) Mola linear (
,
):
-
amortecido em regime transitório e senoidal estacionário em torno de amplitude constante. -
análogo a com defasagem de . - Plano de fase: elipse amortecida até uma órbita fechada elíptica.
(b) Mola não linear (
,
):
-
e transitório semelhante ao linear, mas onda estacionária apresenta distorção (picos mais achatados). - Plano de fase: órbita estacionária distorcida, não elíptica.