Vibración no amortiguada con carga impulsiva senoidal

\begin{equation*} \ddot{u} + \omega^{2} u = \begin{cases} \frac{F_{0}}{m} \sin \bigl( \pi \frac{t}{t_{1}} \bigr) & 0 \leqslant t \leqslant t_{1} \\ 0 & t \geqslant 0 \end{cases} \end{equation*}

La solución es

\begin{equation*} u = \begin{cases} A_{1} \sin \bigl( \omega t \bigr) + B_{1} \cos \bigl( \omega t \bigr) + G_{1} \sin \bigl( \pi \frac{t}{t_{1}} \bigr) & 0 \leqslant t \leqslant t_{1} \\ A_{2} \sin \bigl( \omega t \bigr) + B_{2} \cos \bigl( \omega t \bigr) & t \geqslant 0 \end{cases} \end{equation*}

Coeficientes

\begin{align*} A_{1} &= \frac{\dot{u}_{0}}{\omega} + \frac{F_{0}}{k} \frac{\pi t_{1}}{\pi^{2} - \omega^{2} t_{1}^{2}}\\ B_{1} &= u_{0} \\ G_{1} &= -\frac{F_{0}}{m} \frac{t_{1}^{2}}{\pi^{2} - \omega^{2} t_{1}^{2}} \\ A_{2} &= \frac{\dot{u}(t_{1})}{\omega} \\ B_{2} &= u(t_{1}) \end{align*}

Ejemplo en clase

Usando una traslación de coordenadas en la vibración libre

\begin{equation*} u = \begin{cases} A_{1} \sin \bigl( \omega t \bigr) + B_{1} \cos \bigl( \omega t \bigr) + G_{1} \sin \bigl( \pi \frac{t}{t_{1}} \bigr) & 0 \leqslant t \leqslant t_{1} \\ A_{2} \sin \bigl[ \omega \bigl( t - t_{1} \bigr) \bigr] + B_{2} \cos \bigl[ \omega \bigl( t - t_{1} \bigr) \bigr] & t \geqslant t_{1} \end{cases} \end{equation*}

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

w = 25.0    # Tn
m = w/981.0 # Tn/cm/seg^2
k = 7.0     # Tn/cm
u0 = 0.0    # cm
v0 = 0.0    # cm/seg
F0 = 2.0    # Tn
t1 = 1.0    # segundos

omega = np.power(k/m,1.0/2.0)
A1 = (v0/omega) + (F0/k)*t1/(np.pi**2 - omega**2*t1**2)
B1 = u0
G1 = -(F0/m)*(np.pi*t1**2)/(np.pi**2 - omega**2*t1**2)
A2 = A1*np.cos(omega*t1) - B1*np.sin(omega*t1) + G1*(np.pi/t1*omega)*np.cos(np.pi*t1/t1)
B2 = A1*np.sin(omega*t1) + B1*np.cos(omega*t1) + G1*np.sin(np.pi*t1/t1)

print 'omega =', omega
print 'A1 =', A1
print 'B1 =', B1
print 'G1 =', G1
print 'A2 =', A2
print 'B2 =', B2

omega = 16.5734727803
A1 = -0.0010789390842
B1 = 0.0
G1 = 0.931051784792
A2 = -48.4764672675
B2 = 0.000821531898881
In [20]:
t = np.linspace(0,3,3000)
posicion = np.where(t<=t1, A1*np.sin(omega*t) + B1*np.cos(omega*t) + G1*np.sin(np.pi*t/t1),
                    A2*np.sin(omega*(t-t1)) + B2*np.cos(omega*(t-t1)))  
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()

Cuando $u_{0} = 0$ y $\dot{u}_{0} = 0$

\begin{equation*} u = \begin{cases} A_{1} \sin \bigl( \omega t \bigr) + G_{1} \sin \bigl( \pi \frac{t}{t_{1}} \bigr) & 0 \leqslant t \leqslant t_{1} \\ A_{2} \sin \bigl[ \omega \bigl( t - t_{1} \bigr) \bigr] + B_{2} \cos \bigl[ \omega \bigl( t - t_{1} \bigr) \bigr] & t \geqslant t_{1} \end{cases} \end{equation*}

Coeficientes

\begin{align*} A_{1} &= \frac{F_{0}}{k} \frac{\pi t_{1}}{\pi^{2} - \omega^{2} t_{1}^{2}}\\ G_{1} &= -\frac{F_{0}}{m} \frac{t_{1}^{2}}{\pi^{2} - \omega^{2} t_{1}^{2}} \\ A_{2} &= \frac{\dot{u}(t_{1})}{\omega} \\ B_{2} &= u(t_{1}) \end{align*}

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

w = 25.0    # Tn
m = w/981.0 # Tn/cm/seg^2
k = 7.0     # Tn/cm
F0 = 2.0    # Tn
t1 = 1.0    # segundos

omega = np.power(k/m,1.0/2.0)
A1 = (v0/omega) + (F0/k)*t1/(np.pi**2 - omega**2*t1**2)
G1 = -(F0/m)*t1**2/(np.pi**2 - omega**2*t1**2)
A2 = A1*np.cos(omega*t1) - B1*np.sin(omega*t1) + G1*(np.pi/t1*omega)*np.cos(np.pi*t1/t1)
B2 = A1*np.sin(omega*t1) + B1*np.cos(omega*t1) + G1*np.sin(np.pi*t1/t1)

print 'omega =', omega
print 'A1 =', A1
print 'G1 =', G1
print 'A2 =', A2
print 'B2 =', B2

omega = 16.5734727803
A1 = -0.0010789390842
G1 = 0.296362987648
A2 = -15.4300619876
B2 = 0.000821531898881
In [22]:
t = np.linspace(0,3,3000)
posicion = np.where(t<=t1, A1*np.sin(omega*t) + G1*np.sin(np.pi*t/t1), A2*np.sin(omega*(t-t1)) + B2*np.cos(omega*(t-t1)))  
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} + \omega^{2} u &= \begin{cases} \frac{F_{0}}{m} \sin \bigl( \pi \frac{t}{t_{1}} \bigr) & 0 \leqslant t \leqslant 1 \\\ 0 & t \geqslant 1 \end{cases} \end{align*}

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

In [23]:
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)
    F0 = 2.0
    t1 = 0.5 # segundos
    if x <= t1:
        return [U[1],-omega**2*U[0] + (F0/m)*np.sin(np.pi*x/t1)]
    else:
        return [U[1],-omega**2*U[0]]
U0 = [0.0,0.0]
xs = np.linspace(0, 3, 3000)
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

\begin{equation*} T = \frac{2 \pi}{\omega} \end{equation*}

In [2]:
T = (2.0*np.pi)/omega
print 'T =', T

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-2-727522ecbefe> in <module>()
----> 1 T = (2.0*np.pi)/omega
      2 print 'T =', T

NameError: name 'omega' is not defined
In [7]:
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.97332444148
valor pico = 2.23681695601

Gráfico interactivo

In [8]:
import plotly.offline as offline
import plotly.graph_objs as go
import numpy as np

offline.init_notebook_mode()

t = np.linspace(0,3,3000)
posicion = np.where(t<=t1, A1*np.sin(omega*t) + B1*np.cos(omega*t) + G1*np.sin(omega_bar*t),
                    A2*np.sin(omega*(t-t1)) + B2*np.cos(omega*(t-t1)))
trace = go.Scatter(x = t, y = posicion)
data = [trace]

offline.iplot(data)

Desplazamientos máximos

Usando la transformada de Hilbert

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

t = np.linspace(0,3,3000)
posicion = np.where(t<=t1, A1*np.sin(omega*t) + B1*np.cos(omega*t) + G1*np.sin(omega_bar*t),
                    A2*np.sin(omega*(t-t1)) + B2*np.cos(omega*(t-t1)))
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 [10]:
import numpy as np
from scipy.signal import argrelextrema
import matplotlib.pyplot as plt

t = np.linspace(0,3,3000)
posicion = np.where(t<=t1, A1*np.sin(omega*t) + B1*np.cos(omega*t) + G1*np.sin(omega_bar*t),
                    A2*np.sin(omega*(t-t1)) + B2*np.cos(omega*(t-t1)))

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

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

0 	0.195 	0.46
1 	0.584 	1.368
2 	0.973 	2.237
3 	1.352 	2.233
4 	1.732 	2.233
5 	2.111 	2.233
6 	2.49 	2.233
7 	2.869 	2.233

Periodo aproximado

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

1 	0.389129709903
2 	0.389129709903
3 	0.379126375458
4 	0.379126375458
5 	0.379126375458
6 	0.379126375458
7 	0.379126375458

Espectro de respuesta

Desplazamientos

\begin{equation*} u = \begin{cases} \frac{F_{0}}{k} \Bigl( \frac{1}{1 - \gamma^{2}} \Bigr) \Bigl[ \sin \bigl( \bar{\omega} t \bigr) - \gamma \sin \bigl( \omega t \bigr) \Bigr] & 0 \leqslant t \leqslant t_{1} \\ A_{2} \sin \bigl( \omega t \bigr) + B_{2} \cos \bigl( \omega t \bigr) & t \geqslant 0 \end{cases} \end{equation*}

Derivando la primera ecuación

\begin{equation*} \cos \bigl( \bar{\omega} t \bigr) - \cos \bigl( \omega t \bigr) = 0 \end{equation*}

La solución es

\begin{equation*} \bar{\omega} t = 2 \pi n \pm \omega t \end{equation*}

Despejando $t$

\begin{equation*} t = \frac{2 \pi n}{\bar{\omega} + \omega} \end{equation*}

Multiplicando por $\bar{\omega}$ y $\omega$

\begin{align*} \bar{\omega} t &= \bar{\omega} \frac{2 \pi n}{\bar{\omega} + \omega} = \frac{2 \pi n}{1 + \frac{1}{\gamma}} = \frac{2 \pi n \gamma}{\gamma + 1} \\ \omega t &= \omega \frac{2 \pi n}{\bar{\omega} + \omega} = \frac{2 \pi n}{\gamma + 1} \end{align*}

Reemplazando

\begin{equation*} u = \begin{cases} \frac{F_{0}}{k} \Bigl( \frac{1}{1 - \gamma^{2}} \Bigr) \Bigl[ \sin \bigl( \frac{2 \pi n \gamma}{1 + \gamma} \bigr) - \gamma \sin \bigl( \frac{2 \pi n}{1 + \gamma} \bigr) \Bigr] & 0 \leqslant t \leqslant t_{1} \\ A_{2} \sin \bigl( \omega t \bigr) + B_{2} \cos \bigl( \omega t \bigr) & t \geqslant 0 \end{cases} \end{equation*}

Valores conocidos

\begin{align*} \frac{F_{0}}{k} \Bigl( \frac{1}{1 - \gamma^{2}} \Bigr) \Bigl[ \sin \bigl( \frac{2 \pi n \gamma}{1 + \gamma} \bigr) - \gamma \sin \bigl( \frac{2 \pi n}{1 + \gamma} \bigr) \Bigr] &= \frac{F_{0}}{k} \bigl( \frac{1}{1 - \gamma} \bigr) \sin \bigl( \frac{2 \pi n \gamma}{1 + \gamma} \bigr) \\ \sqrt{A_{2}^{2} + B_{2}^{2}} &= \frac{F_{0}}{k} \sqrt{ \Bigl( \frac{1}{1 - \gamma^{2}} \Bigr)^{2} \bigl[ \bar{\omega} \cos \bigl( \frac{2 \pi \gamma}{1 + \gamma} \bigr) - \bar{\omega} \cos \bigl( \frac{2 \pi}{1 + \gamma} \bigr) \bigr]^{2} + \Bigl( \frac{1}{1 - \gamma^{2}} \Bigr)^{2} \bigl[ \sin \bigl( \frac{2 \pi \gamma}{1 + \gamma} \bigr) - \gamma \sin \bigl( \frac{2 \pi}{1 + \gamma} \bigr) \bigr]^{2}} \\ &= \frac{F_{0}}{k} \big | 2 \bigl( \frac{\gamma}{1 - \gamma^{2}} \bigr) \cos \bigl( \frac{\pi}{2 \gamma} \bigr) \big | \\ &= \frac{2 F_{0}}{k} \big | \frac{\gamma}{1 - \gamma^{2}} \, \cos \bigl( \frac{\pi}{2 \gamma} \bigr) \big | \end{align*}

Reescribiendo

\begin{equation*} u = \begin{cases} \frac{F_{0}}{k} \bigl( \frac{1}{1 - \gamma} \bigr) \sin \bigl( \frac{2 \pi n \gamma}{1 + \gamma} \bigr) & \frac{t_{1}}{T} \geqslant 0 \\ \frac{2 F_{0}}{k} \big | \frac{\gamma}{1 - \gamma^{2}} \, \cos \bigl( \frac{\pi}{2 \gamma} \bigr) \big | & \frac{t_{1}}{T} \geqslant 0 \end{cases} \end{equation*}

Para encontrar $\frac{t_{1}}{T}$ en función $\gamma$

\begin{equation*} \frac{t_{1}}{T} = \frac{\bar{T}}{2 T} \end{equation*}

Relación de frecuencias

\begin{equation*} \gamma = \frac{\bar{\omega}}{\omega} = \frac{\frac{2 \pi}{\bar{T}}}{\frac{2 \pi}{T}} = \frac{T}{\bar{T}} \end{equation*}

Reemplazando

\begin{equation*} \frac{t_{1}}{T} = \frac{1}{2 \gamma} \end{equation*}

In [13]:
t = np.linspace(0.01,8,8000)
gamma = 1/(2.0*t)
u0 = 2*np.abs((2*t/((2*t)**2 - 1))*np.cos(np.pi*t))
u1 = 1/(1-gamma)*np.sin(2*1*np.pi/(1+2*t))
u2 = np.where(t>=1.5, 1/(1-gamma)*np.sin(2*2*np.pi/(1+2*t)), np.nan)
u3 = np.where(t>=2.5, 1/(1-gamma)*np.sin(2*3*np.pi/(1+2*t)), np.nan)
u4 = np.where(t>=3.5, 1/(1-gamma)*np.sin(2*4*np.pi/(1+2*t)), np.nan)

plt.figure(figsize=(11,8.5))
plt.plot(t,u0, '--', label='respuesta libre')
plt.plot(t,u1, label='respuesta forzada $n=1$')
plt.plot(t,u2, label='respuesta forzada $n=2$')
plt.plot(t,u3, label='respuesta forzada $n=3$')
plt.plot(t,u4, label='respuesta forzada $n=4$')
plt.legend(loc='upper right')
plt.xlabel(r'$\frac{t_{1}}{T}$')
plt.ylabel(r'$R_{d}$')
plt.grid(True)
plt.show()

El espectro de respuesta es la envolvente

\begin{equation*} u_{\max} = \begin{cases} \frac{F_{0}}{k} \frac{1}{1 - \frac{T}{2 t_{1}}} \sin \bigl( \frac{2 n \pi}{1 + 2 \frac{t_{1}}{T}} \bigr) & \frac{t_{1}}{T} \geqslant 0 & \mbox{respuesta forzada para } n = 1, 2, 3, 4, \ldots \\ \frac{2 F_{0}}{k} \big | \frac{2 \frac{t_{1}}{T}}{\bigl( 2 \frac{t_{1}}{T} \bigr)^{2} - 1} \cos \bigl( \pi \frac{t_{1}}{T} \bigr) \big | & \frac{t_{1}}{T} \geqslant 0 & \mbox{respuesta libre} \end{cases} \end{equation*}

In [14]:
t = np.linspace(0.01,8,8000)
gamma = 1/(2.0*t)
u0 = np.where(np.logical_and(0.01 < t, t <= 0.5),2*np.abs((2*t/((2*t)**2 - 1))*np.cos(np.pi*t)), 0)
u1 = np.where(np.logical_and(0.5 < t, t <= 2.5), 1/(1-gamma)*np.sin(2*1*np.pi/(1+2*t)), 0)
u2 = np.where(np.logical_and(2.5 < t, t <= 4.5), 1/(1-gamma)*np.sin(2*2*np.pi/(1+2*t)), 0)
u3 = np.where(np.logical_and(4.5 < t, t <= 6.5), 1/(1-gamma)*np.sin(2*3*np.pi/(1+2*t)), 0)
u4 = np.where(t>=6.5, 1/(1-gamma)*np.sin(2*4*np.pi/(1+2*t)), 0)
u = u0 + u1 + u2 + u3 + u4

plt.figure(figsize=(11,8.5))
plt.plot(t,u)
plt.xlabel(r'$\frac{t_{1}}{T}$')
plt.ylabel(r'$R_{d}$')
plt.grid(True)
plt.show()

Factor de amplificación máximo

In [15]:
valor_pico = np.where(u>= max(u))[0]
print 'gamma =', 1/(2.0*t[valor_pico[0]])
print 'Rd =', u[valor_pico[0]]

gamma = 0.617208695089
Rd = 1.76845761979