In [1]:
%matplotlib inline
from scipy import *
from IPython.display import SVG, display, Latex
import matplotlib
import matplotlib.pylab as pl
In [2]:
%config InlineBackend.figure_format = 'svg'
import json ; s = json.load( open("mplrc.json") ) ; del json
matplotlib.rcParams.update(s)
matplotlib.rcParams['figure.figsize'] = 10,4
In [3]:
from IPython.core.display import HTML
def css_styling():
styles = open("custom.css", "r").read()
return HTML(styles)
css_styling()
Out[3]:
We have to study the structure in the figure below (btw, one mass is yellow and other is blue because they're M&M's).
In [4]:
display(SVG(filename="mdof.svg"))
We'll use the Principle of Virtual Displacements, and we start drawing the diagrams of bending moments for the three different load conditions we have loading the structure with a unit force in the location and direction of each one of the dynamic degrees of freedom.
In [5]:
display(SVG(filename="bending.svg"))
The bending moments on each stretch of beam can be represented using polynomials,
so we give a short name to the polynomial object that is defined by the scipy
library, and then fill a table m with NDOF rows and as many columns as the number of different intervals of definition of the bending moments (here 6).
Note that the single argument needed to create a
polynomial object is a sequence, containing all the coefficients from the highest power
to the lowest one --- e.g., in p((1,0)) the sequence (1,0) that is used to create
the polynomial object has length 2, so the polynomial is 1*x**1 + 0*x**0.
In [6]:
p = poly1d
m = [[p((1,0)),p((0.0,-1)),p((-0.50,0)),p((-0.5, -0.5)),p((0.50,0)),p((1.0,0))],
[p((0,0)),p((1.0,+0)),p((+0.50,0)),p(( 0.5, 0.5)),p((0.50,0)),p((1.0,0))],
[p((0,0)),p((0.5,+0)),p((+0.75,0)),p((-0.25,0.75)),p((0.25,0)),p((0.5,0))]]
Now that we have a data structure with a representation of the bending moments we
i and j terms of the flexibility matrix
- i and j are the indices used to access the rows of msms accordinglys to the element i,j of the flexibility matrix.
In [7]:
F = zeros((3,3))
l = [1,1,1,1,2,1]
for i in range(3):
for j in range(3):
s = 0
for n in range(6):
integrand = m[i][n]*m[j][n]
integral = polyint(integrand)
s = s + integral(l[n]) - integral(0)
F[i,j] = s
print F*12,"/ 12"
Now we have F, that is an array object in the speech of scipy,
but it's easier to do what we want to do if we transform it into a matrix object,
that has a lot of useful properties that directly mimics the properties of a matrix
as we know them from algebra.
In [8]:
F = matrix(F)
The three most important properties are, the .I attribute that gives
the matrix inverse, the .T attribute that gives the matrix transpose and last but not least the overloading of the * operator to have matrix product between two matrices.
So we compute K by inversion of F, we construct the mass matrix M directly and
at his point we can import a library function and compute the eigenvalues' and eigenvectors' arrays.
In [9]:
K = F.I
print K*304/3
M = matrix("1 0 0;0 1 0;0 0 1")
from scipy.linalg import eigh
evals, evecs = eigh(K,M)
evecs = matrix(evecs)
evecs[:,0] *= -1
evecs[:,2] *= -1
Are the eignvalues correctly sorted? Are the eigenvectors orthonormal with respect to the structural matrices?
In [10]:
print evals
print
print evecs
print
print evecs.T*M*evecs
print
print evecs.T*K*evecs
The largest absolute off-diagonal term in evecs.T*M*evecs is rather small:
if 1 were an AU (iu.e., the mean Earth-Sun distance) then 1.39E-16*would
be about 21 micron.
Our plan is to use an adimensional time, $\alpha=\omega_ot$, and an adimensional load, $\bm\rho(\alpha)=p(\alpha)/p_o\begin{Bmatrix}1\\0\\0\end{Bmatrix}$.
In [11]:
def r(a):
return where(a<4*pi, sin(0.5*a), 0.0)
In [12]:
a = linspace(0,6*pi, 2401)
pl.plot(a, r(a))
pl.xlabel(r'$\alpha = \omega_ot$')
pl.ylabel(r'$p(t)/p_o$')
n = 200
pl.xticks(a[::n], [r"$%5.2f\pi$"%(x/pi,) for x in a[::n]])
pl.xlim(0,6*pi)
pl.ylim(-1.05,1.05);
Note that the circular frq. of the excitation, $0.5\omega_o$, is quite close to the circular frequencies of the first two modes, especially the 1st one.
In [13]:
print " modal circular frequencies:", sqrt(evals)
print "dynamic amplification factors:", 1/(1-0.25/evals)
With $M_i = m$, $K_i = \omega_i^2 m = \Lambda_i^2\omega_o^2 m$ and $p_o = \delta_o \tfrac{EJ}{L^3}$, the modal eom is
$$ m\ddot q_i(t) + \Lambda^2\omega_o^2 m q_i(t) = \Gamma_i \frac{EJ}{L^3} \delta_o \sin (0.5 \omega_o t)$$dividing both members by $m$, as by definition it is $\omega = \tfrac{EJ}{mL^3}$ we have
$$ \ddot q_i(t) + \Lambda_i^2\omega_o^2 q_i(t) = \Gamma_i \omega_o^2 \delta_o \sin (0.5 \omega_o t).$$The coefficients $\Gamma_i$ are given by
$$\bm\Gamma = \bm\Psi^\text{T}\begin{Bmatrix}1\\0\\0\end{Bmatrix}.$$
In [14]:
Gamma = evecs.T*matrix('1;0;0')
print evecs
print
print Gamma
For our undamped systems, subjected to a sine excitation not in resonance, with $\Lambda_o=0.5$, the general integrals can be written
$$ q_i(a) = A_i\sin(\Lambda_ia) + B_i\cos(\Lambda_ia) + C_i\sin(\Lambda_oa).$$With $\xi_i = C_i \sin(\Lambda_o a)$ and $\ddot\xi_i = - \omega_0^2\Lambda_o^2 C_i \sin(\Lambda_o a)$, substituting in the eom gives
$$ C_i\,\omega_o^2\left(\Lambda_i^2-\Lambda_o^2\right) \sin(\Lambda_o a) = \Gamma_i \omega_o^2 \delta_o \sin(\Lambda_o a) \qquad\Rightarrow\qquad C_i = \frac{\Gamma_i}{\Lambda_i^2-\Lambda_o^2}\delta_o.$$
In [15]:
L_o = 0.50
L = sqrt(evals)
C = ravel(Gamma)/(evals-L_o**2)
print C
We still miss the coefficients $\beta_i = \frac{\Lambda_o}{\Lambda_i}$, so we compute them and immediately define two arrays of functions, the modal velocities and the modal velocities. Note that the modal displacements are normalized with respect to $\delta_o$, while the velocities are normalized with respect to the unit velocity $\delta_o\omega_o$.
In [16]:
B = L_o/sqrt(evals)
q0_f = [lambda a,i=i: C[i] * (sin(L_o*a) - B[i]*sin(L[i]*a)) for i in range (3)]
q1_f = [lambda a,i=i: C[i] * L_o * (cos(L_o*a) - cos(L[i]*a)) for i in range (3)]
It's about time to display our results, first in a textual representation
In [17]:
format_q = r'$\displaystyle\frac{q_%d(a)}{\delta_o} = %+10G \sin( 0.500 a) %+10G \sin(%9G a)$'
format_v = r'$\displaystyle\frac{\dot{q}_%d(a)}{\delta_o\omega_o} = %+10G \left(\cos( 0.500 a) - \cos(%9G a)\right)$'
for i in range(3):
display(Latex(format_q % (i+1, C[i], -C[i]*B[i], L[i])))
display(Latex(format_v % (i+1, C[i]*L_o, L[i])))
print
and eventually in a graphical rendering.
Here I have chosen to display separately, mode by mode, the modal displacement
and the modal velocity, so the plotting commands are inside a loop, as well as a final
pl.show() command that is necessary to obtain three separate plots (try to comment
the last line and re-execute).
In [18]:
a = linspace(0,4*pi, 1601)
n = 200
for i in range(3):
pl.plot(a,q0_f[i](a), label=r'$q_%d(a)$'%(i+1,))
pl.plot(a,q1_f[i](a), label=r'$\dot{q}_%d(a)$'%(i+1,))
pl.legend(loc=0)
pl.xlabel(r'$\alpha = \omega_ot$')
pl.ylabel(r'$q_%d(t)/\delta_o,\quad \dot q_%d(t)/(\delta_o\omega_o)}$'%(i+1,i+1))
pl.xticks(a[::n], [r"$%5.2f\pi$"%(x/pi,) for x in a[::n]])
pl.xlim(0,4*pi)
pl.show();
For plotting the DOF displacements, I put inside the loop only the pl.plot command,
so that the three components can be appreciated side by side.
In [19]:
for i in range(3):
pl.plot(a,
evecs[i,0]*q0_f[0](a)+
evecs[i,1]*q0_f[1](a)+
evecs[i,2]*q0_f[2](a),
label=r'$x_%d(a)$'%(i+1,))
pl.xlabel(r'$\alpha = \omega_ot$')
pl.ylabel(r'$x_i(t)/\delta_o$')
pl.xticks(a[::n], [r"$%5.2f\pi$"%(x/pi,) for x in a[::n]])
pl.xlim(0,4*pi)
pl.legend(loc=8);
In [20]:
a_a = 4*pi # from a in forced range to a in free response range
for i in range(3):
print "%d\t%+12G\t%+12G" % (i+1, q0_f[i](a_a), q1_f[i](a_a))
In [21]:
q = [] ; qf = [] ; pp = []
for i in range(3):
l = L[i]
M = matrix((( sin(l*a_a), + cos(l*a_a)),
( l*cos(l*a_a), - l*sin(l*a_a))))
xvT = matrix(( (q0_f[i](a_a),), (q1_f[i](a_a),)))
Af, Bf = ravel(M.I*xvT)
def qf0(a, i=i, l=l, Af=Af, Bf=Bf):
return Af*sin(l*a)+Bf*cos(l*a)
qf.append(qf0)
qf1 = lambda a: l*Af*cos(l*a) - l*Bf*sin(l*a)
print "%d\t%+12G\t%+12G" % (i+1, qf0(4*pi), qf1(4*pi))
format_qf = r'\frac{q_%d(a)}{\delta_o} &= %+10G \sin(%10G a) %+10G \cos(%10G a)'
pp.append(format_qf%(i+1, Af, L[i], Bf, L[i]))
def qi(a, i=i):
return where(a<a_a, q0_f[i](a), qf[i](a))
q.append(qi)
print
display(Latex(r'\begin{align}'+r'\\'.join(pp)+r'\end{align}'))
To close this soo loong exercise, here are the plots for the modal displacements (note the scaling factor applied to $q_3$, that would be otherwise represented as a horizontal line)
In [22]:
a = linspace(0,8*pi,3201)
for i in range(3):
pl.plot(a,(1,1,100)[i]*q[i](a), label=r'$%sq_%d(a)$'%(('','','100\\,')[i],i+1))
pl.xlabel(r'$\alpha = \omega_ot$')
pl.ylabel(r'$q_i(t)/\delta_o$')
pl.xticks(a[::n], [r"$%5.2f\pi$"%(x/pi,) for x in a[::n]])
pl.xlim(0,8*pi)
pl.legend(loc=0,ncol=2);
and for the DOF displacements.
In [23]:
for i in range(3):
pl.plot(a,evecs[i,0]*q[0](a)+evecs[i,1]*q[1](a)+evecs[i,2]*q[2](a),
label=r'$x_%d(a)$'%(i+1,))
pl.xlabel(r'$\alpha = \omega_ot$')
pl.ylabel(r'$x_i(t)/\delta_o$')
pl.xticks(a[::n], [r"$%5.2f\pi$"%(x/pi,) for x in a[::n]])
pl.xlim(0,8*pi)
pl.legend(loc=3,ncol=2);