Vibración forzada con amortiguamiento

\begin{equation*} m \ddot{u} + c \dot{u} + k u = F_{0} \sin \bigl( \bar{\omega} t \bigr) \end{equation*}

reescribiendo la ecuación

\begin{equation} \ddot{u} + 2 \beta \omega \dot{u} + \omega^{2} u = \frac{F_{0}}{m} \sin \bigl( \bar{\omega} t \bigr) \end{equation}

La solución es

\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] + G_{1} \sin \bigl( \bar{\omega} t \bigr) + G_{2} \cos \bigl( \bar{\omega} t \bigr) \end{equation*}

Coeficientes

\begin{align*} \gamma &= \frac{\bar{\omega}}{\omega} \\ G_{1} &= \frac{F_{0}}{k} \frac{1 - \gamma^{2}}{\bigl( 1 - \gamma^{2} \bigr)^{2} + \bigl( 2 \beta \gamma \bigr)^{2}} \\ G_{2} &= \frac{F_{0}}{k} \frac{-2 \beta \gamma}{\bigl( 1 - \gamma^{2} \bigr)^{2} + \bigl( 2 \beta \gamma \bigr)^{2}} \\ A_{1} &= -\frac{G_{2} \beta \omega + G_{1} \bar{\omega}}{\omega_{d}} = -\frac{G_{2} \beta + G_{1} \gamma}{\sqrt{1 - \beta^{2}}} \\ A_{2} &= -G_{2} \end{align*}

Ejemplo en clase

Usando

\begin{equation*} u = e^{-\beta \omega t} \bigl[ A \sin \bigl( \omega_{d} t \bigr) + B \cos \bigl( \omega_{d} t \bigr) \bigr] + G_{1} \sin \bigl( \bar{\omega} t \bigr) + G_{2} \cos \bigl( \bar{\omega} t \bigr) \end{equation*}

In [1]:
import numpy as np
import matplotlib.pyplot as plt

w = 25.0  # Tn
m = w/981 # Tn/cm/seg^2
k = 7.0  # Tn/cm
u0 = 0.0  # cm
v0 = 0.0 # cm/seg
beta = 0.07
F0 = 2.0
omega_bar = 21.1

omega = np.power(k/m,1.0/2.0)
omega_damping = omega*np.power(1-beta**2,1.0/2.0)
gamma = omega_bar/omega
G1 = (F0/k)*(1 - gamma**2)/((1 - gamma**2)**2 + (2*beta*gamma)**2)
G2 = (F0/k)*(-2*beta*gamma)/((1 - gamma**2)**2 + (2*beta*gamma)**2)
A = -(G2*beta*omega + G1*omega_bar)/omega_damping
B = -G2
theta = np.arctan(-A/B)
a = np.sqrt(A**2 + B**2)
T = (2.0*np.pi)/omega
T_damping = (2.0*np.pi)/omega_damping

print 'omega =', omega
print 'omega_damping =',omega_damping
print 'gamma =',gamma
print 'G1 =', G1
print 'G2 =', G2
print 'A =', A
print 'B =', B
print 'theta =', theta
print 'a =', a
print 'T =', T
print 'T_damping =', T_damping

omega = 16.5734727803
omega_damping = 16.5328179086
gamma = 1.2731188134
G1 = -0.42516875494
G2 = -0.12206314615
A = 0.551186827056
B = 0.12206314615
theta = -1.35285850074
a = 0.564540813377
T = 0.379110967898
T_damping = 0.380043217188
In [2]:
t = np.linspace(0,6,6000)
posicion = np.exp(-beta*omega*t)*(A*np.sin(omega_damping*t) + B*np.cos(omega_damping*t)) + G1*np.sin(omega_bar*t) + G2*np.cos(omega_bar*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 odeint

\begin{align*} u^{\prime} &= z \\ z^{\prime} + 2 \beta \omega z + \omega^{2} u &= \frac{F_{0}}{m} \sin \bigl( \bar{\omega} t \bigr) \end{align*}

sujeto a $u(0) = 0$ y $u^{\prime}(0) = 0$

In [3]:
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 = 25.0/981.0
    k = 7.0
    omega = np.power(k/m,1.0/2.0)
    beta = 0.07
    F0 = 2.0
    omega_bar = 21.1 
    return [U[1],-2.0*beta*omega*U[1] - omega**2*U[0] + (F0/m)*np.sin(omega_bar*x)]
U0 = [0.0,0.0]
#U0 = [6.5,32.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()

Valor máximo aproximado

In [4]:
valor_pico = np.where(posicion >= max(posicion))[0]
print 't pico =', t[valor_pico[0]]
print 'valor pico =', posicion[valor_pico[0]]

t pico = 0.493082180363
valor pico = 0.698120121439

Gráfico interactivo

In [5]:
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 = np.exp(-beta*omega*t)*(A*np.sin(omega_damping*t) + B*np.cos(omega_damping*t)) + G1*np.sin(omega_bar*t) + G2*np.cos(omega_bar*t)
trace = go.Scatter(x = t, y = posicion)
data = [trace]

offline.iplot(data)

Desplazamientos máximos

Usando la transformada de Hilbert

In [6]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert

t = np.linspace(0,6,6000)
posicion = np.exp(-beta*omega*t)*(A*np.sin(omega_damping*t) + B*np.cos(omega_damping*t)) + G1*np.sin(omega_bar*t) + G2*np.cos(omega_bar*t)
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

In [7]:
import numpy as np
from scipy.signal import argrelextrema
import matplotlib.pyplot as plt

t = np.linspace(0,6,6000)
posicion = np.exp(-beta*omega*t)*(A*np.sin(omega_damping*t) + B*np.cos(omega_damping*t)) + G1*np.sin(omega_bar*t) + G2*np.cos(omega_bar*t)

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

In [8]:
for i in range(np.size(x)):
    print i, '\t', round(x[i],3), '\t', round(y[i],3)

0 	0.164 	0.348
1 	0.493 	0.698
2 	0.813 	0.632
3 	1.116 	0.399
4 	1.4 	0.333
5 	1.692 	0.444
6 	1.995 	0.497
7 	2.297 	0.46
8 	2.594 	0.421
9 	2.889 	0.427
10 	3.187 	0.449
11 	3.486 	0.452
12 	3.784 	0.442
13 	4.082 	0.437
14 	4.379 	0.441
15 	4.677 	0.444
16 	4.975 	0.444
17 	5.273 	0.442
18 	5.57 	0.442
19 	5.868 	0.442

Periodo aproximado

In [9]:
for i in range(np.size(x) - 1):
    print i+1, '\t', x[i+1]-x[i]

1 	0.329054842474
2 	0.320053342224
3 	0.303050508418
4 	0.284047341224
5 	0.292048674779
6 	0.303050508418
7 	0.302050341724
8 	0.297049508251
9 	0.295049174862
10 	0.297049508251
11 	0.29904984164
12 	0.298049674946
13 	0.298049674946
14 	0.297049508251
15 	0.298049674946
16 	0.298049674946
17 	0.298049674946
18 	0.297049508251
19 	0.298049674946

Resonancia

Se produce cuando $\gamma = 1$

In [10]:
import numpy as np
import matplotlib.pyplot as plt

w = 25.0  # Tn
m = w/981 # Tn/cm/seg^2
k = 7.0  # Tn/cm
u0 = 0.0  # cm
v0 = 0.0 # cm/seg
beta = 0.07
F0 = 2.0
omega_bar = np.power(k/m,1.0/2.0)

omega = np.power(k/m,1.0/2.0)
omega_damping = omega*np.power(1-beta**2,1.0/2.0)
gamma = omega_bar/omega
G1 = (F0/k)*(1 - gamma**2)/((1 - gamma**2)**2 + (2*beta*gamma)**2)
G2 = (F0/k)*(-2*beta*gamma)/((1 - gamma**2)**2 + (2*beta*gamma)**2)
A = -(G2*beta*omega + G1*omega_bar)/omega_damping
B = -G2
theta = np.arctan(-A/B)
a = np.sqrt(A**2 + B**2)
T = (2.0*np.pi)/omega
T_damping = (2.0*np.pi)/omega_damping

print 'omega =', omega
print 'omega_damping =',omega_damping
print 'gamma =',gamma
print 'G1 =', G1
print 'G2 =', G2
print 'A =', A
print 'B =', B
print 'theta =', theta
print 'a =', a
print 'T =', T
print 'T_damping =', T_damping

omega = 16.5734727803
omega_damping = 16.5328179086
gamma = 1.0
G1 = 0.0
G2 = -2.04081632653
A = 0.143208434382
B = 2.04081632653
theta = -0.0700572930881
a = 2.04583477688
T = 0.379110967898
T_damping = 0.380043217188
In [11]:
t = np.linspace(0,6,6000)
posicion = np.exp(-beta*omega*t)*(A*np.sin(omega_damping*t) + B*np.cos(omega_damping*t)) + G1*np.sin(omega_bar*t) + G2*np.cos(omega_bar*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()

Amplitud

\begin{equation*} \sqrt{G_{1} + G_{2}} = \frac{F_{0}}{k} \frac{1}{\sqrt{ \bigl( 1 - \gamma^{2} \bigr)^{2} + \bigl( 2 \beta \gamma \bigr)^{2}}} \end{equation*}

Transmisibilidad

\begin{equation*} T = \frac{\sqrt{G_{1} + G_{2}}}{\frac{F_{0}}{k}} = \frac{1}{\sqrt{ \bigl( 1 - \gamma^{2} \bigr)^{2} + \bigl( 2 \beta \gamma \bigr)^{2}}} \end{equation*}

In [12]:
gamma = np.linspace(0,6,6000)
T1 = 1/np.sqrt((1-gamma**2)**2 + (2*0.05*gamma)**2)
T2 = 1/np.sqrt((1-gamma**2)**2 + (2*0.1*gamma)**2)
T3 = 1/np.sqrt((1-gamma**2)**2 + (2*0.2*gamma)**2)
T4 = 1/np.sqrt((1-gamma**2)**2 + (2*0.3*gamma)**2)

plt.figure(figsize=(11,8.5))
plt.plot(gamma,T1, label=r'$\beta = 5\%$')
plt.plot(gamma,T2, label=r'$\beta = 10\%$')
plt.plot(gamma,T3, label=r'$\beta = 20\%$')
plt.plot(gamma,T4, label=r'$\beta = 30\%$')
plt.legend()
plt.xlabel(r'$\gamma$')
plt.ylabel(r'$\frac{\sqrt{G_{1} + G_{2}}}{\frac{F_{0}}{k}}$')
plt.grid(True)
plt.show()