The data that describes our system and its loading,
$$\boldsymbol p(t) = p_0\,\begin{Bmatrix}0\\1\end{Bmatrix}\sin 2\omega_0\,t$$
In [2]:
M = np.array(((2.0, 0.0), ( 0.0, 1.0)))
K = np.array(((3.0,-2.0), (-2.0, 2.0)))
p = np.array(( 0.0, 1.0))
w = 2.0
In [3]:
evals, Psi = eigh(K, M)
Mstar = Psi.T@M@Psi
Kstar = Psi.T@K@Psi
pstar = Psi.T@p
The @
operator stands, in this context, for matrix multiplication.
We print the computed values.
In [4]:
print(evals,end='\n\n')
print(Psi,end='\n\n')
print(Mstar,end='\n\n')
print(Kstar,end='\n\n')
print(pstar,end='\n\n')
print(Mstar[0,1]*150*10E6*1000)
The modal equation of motion, divided by $m$ $$ \ddot q_i + \Lambda^2_i\omega_0^2 q_i = \frac{p_0}{k}\frac{k}{m}\,p^\star_i \sin2\omega_0t.$$
The particular integral $\xi_i(t)$ is given by
$$\xi_i = \Delta_{st} \frac{\omega_0^2}{(\Lambda^2_i-4)\omega_0^2}p^\star_i\sin2\omega_0t.$$Eventually the modal response, computed in the hypotesis of initial rest conditions is
$$q_i = \frac{\Delta_{st}}{\Lambda^2_i-4}p^\star_i(\sin2\omega_0t-\beta\sin\Lambda_i\omega_0t).$$
In [6]:
L = np.sqrt(evals)
DAF = 1.0/(L**2-w**2)
beta = w/L
t = np.linspace(0,60,601)[:,None]
q = pstar*DAF*(np.sin(w*t)-beta*np.sin(L*t))
The definition of time vectoris a bit complicated...
In [7]:
curves = plt.plot(t,q)
plt.legend(curves,['q1', 'q2'])
plt.title('Modal Response')
plt.xlabel('$\omega_0t$')
plt.ylabel('$q_i/\Delta_{st}$');
In [8]:
x = (Psi@q.T).T
curves = plt.plot(t, x)
plt.legend(curves,['x1', 'x2'])
plt.title('Structural Response')
plt.xlabel('$\omega_0t$')
plt.ylabel('$X_i/\Delta_{st}$');
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt ; plt.style.use(['fivethirtyeight', '00_mplrc'])
import numpy as np
from scipy.linalg import eigh
np.set_printoptions(suppress=False, linewidth=120)
from IPython.display import HTML
HTML(open('00_custom.css').read())
Out[1]: