In [1]:
from sympy import *
%load_ext sympyprinting
%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$
  • display the full expression of $x_3
In [14]:
display(Psi[2,:]*Matrix(symbols('q_1 q_2 q_3')))
x3 = (Psi[2,:]*Matrix(q123))[0].simplify()
display(Math(r"\frac{x_3(t)}{\delta} = %s"%(latex(x3.evalf(n=6)),)))
$$\left(\begin{smallmatrix}0.431724130356588 q_{1} - 0.451279872797338 q_{2} - 0.780999841021665 q_{3}\end{smallmatrix}\right)$$
$$\frac{x_3(t)}{\delta} = 2.39666 \operatorname{sin}\left(0.70245152 \omega_{0} t\right) + 8.17774 \operatorname{sin}\left(2.7227407 \omega_{0} t\right) - 8.55026 \operatorname{sin}\left(2.75 \omega_{0} t\right) - 0.0967085 \operatorname{sin}\left(4.5102161 \omega_{0} t\right)$$

The x3 expression is not suitable for numerical evaluation, as in plotting, so we have to take an extra step to transform the symbolic expression in a numerical function.

In [15]:
x3t = lambdify(t, x3.subs(w0,1), modules='numpy')

It's time to plot our results

In [16]:
time = linspace(0.0,10.0,201)
plot(time,x3t(time))
xlabel(r'$\omega_0t$', fontsize='x-large')
ylabel(r'$\frac{x_3}{\delta}$', fontsize='xx-large', rotation=0)
axis(xmin=0, xmax=10, ymin=-4, ymax=+4)
grid()

To compute the reaction $f_4$

  • we compute the displacements $x=\Psi q$
  • we compute the nodal forces $f_s=Kx$
  • we use the linear relation between $f_4$ and the $f_S$
  • we display the symbolic expression of $f_4$
In [17]:
x = Psi*Matrix(q123)
fs = K*x
R4 = R4_from_fs*fs ; R4 = R4[0]
display(Math('f_4 = k\delta\;(%s)'%(latex(R4.evalf(n=5)),)))
$$f_4 = k\delta\;(- 1.405 \operatorname{sin}\left(0.70245152 \omega_{0} t\right) + 52.01 \operatorname{sin}\left(2.7227407 \omega_{0} t\right) - 51.623 \operatorname{sin}\left(2.75 \omega_{0} t\right) + 0.29727 \operatorname{sin}\left(4.5102161 \omega_{0} t\right))$$

As before, the sybolic expression must be transformed in a numerical function

In [18]:
R4t = lambdify(t, R4.subs(w0,1), modules='numpy')

In this case the resonant behaviour of the system is particularly apparent.

In [19]:
plot(time,R4t(time))
xlabel(r'$\omega_0t$', fontsize='x-large')
ylabel(r'$f_4/(k\delta)$', fontsize='x-large')
grid()

A partial verification

As a partial verification the numerical response is computed using the linear acceleration algorithm with an adimensional time step $h=0.002$.

In [20]:
k = mat(((912,-684,-1482),(-684,2064,864),(-1482,864,2928)))/209.
m = mat(eye(3)*1.0)
r = mat('1;0;1')

h = 0.002 ; dur = 10.0

w = 2.75 ; w2 = w*w

def load(t):
    return -r*float(-w2*sin(w*t))

A = m*3 ; V = m*6/h
K = k + m*6/h/h
KI = K.I
MI = m # m being the unit matrix

t = 0
x = mat('0;0;0')
v = mat('0;0;0')
p = load(t)
cont = []

while t<dur+h/2:
    a = p-k.dot(x)
    cont.append([t, x[2,0]])
    t = t+h
    dp = load(t)-p
    dp_star = dp+A*a+V*v
    dx = KI*dp_star
    dv = (dx/h-v)*3-a*h/2
    x, v, p = x+dx, v+dv, p+dp

cont = array(cont)

Our last task will be to plot, with a continuous line, the numerically computed response and then superimpose on the plot dots obtained by sampling the analytical response function,

In [21]:
time = linspace(0,10,101)
plot(cont[:,0], cont[:,1], '-k', label='Numerical Response',
     linewidth=0.5)
plot(time,      x3t(time), 'ok', label='Analytical Response',
     alpha=0.5, markersize=2)
plot((0,10),(0,0),'-', color='lightgray', linewidth=1)
axis(xmin=0, xmax=10, ymin=-4, ymax=4)

title('Comparison of Numerical and Analytical Solutions')
xlabel(r'$\omega_0 t$', fontsize='x-large')
ylabel(r'$\frac{x_3}{\delta}$', fontsize='xx-large', rotation=0)

legend(fancybox=True, shadow=True, numpoints=5)
setp(gca().get_legend().get_texts(), fontsize='small')
setp(gca().get_legend().get_frame(), color='#FFF040')
setp(gca().get_legend().get_frame(), edgecolor='black')
grid()

The agreement between the two solutions is indeed very good.

In [21]: