Espectro de respuesta no amortiguado

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

Carga rectangular

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

F0 = 4.0
k = 3.73

Rd = np.zeros(1000)
tT = np.linspace(0.01,6,1000)

def dU_dx(U, x, n):
    # Here U is a vector such that y=U[0] and z=U[1]. This function should return [z', u']
    m = 0.0236205
    k = 3.73
    omega = np.power(k/m,1.0/2.0)
    T = (2.0*np.pi)/omega
    F0 = 4.0
    t1 = T*n
    if x <= t1:
        return [U[1],-omega**2*U[0] + F0/m]
    else:
        return [U[1],-omega**2*U[0]]
    
U0 = [0.0,0.0]
xs = np.linspace(0, 3, 3000)
plt.figure(figsize=(19,8.5))

for i in range(tT.size):
    Us = odeint(dU_dx, U0, xs,args=(tT[i],))
    ys = np.abs(Us[:,0])
    plt.plot(xs, ys)
    Rd[i] = max(ys)*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.99999941937

Carga Senoidal

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

F0 = 4.0
k = 3.73

Rd = np.zeros(1000)
tT = np.linspace(0.01,6,1000)

def dU_dx(U, x, n):
    # Here U is a vector such that y=U[0] and z=U[1]. This function should return [z', u']
    m = 0.0236205
    k = 3.73
    omega = np.power(k/m,1.0/2.0)
    T = (2.0*np.pi)/omega
    F0 = 4.0
    t1 = T*n # 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)
plt.figure(figsize=(19,8.5))

for i in range(tT.size):
    Us = odeint(dU_dx, U0, xs,args=(tT[i],))
    ys = np.abs(Us[:,0])
    plt.plot(xs, ys)
    Rd[i] = max(ys)*k/F0
    
plt.xlabel(r'$t$')
plt.ylabel(r'$u$')
plt.grid(True)
plt.show()
In [5]:
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 [6]:
print 'Rd max =', max(Rd)

Rd max = 1.76844123466

Carga triangular simétrica

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

F0 = 4.0
k = 3.73

Rd = np.zeros(1000)
tT = np.linspace(0.01,6,1000)
    
def dU_dx(U, x, n):
    # Here U is a vector such that y=U[0] and z=U[1]. This function should return [z', u']
    m = 0.0236205
    k = 3.73
    omega = np.power(k/m,1.0/2.0)
    T = (2.0*np.pi)/omega
    F0 = 4.0
    t1 = T*n*0.5 # segundos
    t2 = T*n     # 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)
plt.figure(figsize=(19,8.5))

for i in range(tT.size):
    Us = odeint(dU_dx, U0, xs,args=(tT[i],))
    ys = np.abs(Us[:,0])
    plt.plot(xs, ys)
    Rd[i] = max(ys)*k/F0
    
plt.xlabel(r'$t$')
plt.ylabel(r'$u$')
plt.grid(True)
plt.show()
In [8]:
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 [9]:
print 'Rd max =', max(Rd) 

Rd max = 1.51716174941