In [1]:
%matplotlib inline
from IPython.display import HTML
HTML(open("00_custom.css", "r").read())
Out[1]:
In [2]:
%matplotlib inline
import numpy as np
from numpy import poly1d as p, polyint as i
from scipy.linalg import eigh
from IPython.display import Latex, display
import matplotlib.pyplot as plt
plt.style.use(['fivethirtyeight', './00_mplrc'])
def P(a,b): return p((a,b))
def print_mat(mat, dl='b', fmt='%f', nl='\\\\'):
return (r'\begin{'+dl+'matrix}' +
nl.join('&'.join(fmt%n for n in row) for row in mat) +
r'\end{'+dl+'matrix}')
We remove a constraint, permitting the vertical displacement of the support in $C$, then we compute the reactions and the ending moments in equilibrium, in turns, with a unit force applied correspondingly a DOF.
We will need
ls with the lengths of the intervals and ms, one list for each unit loading containing the polynomial representation of the bending moments.
In [3]:
ls = (1, 2, 2, 2, 2, 2, 1)
ms = [[P(-2,0),P(0,-2),P(1,-2),P(0,0),P(1,0),P(1,2),P(4,0)],
[P( 0,0),P(0, 0),P(0, 0),P(0,0),P(1,0),P(0,2),P(2,0)],
[P( 2,0),P(0, 2),P(0, 2),P(1,0),P(1,0),P(1,2),P(4,0)]]
Just for peace of conscience, let's print the bending moments.
In [4]:
for nf, force in enumerate(ms, 1):
for ni, interval in enumerate(force, 1):
print ("M_%d, x_%d: %s" % (nf, ni, interval)).replace('\n',' ')
if nf<3 : print
Let's compute the flexibility coefficients of the 3 DOF system... knowing already the results I know that it is convenient to multiply each term by 3
In [5]:
for n1 in (0, 1, 2):
for n2 in (0, 1, 2):
print n1+1, n2+1,
print 3*sum(i(m[0]*m[1])(ls[k])
for k,m in enumerate(zip(ms[n1],ms[n2])))
From the previous results, we can define the flexibility matrix (note that the terms listed in parentheses are eventually divided by 3),
In [6]:
F = np.array(((116,52,40),
(52,36,52),
(40,52,140)))/3.
Latex(r'$$F=\frac{L^3}{3EJ}' +
print_mat(F*3., fmt='%05.1f') + '$$')
Out[6]:
The stiffness matrix follows immediately by inversion
In [7]:
K = np.linalg.inv(F)
Latex(r'$$K=\frac{3EJ}{3196L^3}' +
print_mat(K*3196./3., fmt='%+08.3f') + '$$')
Out[7]:
The 3x3 stiffness matrix must be partitioned accordingly to the structural degrees of freedom and the DOF associated with the imposed displacement.
In [8]:
Kss = K[:2,:2] ; Ksg = K[:2,2:]
Kgs = K[2:,:2] ; Kgg = K[2:,2:]
Mss = np.array((((1.0, 0.0),
(0.0, 1.0))))
display(Latex('$$K_{ss} =\\frac{3EJ}{3196L^3}' +
print_mat(Kss*3196./3.) + '$$'))
display(Latex('$$M_{ss} = m' +
print_mat(Mss) + '$$'))
The influence matrix is simply $E = - K_{ss}^{-1}\,K_{sg}$
In [9]:
E = -np.dot(np.linalg.inv(Kss), Ksg)
Latex('$$E=\\frac1{35}'+print_mat(35*E)+'$$')
Out[9]:
We represent the eigenvalues as multiples of a reference squared frequency $\omega_0^2$.
In [10]:
evl, evc = eigh(Kss, Mss)
display(Latex(r'$$\omega_0^2 = \frac{EJ}{mL^3}.$$'))
display(Latex(r'$$\omega_1^2 = %f\,\omega_0^2,\quad'%evl[0] +
r'\omega_2^2 = %f\,\omega_0^2.$$'%evl[1]))
Latex(r'$$\psi_1 = ' + print_mat(evc[:,0,None], dl='B') + r',\quad' +
r' \psi_2 = ' + print_mat(evc[:,1,None], dl='B', fmt='%+f') + r'.$$')
Out[10]:
The equations of motion can be written in terms of the structural deflections $\boldsymbol x$
$$M\ddot{x} + Kx = -ME\ddot{v}_c = -mE\ddot{v}_c(t)$$Introducing the modal coordinates $\boldsymbol q$, from $x = \Psi q$ and premultiplying term by term by $\Psi^T$ we have
$$ m \ddot{q}_i + m \Lambda_i^2\omega_o^2 q_i = -m\, \psi_i^TE\,\ddot{v}_c(t),\qquad i=1,2.$$Note that we can simplify the unit mass $m$ from the above equation.
In [11]:
print (r'\begin{align}'+
r'\\'.join(r'\ddot{q}_%d + %f\omega_0^2q_%d &= %+f\,\ddot{v}_c(t)'%
(i+1, evl[i], i+1, -np.dot(evc.T,E)[i,0])
for i in (0,1)) + r'\end{align}')
In [12]:
Latex(r'\begin{align}'+
r'\\'.join(r'\ddot{q}_%d + %f\omega_0^2q_%d &= %+f\,\ddot{v}_c(t)'%
(i+1, evl[i], i+1, -np.dot(evc.T,E)[i,0])
for i in (0,1)) + r'\end{align}')
Out[12]:
We have $v_c=v_c(\tau)$, so we have to see what happens when we write it as a function of $t$, but we want to start by plotting $v,\;\dot v$ and $\ddot v$ as a function of $\tau$. First, define our functions
In [13]:
def f0(tau): return np.where(tau<=0.25, tau**3*(tau*(tau*6144 - 3840) + 640), 1)
def f1(tau): return np.where(tau<=0.25, tau**2*(tau*(tau*30720 - 15360) + 1920), 0)
def f2(tau): return np.where(tau<=0.25, tau**1*(tau*(tau*122880 - 46080) + 3840), 0)
then we plot the support displacement, the support velocity and the support acceleration (note that we applied different scale factors to have the different functions in the same plot, and note that deriving with respect to $\tau$ is different than deriving with respect to $t$).
In [14]:
t = np.linspace(0.0, 0.4, 1001)
plt.plot(t, f0(t), label='$x_c$')
plt.plot(t, f1(t)/8, label='$v_c/8$')
plt.plot(t, f2(t)/100, label='$a_c/100$')
plt.legend(loc=4)
plt.xlim((0,0.40))
plt.ylim((-1.05,1.05))
plt.grid(1);
In [15]:
def particular(m, k, c1, c2, c3, tmax):
d = c3 / k
c = c2 / k
b = c1 / k - 6*(c3*m) / k**2
a = -2*(c2*m) / k**2
def xi(t): return t*(t*(d*t+c)+b)+a if t<= tmax else 0
d3 = d*3 ; c2 = c*2 ; b1 = b*1 ; a0 =a*0
def dotxi(t): return t*(d3*t+c2)+b1 if t<= tmax else 0
return xi, dotxi
But $$\frac1{T_1^2} = \frac{\Lambda_1^2}{(2\pi)^2}\,\omega_0^2$$ and the modal equation of motion is hence \begin{align} m \ddot{q}_i(t) + m \Lambda_i^2\omega_o^2 q_i(t) &= -m\,\delta\, \psi_i^TE\, \frac{\Lambda_1^2}{(2\pi)^2}\,\omega_0^2 \left(3840 \frac{t}{T_1} - 46080 \frac{t^2}{T_1^2} + 122880 \frac{t^3}{T_1^3}\right) \end{align}
Using the adimensional time coordinate $\tau = t/T_1$ and taking into account that
\begin{align} \dot f(t) &= \dot f(\tau) d\tau/dt = \dot f(\tau)/T_1 \\ &= \Lambda_1\omega_0\,\dot f(\tau)/2\pi \end{align}eventually we have
\begin{align} m \frac{\Lambda_1^2}{(2\pi)^2}\omega_0^2\ddot{q}_i(\tau) + m \Lambda_i^2\omega_o^2 q_i(\tau) &= -m\,\delta\, \psi_i^TE\, \frac{\Lambda_1^2}{(2\pi)^2}\,\omega_0^2 \left(3840 \tau - 46080 \tau^2 + 122880 \tau^3\right), \end{align}simplifying
\begin{align} \ddot{q}_i(\tau) + (2\pi)^2\frac{\Lambda_i^2}{\Lambda_1^2}q_i(\tau) = \ddot{q}_i(\tau) + a_i^2q_i(\tau) &= -\delta\,\psi_i^TE\,\left(3840 \tau - 46080 \tau^2 + 122880 \tau^3\right),\quad a_i = 2\pi\frac{\Lambda_i}{\Lambda_1}. \end{align}The integral is
\begin{align} q_i(\tau) &= A_i\sin a_i\tau + B_i\cos a_i\tau + \xi_i(\tau),\\ \dot q_i(\tau) &= a_i A_i \cos a_i\tau - a_i B_i \sin a_i\tau + \dot\xi_i(\tau). \end{align}Because the system starts from rest
\begin{align} B_i &= -\xi(0),&A_i &= -\dot\xi(0)/a_i. \end{align}Eventually, the particular integral is given by
\begin{align} \xi_i(\tau) = -\delta\,\psi_i^TE\,\left( \frac{92160}{a_i^4} + \frac{3840a_i^2-737280}{a_i^4} \tau - \frac{46080}{a_i^2} \tau^2 + \frac{122880}{a_i^2} \tau^3 \right). \end{align}
In [16]:
X3, X4, X5 = 640, -3840, 6144
V2, V3, V4 = X3*3, X4*4, X5*5
A1, A2, A3 = V2*2, V3*3, V4*4
A0 = 0
tau_01 = np.linspace(0.00, 0.25, 251)[:,None]
tau_12 = np.linspace(0.25, 1.00, 751)[:,None]
a = 2*np.pi*np.sqrt(evl/evl[0])
print "a =", a
print "a / 2pi =", a/2/np.pi
psiE = evc.T.dot(E).flatten()
print "Psi^T E =", psiE
print -psiE*A3/a**2
print -psiE*A2/a**2
print -psiE*(A1*a**2-6*A3)/a**4
print -psiE*(-2*A2)/a**4
t = tau_01
dxi = A3*t**3/a**2 + A2*t**2/a**2 + t*(A1*a**2-6*A3)/a**4 - 2*A2/a**4
vxi = 3*A3*t**2/a**2 + 2*A2*t/a**2 + (A1*a**2-6*A3)/a**4
axi = 6*A3*t/a**2 + 2*A2/a**2
plt.plot(t, dxi*a**2+axi);
In [17]:
axi = None
dxi = -dxi*psiE ; vxi = -vxi*psiE
B = -dxi[0] ; A = -vxi[0]/a
dq = A*np.sin(a*t) + B*np.cos(a*t) + dxi
vq = a*(A*np.cos(a*t) - B*np.sin(a*t)) + vxi
t2 = np.linspace(0.00, 0.75, 751)[:,None]
B2 = dq[-1] ; A2 = vq[-1]/a
dq2 = A2*np.sin(a*t2) + B2*np.cos(a*t2)
vq2 = a*(A2*np.cos(a*t2)-B2*np.sin(a*t2))
In [18]:
T = np.vstack((t,t2+0.25))
Q = np.vstack((dq,dq2))
plt.plot(T, Q[:,0], label="$q_1(\\tau)$")
plt.plot(T, Q[:,1], label="$q_2(\\tau)$")
plt.legend(loc=1)
plt.xticks(np.linspace(0,1,9));
In [19]:
X = evc.dot(Q.T).T
plt.plot(T, X);