Here we import everything from sympy, load the ability to pretty-print all sympy expressions and also the ability to display external images and LaTeX expressions.



In [1]:

    
from sympy import *
%load_ext sympyprinting
from IPython.core.display import Image, Math

Rayleigh approximations



In [2]:

    
Image(filename='figures/model.png')









    Out[2]:

We need the relative rotation across each one of the 5 hinges in the top part of the figure, so we introduce two fictitious bar so that we can compute the relative rotation as the difference of the rotations of the adjacent bars for all the hinges. Having introduced the two fictitious bar, to compute their rotations we have to take into account four additional degrees of freedom, as you can see in the bottom part of the figure.

We start defining the symbols that we will use



In [3]:

    
xm, x0, x1, x2, x3, x4, x5 = symbols('x_-1 x_0 x_1 x_2 x_3 x_4 x_5')
K, k, L , m, W, w, Zo, phi, theta = symbols('K k L m Omega omega Z_0 phi theta')

We define a vector of displacements



In [4]:

    
disp = (xm, x0, x1, x2, x3, x4, x5)
display(disp)









    





$$\begin{pmatrix}x_{-1}, & x_{0}, & x_{1}, & x_{2}, & x_{3}, & x_{4}, & x_{5}\end{pmatrix}$$

We define a function that takes the difference between two adjacent menbers of a sequence and optionally scales the difference by a factor, note that if the original sequence lenght is $N$, the sequence of differences that's returned by our function has lenght $N-1$.



In [5]:

    
def differences(seq, factor=1):
    return [(b-a)*factor for a, b in zip(seq[:-1],seq[1:])]

Now we use our differences to compute first the rotation of each of the six bars (the scaling factor being $\displaystyle\frac1L$) and later the relative rotations across each of the five hinges



In [36]:

    
rot = differences(disp, factor=1/L)
print "The rotations of the rods"
display(Math(r',\quad'.join([r'\phi_{%d} = %s'%(i+0,latex(r.simplify())) for i, r in enumerate(rot[:3])])+','))
display(Math(r',\quad'.join([r'\phi_{%d} = %s'%(i+3,latex(r.simplify())) for i, r in enumerate(rot[3:])])+'.'))

rrot = differences(rot)
print "\nThe relative rotations across the hinges:"
display(Math(r',\quad'.join([r'\theta_{%d} = %s'%(i+0,latex(r.simplify())) for i, r in enumerate(rrot[:3])])+','))
display(Math(r',\quad'.join([r'\theta_{%d} = %s'%(i+3,latex(r.simplify())) for i, r in enumerate(rrot[3:])])+'.'))
#dummy = [display(Math(r"\theta_%d = {%s}"%(i,latex(r.simplify())))) for i, r in enumerate(rrot)]









    



The rotations of the rods






    





$$\phi_{0} = \frac{- x_{-1} + x_{0}}{L},\quad\phi_{1} = \frac{- x_{0} + x_{1}}{L},\quad\phi_{2} = \frac{- x_{1} + x_{2}}{L},$$






    





$$\phi_{3} = \frac{- x_{2} + x_{3}}{L},\quad\phi_{4} = \frac{- x_{3} + x_{4}}{L},\quad\phi_{5} = \frac{- x_{4} + x_{5}}{L}.$$






    



The relative rotations across the hinges:






    





$$\theta_{0} = \frac{x_{-1} - 2 x_{0} + x_{1}}{L},\quad\theta_{1} = \frac{x_{0} - 2 x_{1} + x_{2}}{L},\quad\theta_{2} = \frac{x_{1} - 2 x_{2} + x_{3}}{L},$$






    





$$\theta_{3} = \frac{x_{2} - 2 x_{3} + x_{4}}{L},\quad\theta_{4} = \frac{x_{3} - 2 x_{4} + x_{5}}{L}.$$

Now, for each element of rrot, we substitute 0 for each of the fictitious displacements, to obtain the relative rotations across the springs in terms of the 3 free coordinates of our problem.



In [28]:

    
rrot = [rr.subs({xm:0,x0:0,x4:0,x5:0}) for rr in rrot]
print "The relative rotations across the hinges after substitution of 0 for fictitious degrees of freedom:"
display(Math(r',\quad'.join([r'\theta_{%d} = %s'%(i,latex(r.simplify())) for i, r in enumerate(rrot)])+'.'))









    



The relative rotations across the hinges after substitution of 0 for fictitious degrees of freedom:






    





$$\theta_{0} = \frac{x_{1}}{L},\quad\theta_{1} = \frac{- 2 x_{1} + x_{2}}{L},\quad\theta_{2} = \frac{x_{1} - 2 x_{2} + x_{3}}{L},\quad\theta_{3} = \frac{x_{2} - 2 x_{3}}{L},\quad\theta_{4} = \frac{x_{3}}{L}.$$

Next, we can compute the (double of the) strain energy in terms of the free coordinates (note that I substituted k for K/L**2, just to have the expression for $\omega^2$ look more familiar...)



In [8]:

    
V2 = sum([stif*rr**2 for stif, rr in zip((2*K, K, K, K, 2*K), rrot)]).expand()
V2 = V2.subs(K,k*L*L)
display(V2)









    





$$7 k x_{1}^{2} - 8 k x_{1} x_{2} + 2 k x_{1} x_{3} + 6 k x_{2}^{2} - 8 k x_{2} x_{3} + 7 k x_{3}^{2}$$

Here we define the stiffness and the mass matrix



In [9]:

    
k_mat=k*Matrix(((7, -4, 1), (-4, 6, -4), (1, -4, 7)))
m_mat = m*Matrix(((1,0,0),(0,1,0),(0,0,1)))
display(k_mat) ; display(m_mat)









    





$$\left(\begin{smallmatrix}7 k & - 4 k & k\\- 4 k & 6 k & - 4 k\\k & - 4 k & 7 k\end{smallmatrix}\right)$$






    





$$\left(\begin{smallmatrix}m & 0 & 0\\0 & m & 0\\0 & 0 & m\end{smallmatrix}\right)$$

In the following we will use many times the duoble matrix multiplication, so we define a shorthand function.



In [10]:

    
def dmp(M,v):
    return (v.transpose()*M*v)[0]

Our first use of dmp is for verifying the correctness of the stiffness matrix



In [11]:

    
V2test = dmp(k_mat,Matrix((x1,x2,x3))).expand()
display(V2test)









    





$$7 k x_{1}^{2} - 8 k x_{1} x_{2} + 2 k x_{1} x_{3} + 6 k x_{2}^{2} - 8 k x_{2} x_{3} + 7 k x_{3}^{2}$$

Having verified the correctness of the stiffness matrix, we define a vector of maximum displacements and a vector of maximum velocities in terms of the unknown $\omega$.



In [12]:

    
d_max = Matrix((1,1,1))*Zo
v_max = w*d_max
display(d_max)
display(v_max)









    





$$\left(\begin{smallmatrix}Z_{0}\\Z_{0}\\Z_{0}\end{smallmatrix}\right)$$






    





$$\left(\begin{smallmatrix}Z_{0} \omega\\Z_{0} \omega\\Z_{0} \omega\end{smallmatrix}\right)$$

The (double of the) strain and kinetic energies are now computed, using dmp.



In [13]:

    
V2 = dmp(k_mat,d_max)
T2 = dmp(m_mat,v_max)
display(Eq(V2, T2))









    





$$6 Z_{0}^{2} k = 3 Z_{0}^{2} m \omega^{2}$$

To solve correctly the equation above we have to make the substitution $\omega^2 = \Omega$, so that the solver is not confused by the unknown being squared



In [14]:

    
R00 = solve(Eq(V2,T2.subs(w,sqrt(W))),W)[0]
display(Eq(w*w, R00))









    





$$\omega^{2} = 2 \frac{k}{m}$$

We compute the maximum values of the forces of inertia multiplying by W${}=\Omega=\omega^2$, then a better approximation of the strain energy , that equated to the previuous value of the kinetic energy gives the new approximation to $\omega^2$.



In [15]:

    
f_in = -W*m_mat*d_max
V2 = dmp(k_mat.inv(),f_in)
R01 = solve(Eq(V2,T2.subs(w,sqrt(W))),W)[0]
display(Eq(w*w, R01))









    





$$\omega^{2} = \frac{4}{3} \frac{k}{m}$$

Computing a better approximation of the kinetic energy gives the third approximation,



In [16]:

    
v_max = -w*k_mat.inv()*f_in
T2 = dmp(m_mat,v_max)
R11 = solve(Eq(V2,T2.subs(w,sqrt(W))),W)[0]
display(Eq(w*w, R11))









    





$$\omega^{2} = \frac{24}{19} \frac{k}{m}$$

Finally, the numerical values



In [17]:

    
display([R.evalf() for R in (R00,R01,R11)])









    





$$\begin{bmatrix}2.0 \frac{k}{m}, & 1.33333333333333 \frac{k}{m}, & 1.26315789473684 \frac{k}{m}\end{bmatrix}$$