Vibración no amortiguada con carga impulsiva triangular simétrica

\begin{equation*} \ddot{u} + \omega^{2} u = \begin{cases} \frac{F_{0}}{m t_{1}} t & 0 \leqslant t \leqslant t_{1} \\ -\frac{F_{0}}{m t_{1}} t + \frac{F_{0}}{m} & 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) + \frac{F_{0}}{k} \frac{t}{t_{1}} & 0 \leqslant t \leqslant t_{1} \\ A_{2} \sin \bigl( \omega t \bigr) + B_{2} \cos \bigl( \omega t \bigr) + \frac{F_{0}}{k} \bigl( 1 - \frac{t}{t_{1}} \bigr) & 0 \leqslant t \leqslant t_{1} \\ A_{3} \sin \bigl( \omega t \bigr) + B_{3} \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 \omega t_{1}} \\ B_{1} &= u_{0} \\ A_{2} &= \frac{\dot{u}(t_{1})}{\omega} + \frac{F_{0}}{k \omega t_{1}} \\ B_{2} &= u(t_{1}) - \frac{F_{0}}{k} \\ A_{3} &= \frac{\dot{u}(2 t_{1})}{\omega} \\ B_{3} &= u(2 t_{1}) \end{align*}

Ejemplo en clase

Usando una traslación de coordenadas en la segunga carga impulsiva y en la vibración libre

\begin{equation*} u = \begin{cases} A_{1} \sin \bigl( \omega t \bigr) + B_{1} \cos \bigl( \omega t \bigr) + \frac{F_{0}}{k} \frac{t}{t_{1}} & 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] + \frac{F_{0}}{k} \bigl( 1 - \frac{t - t_{1}}{t_{1}} \bigr) & t_{1} \leqslant t \leqslant 2 t_{1} \\ A_{3} \sin \bigl[ \omega \bigl( t - 2 t_{1} \bigr) \bigr] + B_{3} \cos \bigl[ \omega \bigl( t - 2 t_{1} \bigr) \bigr] & t \geqslant 2 t_{1} \end{cases} \end{equation*}

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

w = 30.0    # Tn
m = w/981.0 # Tn/cm/seg^2
k = 3.0     # Tn/cm
u0 = 0.0    # cm
v0 = 0.0    # cm/seg
F0 = 4.5   # Tn
t1 = 0.3    # segundos
t2 = 0.3    # segundos

omega = np.power(k/m,1.0/2.0)
A1 = (v0/omega) - F0/(k*omega*t1)
B1 = u0
A2 = A1*np.cos(omega*t1) - B1*np.sin(omega*t1) + F0/(k*omega*t1) + F0/(k*omega*t2)
B2 = A1*np.sin(omega*t1) + B1*np.cos(omega*t1) + (F0/k)*(t1/t1) - F0/k
A3 = A2*np.cos(omega*t2) - B2*np.sin(omega*t2) - F0/(k*omega*t2)
B3 = A2*np.sin(omega*t2) + B2*np.cos(omega*t2) + (F0/k)*(1 - t2/t2)

print 'omega =', omega
print 'A1 =', A1
print 'B1 =', B1
print 'A2 =', A2
print 'B2 =', B2
print 'A3 =', A3
print 'B3 =', B3

omega = 9.90454441153
A1 = -0.504818777346
B1 = 0.0
A2 = 1.50715965234
B2 = -0.0855205256599
A3 = -1.97570596448
B3 = 0.339609865429
In [2]:
t = np.linspace(0,3,3000)
posicion1 = np.where(t<=t1, A1*np.sin(omega*t) + B1*np.cos(omega*t) + (F0/k)*(t/t1), 0)
posicion2 = np.where(np.logical_and(t1<t,t<=t1+t2), A2*np.sin(omega*(t - t1)) + B2*np.cos(omega*(t - t1)) + (F0/k)*(1 - (t - t1)/t2), 0)
posicion3 = np.where(t>t1+t2, A3*np.sin(omega*(t - (t1 + t2))) + B3*np.cos(omega*(t - (t1 + t2))), 0)
posicion = posicion1 + posicion2 + posicion3
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) + \frac{F_{0}}{k} \frac{t}{t_{1}} & 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] + \frac{F_{0}}{k} \bigl( 1 - \frac{t - t_{1}}{t_{1}} \bigr) & t_{1} \leqslant t \leqslant 2 t_{1} \\ A_{3} \sin \bigl[ \omega \bigl( t - 2 t_{1} \bigr) \bigr] + B_{3} \cos \bigl[ \omega \bigl( t - 2 t_{1} \bigr) \bigr] & t \geqslant 2 t_{1} \end{cases} \end{equation*}

Coeficientes

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

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

w = 30.0    # Tn
m = w/981.0 # Tn/cm/seg^2
k = 3.0     # Tn/cm
F0 = 4.5   # Tn
t1 = 0.3    # segundos
t2 = 0.3    # segundos

omega = np.power(k/m,1.0/2.0)
A1 = - F0/(k*omega*t1)
A2 = A1*np.cos(omega*t1) - B1*np.sin(omega*t1) + F0/(k*omega*t1) + F0/(k*omega*t2)
B2 = A1*np.sin(omega*t1) + B1*np.cos(omega*t1) + (F0/k)*(t1/t1) - F0/k
A3 = A2*np.cos(omega*t2) - B2*np.sin(omega*t2) - F0/(k*omega*t2)
B3 = A2*np.sin(omega*t2) + B2*np.cos(omega*t2) + (F0/k)*(1 - t2/t2)

print 'omega =', omega
print 'A1 =', A1
print 'A2 =', A2
print 'B2 =', B2
print 'A3 =', A3
print 'B3 =', B3

omega = 9.90454441153
A1 = -0.504818777346
A2 = 1.50715965234
B2 = -0.0855205256599
A3 = -1.97570596448
B3 = 0.339609865429
In [4]:
t = np.linspace(0,3,3000)
posicion1 = np.where(t<=t1, A1*np.sin(omega*t) + (F0/k)*(t/t1), 0)
posicion2 = np.where(np.logical_and(t1<t,t<=t1+t2), A2*np.sin(omega*(t - t1)) + B2*np.cos(omega*(t - t1)) + (F0/k)*(1 - (t - t1)/t2), 0)
posicion3 = np.where(t>t1+t2, A3*np.sin(omega*(t - (t1 + t2))) + B3*np.cos(omega*(t - (t1 + t2))), 0)
posicion = posicion1 + posicion2 + posicion3
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 t_{1}} t & 0 \leqslant t \leqslant 0.3 \\ -\frac{F_{0}}{m \bigl( t_{2} - t_{1} \bigr)} \bigl( t - t_{1} \bigr) + \frac{F_{0}}{m} & 0.3 \leqslant t \leqslant 0.6 \\ 0 & t > 0.6 \end{cases} \end{align*}

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

In [5]:
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 = 30.0/981.0
    k = 3.0
    omega = np.power(k/m,1.0/2.0)
    F0 = 4.5
    t1 = 0.3 # segundos
    t2 = 0.6 # segundos
    if x <= t1:
        return [U[1],-omega**2*U[0] + (F0*x)/(m*t1)]
    elif t1 < x <= t2:
        return [U[1],-omega**2*U[0] - (F0*(x - t1))/(m*(t2 - t1)) + F0/m]
    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 [6]:
T = (2.0*np.pi)/omega
print 'T =', T

T = 0.634373984922
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.430143381127
valor pico = 2.27322032691

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)
posicion1 = np.where(t<=t1, A1*np.sin(omega*t) + (F0/k)*(t/t1), 0)
posicion2 = np.where(np.logical_and(t1<t,t<=t1+t2), A2*np.sin(omega*(t - t1)) + B2*np.cos(omega*(t - t1)) + (F0/k)*(1 - (t - t1)/t2), 0)
posicion3 = np.where(t>t1+t2, A3*np.sin(omega*(t - (t1 + t2))) + B3*np.cos(omega*(t - (t1 + t2))), 0)
posicion = posicion1 + posicion2 + posicion3
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)
posicion1 = np.where(t<=t1, A1*np.sin(omega*t) + (F0/k)*(t/t1), 0)
posicion2 = np.where(np.logical_and(t1<t,t<=t1+t2), A2*np.sin(omega*(t - t1)) + B2*np.cos(omega*(t - t1)) + (F0/k)*(1 - (t - t1)/t2), 0)
posicion3 = np.where(t>t1+t2, A3*np.sin(omega*(t - (t1 + t2))) + B3*np.cos(omega*(t - (t1 + t2))), 0)
posicion = posicion1 + posicion2 + posicion3
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)
posicion1 = np.where(t<=t1, A1*np.sin(omega*t) + (F0/k)*(t/t1), 0)
posicion2 = np.where(np.logical_and(t1<t,t<=t1+t2), A2*np.sin(omega*(t - t1)) + B2*np.cos(omega*(t - t1)) + (F0/k)*(1 - (t - t1)/t2), 0)
posicion3 = np.where(t>t1+t2, A3*np.sin(omega*(t - (t1 + t2))) + B3*np.cos(omega*(t - (t1 + t2))), 0)
posicion = posicion1 + posicion2 + posicion3

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.43 	2.273
1 	1.093 	2.005
2 	1.728 	2.005
3 	2.362 	2.005
4 	2.996 	2.005

Periodo aproximado

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

1 	0.663221073691
2 	0.634211403801
3 	0.634211403801
4 	0.634211403801

Espectro de respuesta

Desplazamientos

\begin{equation*} u = \begin{cases} \frac{F_{0}}{k} \bigl[ \frac{t}{t_{1}} - \frac{1}{\omega t_{1}} \sin \bigl( \omega t \bigr) \bigr] & 0 \leqslant t \leqslant t_{1} \\ ??? & 0 \leqslant t \leqslant t_{1} \\ ??? & t \geqslant 0 \end{cases} \end{equation*}

Valores conocidos

\begin{align*} \omega t &= 2 \pi \frac{t}{T} \\ \sqrt{A_{3}^{2} + B_{3}^{2}} &= ??? \end{align*}

Reemplazando

\begin{equation*} u = \begin{cases} \frac{F_{0}}{k} \bigl[ 1 - \frac{1}{2 \pi \frac{t_{1}}{T}} \sin \bigl( 2 \pi \frac{t}{T} \bigr) \bigr] & \frac{t_{1}}{T} \geqslant 0 \\ ??? & \frac{t_{1}}{T} \geqslant 0 \\ ??? & \frac{t_{1}}{T} \geqslant 0 \end{cases} \end{equation*}

Reescribiendo

\begin{equation*} u = \begin{cases} \frac{F_{0}}{k} \bigl[ \frac{t}{t_{1}} - \frac{1}{2 \pi \frac{t_{1}}{T}} \sin \bigl( 2 \pi \frac{t}{t_{1}} \frac{t_{1}}{T} \bigr) \bigr] & \frac{t_{1}}{T} \geqslant 0 \\ ??? & \frac{t_{1}}{T} \geqslant 0 \\ ??? & \frac{t_{1}}{T} \geqslant 0 \end{cases} \end{equation*}

In [13]:
t = np.linspace(0.01,3,3000)
u0 = np.where(t<=(t1/T), 1 - np.sin(2*np.pi*t)/(2*np.pi*t), np.nan)

plt.figure(figsize=(11,8.5))
plt.plot(t,u0, '--', label='respuesta forzada')
plt.legend(loc='upper right')
plt.xlabel(r'$\frac{t_{1}}{T}$')
plt.ylabel(r'$R_{d}$')
plt.xlim(0,3)
plt.grid(True)
plt.show()