Espectro de respuesta no amortiguado

Ejemplo 4.1 del libro Dynamics of Structures - Anil K. Chopra

Carga rectangular

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

m = 0.0236205
k = 3.73
omega = np.power(k/m,1.0/2.0)
T = (2.0*np.pi)/omega
F0 = 4.0

def rectangular(tau,F0,t1):
    if tau <= t1:
        return F0
    else:
        return 0
    
def dA(tau,omega,F0,t1):
    return rectangular(tau,F0,t1)*np.cos(omega*tau)

def dB(tau,omega,F0,t1):
    return rectangular(tau,F0,t1)*np.sin(omega*tau)

def A(tau,m,omega,F0,t1):
    return integrate.quad(dA, 0, tau, args=(omega,F0,t1))[0]/(m*omega)

def B(tau,m,omega,F0,t1):
    return -integrate.quad(dB, 0, tau, args=(omega,F0,t1))[0]/(m*omega)

Rd = np.zeros(100)
tT = np.linspace(0.01,6,100)
plt.figure(figsize=(19,8.5))

t = np.linspace(0,3,100)
posicion = np.zeros(100)
    
for j in range(tT.size):
    for i in range(t.size):
        posicion[i] = A(t[i],m,omega,F0,T*tT[j])*np.sin(omega*t[i]) + B(t[i],m,omega,F0,T*tT[j])*np.cos(omega*t[i]) 
    plt.plot(t, np.abs(posicion))
    Rd[j] = max(posicion)*k/F0
    
plt.xlabel(r'$t$')
plt.ylabel(r'$u$')
plt.grid(True)
plt.show()
In [2]:
plt.figure(figsize=(19,8.5))
plt.plot(tT, Rd)
plt.xlabel(r'$\frac{t_{1}}{T}$')
plt.ylabel(r'$R_{d}$')
plt.grid(True)
plt.show()
In [3]:
print 'Rd max =', max(Rd)

Rd max = 1.99547197568

Carga senoidal

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

m = 0.0236205
k = 3.73
omega = np.power(k/m,1.0/2.0)
T = (2.0*np.pi)/omega
F0 = 4.0

def sinoide(tau,F0,t1):
    if tau <= t1:
        return F0*np.sin(np.pi*tau/t1)
    else:
        return 0
    
def dA(tau,omega,F0,t1):
    return sinoide(tau,F0,t1)*np.cos(omega*tau)

def dB(tau,omega,F0,t1):
    return sinoide(tau,F0,t1)*np.sin(omega*tau)

def A(tau,m,omega,F0,t1):
    return integrate.quad(dA, 0, tau, args=(omega,F0, t1))[0]/(m*omega)

def B(tau,m,omega,F0,t1):
    return -integrate.quad(dB, 0, tau, args=(omega,F0, t1))[0]/(m*omega)

Rd = np.zeros(100)
tT = np.linspace(0.01,6,100)
plt.figure(figsize=(19,8.5))

t = np.linspace(0,3,100)
posicion = np.zeros(100)
    
for j in range(tT.size):
    for i in range(t.size):
        posicion[i] = A(t[i],m,omega,F0,T*tT[j])*np.sin(omega*t[i]) + B(t[i],m,omega,F0,T*tT[j])*np.cos(omega*t[i]) 
    plt.plot(t, np.abs(posicion))
    Rd[j] = max(posicion)*k/F0
    
plt.xlabel(r'$t$')
plt.ylabel(r'$u$')
plt.grid(True)
plt.show()
In [2]:
plt.figure(figsize=(19,8.5))
plt.plot(tT, Rd)
plt.xlabel(r'$\frac{t_{1}}{T}$')
plt.ylabel(r'$R_{d}$')
plt.grid(True)
plt.show()
In [3]:
print 'Rd max =', max(Rd)

Rd max = 1.76672035413

Carga triangular simétrica

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

m = 0.0236205
k = 3.73
omega = np.power(k/m,1.0/2.0)
T = (2.0*np.pi)/omega
F0 = 4.0

def triangular(tau,F0,t1,t2):
    if tau <= t1:
        return (F0*tau)/t1
    elif t1 < tau <= t2:
        return -(F0*(tau - t1))/(t2 - t1) + F0
    else:
        return 0
    
def dA(tau,omega,F0,t1,t2):
    return triangular(tau,F0,t1,t2)*np.cos(omega*tau)

def dB(tau,omega,F0,t1,t2):
    return triangular(tau,F0,t1,t2)*np.sin(omega*tau)

def A(tau,m,omega,F0,t1,t2):
    return integrate.quad(dA, 0, tau, args=(omega,F0,t1,t2))[0]/(m*omega)

def B(tau,m,omega,F0,t1,t2):
    return -integrate.quad(dB, 0, tau, args=(omega,F0,t1,t2))[0]/(m*omega)

Rd = np.zeros(100)
tT = np.linspace(0.01,6,100)
plt.figure(figsize=(19,8.5))

t = np.linspace(0,3,100)
posicion = np.zeros(100)
    
for j in range(tT.size):
    for i in range(t.size):
        posicion[i] = A(t[i],m,omega,F0,T*tT[j]*0.5,T*tT[j])*np.sin(omega*t[i]) + B(t[i],m,omega,F0,T*tT[j]*0.5,T*tT[j])*np.cos(omega*t[i]) 
    plt.plot(t, np.abs(posicion))
    Rd[j] = max(posicion)*k/F0
    
plt.xlabel(r'$t$')
plt.ylabel(r'$u$')
plt.grid(True)
plt.show()
In [2]:
plt.figure(figsize=(19,8.5))
plt.plot(tT, Rd)
plt.xlabel(r'$\frac{t_{2}}{T}$')
plt.ylabel(r'$R_{d}$')
plt.grid(True)
plt.show()
In [3]:
print 'Rd max =', max(Rd)

Rd max = 1.51697231587