In [1]:
from sympy import *
%config InlineBackend.figure_format = 'svg'
from IPython.core.display import display, Image, Math, Latex
rcParams['savefig.dpi'] = 120
ndof = 3


# 06 Three Degrees of Freedom System¶

In [2]:
Image(filename='figures/structure.png')

Out[2]:

First we need a number of symbols, to be used in the following.

Note that the symbols function is rather smart when dealing with pedices and greek letters...

In [3]:
x, t, k, m, EJ, L, wn, w0, w, delta, Lambda = symbols(
'x, t, k, m, EJ, L, omega_n, omega_0, omega, delta, Lambda')
w2, w20, w2n = map(lambda xxx: xxx*xxx, (w,w0,wn))
psi = Matrix(symbols('psi_1 psi_2 psi_3'))
A, B, C = symbols('A B C')
beta = symbols('beta_n')
influence = Matrix((1, 0, 1))


Next, using some of the symbols previously defined we assing a name to an algebraic expression, so defining the ground displacement and later, by double differentiation with respect to t, the ground acceleration

In [4]:
ug = delta*sin(w*t) ; ag = diff(ug,t,2)
display(Eq(symbols('u_g'),ug), Eq(symbols('a_g'),ag))

$$u_{g} = \delta \operatorname{sin}\left(\omega t\right)$$
$$a_{g} = - \delta \omega^{2} \operatorname{sin}\left(\omega t\right)$$

The flexibility matrix for the 4 degrees of freedom system is

In [5]:
F0 = Matrix(((  24,  16,   9, -32),
(  16,  16,   5, -33),
(   9,   5,   4, -10),
( -32, -33, -10,  72)))/6


The flexibility matrix for the 3 DOF system can be computed using the "slice" notation, the stiffness matrix is computed inverting the flexibility and the mass matrix is simply the unit matrix. All structural matrices are normalized.

While we are at it, we derive the matrix that gives the reaction $f_4$ as a linear combination of the nodal loads.

In [6]:
F = Matrix(F0[:-1,:-1] - F0[:-1,-1:]*F0[-1:,-1:].inv()*F0[-1:,:-1])
R4_from_fs = F0[-1:,-1:].inv()*F0[-1:,:-1]

K = F.inv() ; M = eye(3)
display(Math(r'\mathbf{F} = \frac1{432k}' + latex(F*432)))
display(Math(r'\mathbf{K} = \frac{k}{209}' + latex(K*209)))
display(Math(r'\mathbf{M} = m' + latex(M)))

$$\mathbf{F} = \frac1{432k}\left(\begin{smallmatrix}704 & 96 & 328\\96 & 63 & 30\\328 & 30 & 188\end{smallmatrix}\right)$$
$$\mathbf{K} = \frac{k}{209}\left(\begin{smallmatrix}912 & -684 & -1482\\-684 & 2064 & 864\\-1482 & 864 & 2928\end{smallmatrix}\right)$$
$$\mathbf{M} = m\left(\begin{smallmatrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{smallmatrix}\right)$$

Now we can write the equation of free vibrations, denoting by Free the matrix of coefficients of the equation, we can write the scalar equations

In [7]:
Free =  K - Lambda * M
dummy = [display(Math(latex(line)+'=0')) for line in Free*psi]

$$\psi_{1} \left(- \Lambda + \frac{48}{11}\right) - \frac{36}{11} \psi_{2} - \frac{78}{11} \psi_{3}=0$$
$$- \frac{36}{11} \psi_{1} + \psi_{2} \left(- \Lambda + \frac{2064}{209}\right) + \frac{864}{209} \psi_{3}=0$$
$$- \frac{78}{11} \psi_{1} + \frac{864}{209} \psi_{2} + \psi_{3} \left(- \Lambda + \frac{2928}{209}\right)=0$$

To have non trivial solution, we need that the determinant of Free is zero

In [8]:
char_eq = Free.det().expand()
display(Eq(char_eq))

$$- \Lambda^{3} + \frac{5904}{209} \Lambda^{2} - \frac{34380}{209} \Lambda + \frac{15552}{209} = 0$$

The eigenvalues are normalized with respect to $\omega_0^2 = k/m$, that is $\omega_i^2 = \Lambda_i × \omega_0^2$.

The eigenvalues are the solutions of the characteristic equation,

In [9]:
Lambdas = [(r).evalf().subs(I,0) for r in solve(char_eq)]
Lambdas.sort()
Omegas = map(sqrt, Lambdas)
Omegas = [om.evalf(n=8) for om in Omegas]

display(Math(r'\begin{align*}' +
',\\\ '.join([
r'\frac{\omega^2_%d}{k/m} &='%(i+1,)+str(Lambdas[i]) for i in range(ndof)]) +
r'.\end{align*}'))

\begin{align*}\frac{\omega^2_1}{k/m} &=0.493438135532393,\\ \frac{\omega^2_2}{k/m} &=7.41331671942034,\\ \frac{\omega^2_3}{k/m} &=20.3420489727985.\end{align*}

Now we solve the equation of free vibration once for each eigenvalue, to obtain the corresponding eigenvectors.

In [10]:
Psi = eye(3)
for i, r in enumerate(Lambdas):
equ = [Eq(r,) for r in (Free.subs(Lambda, r)*psi)[1:]]
print "The linearly indipendent equations of free vibrations for i=%d"%(i+1,)
display(equ)
d = solve(equ,psi[1:]) ; d[psi[0]] = 1
unnorm = psi.subs(d)
factor2 = (unnorm.transpose()*M*unnorm)[0]
print "The non-normalized eigenvector and the modal mass for i=%d"%(i+1,)
display((unnorm, factor2))
Psi[:,i] = (unnorm/sqrt(factor2))
print "---------------------------------------------------------\n"
print "The matrix of normalized eigenvectors:"
display(Psi)

The linearly indipendent equations of free vibrations for i=1


$$\begin{bmatrix}- \frac{36}{11} \psi_{1} + 9.38215995059201 \psi_{2} + \frac{864}{209} \psi_{3} = 0, & - \frac{78}{11} \psi_{1} + \frac{864}{209} \psi_{2} + 13.516131242458 \psi_{3} = 0\end{bmatrix}$$
The non-normalized eigenvector and the modal mass for i=1


$$\begin{pmatrix}\left(\begin{smallmatrix}1\\0.13599061023812\\0.483032288981684\end{smallmatrix}\right), & 1.25181363827182\end{pmatrix}$$
---------------------------------------------------------

The linearly indipendent equations of free vibrations for i=2


$$\begin{bmatrix}- \frac{36}{11} \psi_{1} + 2.46228136670406 \psi_{2} + \frac{864}{209} \psi_{3} = 0, & - \frac{78}{11} \psi_{1} + \frac{864}{209} \psi_{2} + 6.59625265857009 \psi_{3} = 0\end{bmatrix}$$
The non-normalized eigenvector and the modal mass for i=2


$$\begin{pmatrix}\left(\begin{smallmatrix}1\\9.11190262704693\\-4.63557664727324\end{smallmatrix}\right), & 105.51534033753\end{pmatrix}$$
---------------------------------------------------------

The linearly indipendent equations of free vibrations for i=3


$$\begin{bmatrix}- \frac{36}{11} \psi_{1} - 10.4664508866741 \psi_{2} + \frac{864}{209} \psi_{3} = 0, & - \frac{78}{11} \psi_{1} + \frac{864}{209} \psi_{2} - 6.33247959480803 \psi_{3} = 0\end{bmatrix}$$
The non-normalized eigenvector and the modal mass for i=3


$$\begin{pmatrix}\left(\begin{smallmatrix}1\\-1.0172645823143\\-1.78385915048036\end{smallmatrix}\right), & 5.21698069918361\end{pmatrix}$$
---------------------------------------------------------

The matrix of normalized eigenvectors:


$$\left(\begin{smallmatrix}0.893779029279258 & 0.0973513992186476 & 0.437814746086514\\0.121545555609721 & 0.88705647028709 & -0.445373434808739\\0.431724130356588 & -0.451279872797338 & -0.780999841021665\end{smallmatrix}\right)$$
In [11]:
display(Latex(r'''The modal loads
$$\mathbf{p}^* = -\Psi^T M E a_g ='''+latex(Psi.transpose()*M*influence)+"\\;"+latex(-ag)+'$$'))

The modal loads $$\mathbf{p}^* = -\Psi^T M E a_g =\left(\begin{smallmatrix}1.32550315963585\\-0.35392847357869\\-0.343185094935151\end{smallmatrix}\right)\;\delta \omega^{2} \operatorname{sin}\left(\omega t\right)$$
In [12]:
q123 = []
for i in range(3):
wi = symbols('omega_%d'%(i+1,))
betai = symbols('beta_%d'%(i+1,))
qi = symbols('q_%d'%(i+1,))
display(Latex("The particular integral for mode %d, \
$\\xi(t)=C\\times \\sin(\\omega t)$ and the corresponding equation of equilibrium" % (i+1,)))
xi = C*sin(w*t)
p  = -(Psi[:,i].transpose()*influence*ag)[0]
display(Eq(xi*symbols('omega_%d'%(i+1,))**2+xi.diff(t,2),p))
xi = xi.subs(C, solve(xi.diff(t,2)+wn*wn*xi-p,C)[0])
# homogeneous solution, displacements and velocities
h_d = A*sin(wn*t) + B*cos(wn*t)
display(Latex(r'$$q_{%d}(t) = '%(i+1,) + latex((h_d+xi).subs(wn,wi)) + "$$"))
h_v = h_d.diff(t)
d = solve(map(lambda x: x.subs(t,0),(h_d + xi, h_v + xi.diff(t))), A,B)
display(Latex(r'''Imposing the rest initial conditions we have
\begin{align*} A & = %s, & (%s & = %s)\\ B & = %s. \end{align*}''' % (
latex(d[A].subs(w,wn*betai).simplify()), latex(betai), latex(w/wi),
latex(d[B]))))
d[w] = 2.75 ; d[w2n] = wn*wn ; d[wn] = Omegas[i]
q = (h_d+xi).subs(d).simplify()
d[t] = w0*t ; d[delta] = 1
q123.append(q.subs(d))
display(Latex(
'''Substituting the numerical values, expressing the frequencies
as multiples of $\\omega_0$ and'''))
display(Latex('''normalizing with respect to
$\\delta$ we have'''+
r'$$\frac{%s(t)}{\delta} = %s.$$'%(latex(qi),latex(q.subs(d).evalf(n=8)))))
a = display("----+----|"*7+'--') if i<ndof-1 else 1

The particular integral for mode 1, $\xi(t)=C\times \sin(\omega t)$ and the corresponding equation of equilibrium
$$- C \omega^{2} \operatorname{sin}\left(\omega t\right) + C \omega_{1}^{2} \operatorname{sin}\left(\omega t\right) = 1.32550315963585 \delta \omega^{2} \operatorname{sin}\left(\omega t\right)$$
$$q_{1}(t) = A \operatorname{sin}\left(\omega_{1} t\right) + B \operatorname{cos}\left(\omega_{1} t\right) - 1.32550315963585 \frac{\delta \omega^{2} \operatorname{sin}\left(\omega t\right)}{1.0 \omega^{2} - 1.0 \omega_{1}^{2}}$$
Imposing the rest initial conditions we have \begin{align*} A & = 1.32550315963585 \frac{\beta_{1}^{3} \delta}{\beta_{1}^{2} -1}, & (\beta_{1} & = \frac{\omega}{\omega_{1}})\\ B & = 0. \end{align*}
Substituting the numerical values, expressing the frequencies as multiples of $\omega_0$ and
normalizing with respect to $\delta$ we have$$\frac{q_{1}(t)}{\delta} = 5.5513769 \operatorname{sin}\left(0.70245152 \omega_{0} t\right) - 1.4180266 \operatorname{sin}\left(2.75 \omega_{0} t\right).$$
----+----|----+----|----+----|----+----|----+----|----+----|----+----|--

The particular integral for mode 2, $\xi(t)=C\times \sin(\omega t)$ and the corresponding equation of equilibrium
$$- C \omega^{2} \operatorname{sin}\left(\omega t\right) + C \omega_{2}^{2} \operatorname{sin}\left(\omega t\right) = - 0.35392847357869 \delta \omega^{2} \operatorname{sin}\left(\omega t\right)$$
$$q_{2}(t) = A \operatorname{sin}\left(\omega_{2} t\right) + B \operatorname{cos}\left(\omega_{2} t\right) + 0.35392847357869 \frac{\delta \omega^{2} \operatorname{sin}\left(\omega t\right)}{1.0 \omega^{2} - 1.0 \omega_{2}^{2}}$$
Imposing the rest initial conditions we have \begin{align*} A & = 0.35392847357869 \frac{\beta_{2}^{3} \delta}{- \beta_{2}^{2} + 1}, & (\beta_{2} & = \frac{\omega}{\omega_{2}})\\ B & = 0. \end{align*}
Substituting the numerical values, expressing the frequencies as multiples of $\omega_0$ and
normalizing with respect to $\delta$ we have$$\frac{q_{2}(t)}{\delta} = - 18.121208 \operatorname{sin}\left(2.7227407 \omega_{0} t\right) + 17.941582 \operatorname{sin}\left(2.75 \omega_{0} t\right).$$
----+----|----+----|----+----|----+----|----+----|----+----|----+----|--

The particular integral for mode 3, $\xi(t)=C\times \sin(\omega t)$ and the corresponding equation of equilibrium
$$- C \omega^{2} \operatorname{sin}\left(\omega t\right) + C \omega_{3}^{2} \operatorname{sin}\left(\omega t\right) = - 0.343185094935151 \delta \omega^{2} \operatorname{sin}\left(\omega t\right)$$
$$q_{3}(t) = A \operatorname{sin}\left(\omega_{3} t\right) + B \operatorname{cos}\left(\omega_{3} t\right) + 0.343185094935151 \frac{\delta \omega^{2} \operatorname{sin}\left(\omega t\right)}{1.0 \omega^{2} - 1.0 \omega_{3}^{2}}$$
Imposing the rest initial conditions we have \begin{align*} A & = 0.343185094935151 \frac{\beta_{3}^{3} \delta}{- \beta_{3}^{2} + 1}, & (\beta_{3} & = \frac{\omega}{\omega_{3}})\\ B & = 0. \end{align*}
Substituting the numerical values, expressing the frequencies as multiples of $\omega_0$ and
normalizing with respect to $\delta$ we have$$\frac{q_{3}(t)}{\delta} = - 0.2030852 \operatorname{sin}\left(2.75 \omega_{0} t\right) + 0.12382651 \operatorname{sin}\left(4.5102161 \omega_{0} t\right).$$

What is the definition of $x_3$?

In [13]:
display(Math(r'x_3(t) = %s + %s + %s' % tuple(r'\psi_{3%d}\,q_{%d}(t)'%(i,i) for i in range(1,ndof+1))))

$$x_3(t) = \psi_{31}\,q_{1}(t) + \psi_{32}\,q_{2}(t) + \psi_{33}\,q_{3}(t)$$

We

• display the numerical definition of $x_3$ in terms of the $q_i$,
• compute the expression of $x_3$