Vibración libre amortiguada¶
\begin{equation*} m \ddot{u} + c \dot{u} + k u = 0 \end{equation*}
Frecuencia y coeficiente de amortiguación
\begin{align*} \omega &= \sqrt{\frac{k}{m}} \\ \beta &= \frac{c}{c_{\text{crit}}} = \frac{c}{2 m \omega} \end{align*}
Reordenando
\begin{equation} \ddot{u} + 2 \beta \omega \dot{u} + \omega^{2} u = 0 \end{equation}
Suponiendo que la solución particular es
\begin{equation*} u = r \, e^{s t} \end{equation*}
Reemplazando
\begin{equation*} r \, s^{2} \, e^{s t} + 2 \beta \omega \bigl( r \, s \, e^{s t} \bigr) + \omega^{2} \bigl( r \, e^{s t} \bigr) = 0 \end{equation*}
Simplificando
\begin{equation*} s^{2} + 2 \beta \omega s + \omega^{2} = 0 \end{equation*}
Despejando
\begin{equation*} s = \begin{cases} -\beta \omega + \omega \sqrt{1 - \beta^{2}} \\ -\beta \omega - \omega \sqrt{1 - \beta^{2}} \end{cases} \end{equation*}
Si llamamos frecuencia amortiguada a
\begin{equation*} \omega_{d} = \omega \sqrt{1 - \beta^{2}} \end{equation*}
Reescribiendo la solución
\begin{equation*} s = \begin{cases} -\beta \omega + \omega_{d} \\ -\beta \omega - \omega_{d} \end{cases} \end{equation*}
Reemplazando la primera solución
\begin{equation*} u = r \, e^{-\beta \omega t} \, e^{\omega_{d} t} \end{equation*}
Usando la solución de vibración libre sin amortiguación
\begin{equation} u = e^{-\beta \omega t} \bigl[ A_{1} \sin \bigl( \omega_{d} t \bigr) + A_{2} \cos \bigl( \omega_{d} t \bigr) \bigr] \end{equation}
Puede ser reescrito como
\begin{equation} u = a \, e^{-\beta \omega t} \cos \bigl( \omega_{d} t + \phi \bigr) \end{equation}
La amplitud máxima es
\begin{equation*} a = \sqrt{A_{1}^{2} + A_{2}^{2}} \end{equation*}
El ángulo de fase es
\begin{equation*} \phi = \arctan \biggl( -\frac{A_{1}}{A_{2}} \biggr) \end{equation*}
Coeficientes
\begin{align*} A_{1} &= \frac{\dot{u}_{0} + u_{0} \beta \omega}{\omega_{d}} \\ A_{2} &= u_{0} \end{align*}
Ejemplo en clase¶
import numpy as np
import matplotlib.pyplot as plt
w = 42.0 # Tn
m = w/981 # Tn/cm/seg^2
k = 10.5 # Tn/cm
u0 = 6.2 # cm
v0 = 30.0 # cm/seg
beta = 0.08
omega = np.power(k/m,1.0/2.0)
omega_damping = omega*np.power(1-beta**2,1.0/2.0)
A = (v0 + u0*beta*omega)/omega_damping
B = u0
theta = np.arctan(-A/B)
a = np.sqrt(A**2 + B**2)
print 'omega =', omega
print 'omega_damping =',omega_damping
print 'A =', A
print 'B =', B
print 'theta =', theta
print 'a =', a
Usando
\begin{equation*} u = e^{-\beta \omega t} \bigl[ A_{1} \sin \bigl( \omega_{d} t \bigr) + A_{2} \cos \bigl( \omega_{d} t \bigr) \bigr] \end{equation*}
t = np.linspace(0,6,6000)
posicion = np.exp(-beta*omega*t)*(A*np.sin(omega_damping*t) + B*np.cos(omega_damping*t))
plt.figure(figsize=(19,8.5))
plt.plot(t,posicion)
plt.xlabel(r'$t$ seg.')
plt.ylabel(r'$u$ cm.')
plt.grid(True)
plt.show()
Usando
\begin{equation*} u = a \, e^{-\beta \omega t} \cos \bigl( \omega_{d} t + \phi \bigr) \end{equation*}
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(0,6,6000)
posicion = a*np.exp(-beta*omega*t)*np.cos(omega_damping*t + theta)
plt.figure(figsize=(19,8.5))
plt.plot(t,posicion)
plt.xlabel(r'$t$ seg.')
plt.ylabel(r'$u$ cm.')
plt.grid(True)
plt.show()
Usando odeint
\begin{align*} u^{\prime} &= z \\ z^{\prime} + 2 \beta \omega z + \omega^{2} u &= 0 \end{align*}
sujeto a $u(0) = 6.2$ y $u^{\prime}(0) = 30$
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
def dU_dx(U, x):
# Here U is a vector such that y=U[0] and z=U[1]. This function should return [z', u']
m = 42.0/981.0
k = 10.5
beta = 0.08
omega = np.power(k/m,1.0/2.0)
return [U[1],-2*beta*omega*U[1] - omega**2*U[0]]
U0 = [6.2,30.0]
xs = np.linspace(0, 6, 6000)
Us = odeint(dU_dx, U0, xs)
ys = Us[:,0]
plt.figure(figsize=(19,8.5))
plt.plot(xs, ys)
plt.xlabel(r'$t$ seg.')
plt.ylabel(r'$u$ cm.')
plt.grid(True)
plt.show()
Periodo y periodo amortiguado
\begin{align*} T &= \frac{2 \pi}{\omega} \\ T_{d} &= \frac{2 \pi}{\omega_{d}} \end{align*}
T = (2.0*np.pi)/omega
T_damping = (2.0*np.pi)/omega_damping
print 'T =', T
print 'T_damping =', T_damping
Tiempo hasta el ciclo 0
\begin{equation*} t_{0} = \frac{-\phi - \arctan \bigl( \frac{\beta \omega}{\omega_{d}} \bigr)}{\omega_{d}} \end{equation*}
t0 = (-theta-np.arctan(beta*omega/omega_damping))/omega_damping
print 't0 =', t0
tiempo = np.zeros(14)
tiempo[0] = t0
for i in range(1,14):
tiempo[i] = tiempo[i-1] + T_damping
amplitud_maxima = a*np.exp(-beta*omega*tiempo)*np.cos(omega_damping*tiempo + theta)
for i in range(14):
if amplitud_maxima[i] >= 0.009:
print i, '\t', round(tiempo[i],3), '\t', round(amplitud_maxima[i],3)
Valor máximo aproximado
valor_pico = np.where(posicion >= max(posicion))[0]
print 't pico =', t[valor_pico[0]]
print 'valor pico =', posicion[valor_pico[0]]
Gráfico interactivo
import plotly.offline as offline
import plotly.graph_objs as go
import numpy as np
offline.init_notebook_mode()
t = np.linspace(0,6,6000)
posicion = a*np.exp(-beta*omega*t)*np.cos(omega_damping*t + theta)
trace = go.Scatter(x = t, y = posicion)
data = [trace]
offline.iplot(data)
Desplazamientos máximos¶
Usando la transformada de Hilbert
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
t = np.linspace(0,6,6000)
posicion = a*np.exp(-beta*omega*t)*np.cos(omega_damping*t + theta)
analytic_signal = hilbert(posicion)
amplitude_envelope = np.abs(analytic_signal)
instantaneous_phase = np.unwrap(np.angle(analytic_signal))
instantaneous_frequency = (np.diff(instantaneous_phase)/(2.0*np.pi))
fig = plt.figure(figsize=(19,8.5))
ax0 = fig.add_subplot(211)
ax0.plot(t, posicion, label='posicion')
ax0.plot(t, amplitude_envelope, label='envolvente')
ax0.set_xlabel("tiempo en segundos")
ax0.legend()
ax1 = fig.add_subplot(212)
ax1.plot(t[1:], instantaneous_frequency)
ax1.set_xlabel("tiempo en segundos")
plt.show()
Usando argrelextrema
import numpy as np
from scipy.signal import argrelextrema
import matplotlib.pyplot as plt
t = np.linspace(0,6,6000)
posicion = a*np.exp(-beta*omega*t)*np.cos(omega_damping*t + theta)
maximo = argrelextrema(posicion, np.greater)[0] #array of indexes of the locals maxima
x = [t[maximo] for i in maximo][0]
y = [posicion[maximo] for i in maximo][0]
plt.figure(figsize=(19,8.5))
plt.plot(t,posicion)
plt.plot(x, y, 'ro')
plt.xlabel(r'$t$ seg.')
plt.ylabel(r'$u$ cm.')
plt.show()
Desplazamiento máximo aproximado
for i in range(np.size(x)):
print i, '\t', round(x[i],3), '\t', round(y[i],3)
Periodo aproximado
for i in range(np.size(x) - 1):
print i+1, '\t', x[i+1]-x[i]