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.


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

# Rayleigh approximations


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

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


```
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


```
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$.


```
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


```
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.


```
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...)


```
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)
```

Here we define the stiffness and the mass matrix


```
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.


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

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


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

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$.


```
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`.


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

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


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

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$.


```
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))
```

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


```
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))
```

Finally, the numerical values


```
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}$$


09529901561d252eda6ac2927a8d568bd92f8e34
ipynb
Jupyter Notebook
nb/05_Rayleigh.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
nb/05_Rayleigh.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
nb/05_Rayleigh.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
2019-06-23T12:32:39.000Z
2021-08-15T18:33:55.000Z
Qwen/Qwen-72B
1. YES
2. YES