Generalized coordinates are \(y_{3},y_{2},y_{1}\). Kinetic energy is \(T=\frac {1}{2}m_{3}\left (y_{3}^{\prime }\right ) ^{2}+\frac {1}{2}m_{2}\left ( y_{2}^{\prime }\right ) ^{2}+\frac {1}{2}m_{1}\left (y_{1}^{\prime }\right ) ^{2}.\) Potential energy due to springs is \(V_{spring}=\frac {1}{2}k_{3}y_{3}^{2}+\frac {1}{2}k_{2}\left (y_{2}-y_{3}\right ) ^{2}+\frac {1}{2}k_{1}\left (y_{1}-y_{2}\right ) ^{2}\). Therefore\begin {align*} V_{spring} & =\frac {1}{2}k_{3}y_{3}^{2}+\frac {1}{2}k_{2}\left (y_{2}^{2}+y_{3}^{2}-2y_{2}y_{3}\right ) +\frac {1}{2}k_{1}\left (y_{1}^{2}+y_{2}^{2}-2y_{1}y_{2}\right ) \\ & =y_{3}^{2}\left (\frac {1}{2}k_{3}+\frac {1}{2}k_{2}\right ) +y_{2}^{2}\left (\frac {1}{2}k_{2}+\frac {1}{2}k_{1}\right ) +y_{1}^{2}\left ( \frac {1}{2}k_{1}\right ) +y_{1}y_{2}\left (-k_{1}\right ) +y_{1}y_{3}\left ( 0\right ) +y_{2}y_{3}\left (-k_{2}\right ) \end {align*}
The EOM is\[\begin {bmatrix} m_{1} & 0 & 0\\ 0 & m_{2} & 0\\ 0 & 0 & m_{3}\end {bmatrix}\begin {bmatrix} y_{1}^{\prime \prime }\\ y_{2}^{\prime \prime }\\ y_{3}^{\prime \prime }\end {bmatrix} +\begin {bmatrix} k_{1} & -k_{1} & 0\\ -k_{1} & k_{2}+k_{1} & -k_{2}\\ 0 & -k_{2} & k_{3}+k_{2}\end {bmatrix}\begin {bmatrix} y_{1}\\ y_{2}\\ y_{3}\end {bmatrix} =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \]
Following values are for mass (units in kg) \(m_{1}=100,m_{2}=200,m_{3}=300\). Following values are for spring constants (units in N/m) \(k_{1}=40^{2}\left ( 100\right ) ,k_{2}=50^{2}\left (200\right ) ,k_{3}=60^{2}\left (300\right ) .\) EOM becomes\[\begin {bmatrix} 100 & 0 & 0\\ 0 & 200 & 0\\ 0 & 0 & 300 \end {bmatrix}\begin {bmatrix} y_{1}^{\prime \prime }\\ y_{2}^{\prime \prime }\\ y_{3}^{\prime \prime }\end {bmatrix} +\begin {bmatrix} 160000 & -160000 & 0\\ -160000 & 660000 & -500000\\ 0 & -500000 & 1580000 \end {bmatrix}\begin {bmatrix} y_{1}\\ y_{2}\\ y_{3}\end {bmatrix} =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \]
Characteristic equation is\begin {align*} \det \left (\left [ K\right ] -\omega ^{2}\left [ M\right ] \right ) & =0\\ \det \left ( \begin {bmatrix} 160000 & -160000 & 0\\ -160000 & 660000 & -500000\\ 0 & -500000 & 1580000 \end {bmatrix} -\omega ^{2}\begin {bmatrix} 100 & 0 & 0\\ 0 & 200 & 0\\ 0 & 0 & 300 \end {bmatrix} \right ) & =0\\ \det \begin {bmatrix} 160000-100\omega ^{2} & -160000 & 0\\ -160000 & 660000-200\omega ^{2} & -500000\\ 0 & -500000 & 1580000-300\omega ^{2}\end {bmatrix} & =0\\ -6\times 10^{6}\omega ^{6}+6.1\times 10^{10}\omega ^{4}-1.54\times 10^{14}\omega ^{2}+8.64\times 10^{16} & =0 \end {align*}
Positive roots of the above polynomial are the natural frequencies (units in rad/sec). \begin {align*} \omega _{1} & =28.1\\ \omega _{2} & =52.6\\ \omega _{3} & =81.3 \end {align*}
To obtain mode shapes, the eigenvector associated with each eigenvalue is found. Starting with \(\omega _{1}=28.1\)\[ \left ( \begin {bmatrix} 160000 & -160000 & 0\\ -160000 & 660000 & -500000\\ 0 & -500000 & 1580000 \end {bmatrix} -28.1^{2}\begin {bmatrix} 100 & 0 & 0\\ 0 & 200 & 0\\ 0 & 0 & 300 \end {bmatrix} \right ) \begin {bmatrix} 1\\ \varphi _{21}\\ \varphi _{31}\end {bmatrix} =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \]
Hence\begin {align*} \begin {bmatrix} 8.1\times 10^{4} & -160000 & 0\\ -160000 & 5.02\times 10^{5} & -500000\\ 0 & -500000 & 1.34\times 10^{6}\end {bmatrix}\begin {bmatrix} 1\\ \varphi _{21}\\ \varphi _{31}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \\\begin {bmatrix} 8.1\times 10^{4}-1.6\times 10^{5}\varphi _{21}\\ 5.02\times 10^{5}\varphi _{21}-5.0\times 10^{5}\varphi _{31}-1.6\times 10^{5}\\ 1.34\times 10^{6}\varphi _{31}-5\times 10^{5}\varphi _{21}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \end {align*}
Solving gives \(\varphi _{21}=0.506\,\) and \(\varphi _{31}=0.188.\) First eigenvector is \[ \varphi _{1}=\begin {bmatrix} 1\\ 0.506\\ 0.188 \end {bmatrix} \]
For \(\omega _{2}=52.6\),\[ \left ( \begin {bmatrix} 160000 & -160000 & 0\\ -160000 & 660000 & -500000\\ 0 & -500000 & 1580000 \end {bmatrix} -52.6^{2}\begin {bmatrix} 100 & 0 & 0\\ 0 & 200 & 0\\ 0 & 0 & 300 \end {bmatrix} \right ) \begin {bmatrix} 1\\ \varphi _{22}\\ \varphi _{32}\end {bmatrix} =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \]
Hence\begin {align*} \begin {bmatrix} -1.17\times 10^{5} & -1.6\times 10^{5} & 0\\ -1.6\times 10^{5} & 1.07\times 10^{5} & -5.0\times 10^{5}\\ 0 & -5.0\times 10^{5} & 7.50\times 10^{5}\end {bmatrix}\begin {bmatrix} 1\\ \varphi _{22}\\ \varphi _{32}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \\\begin {bmatrix} -1.6\times 10^{5}\varphi _{22}-1.17\times 10^{5}\\ 1.07\times 10^{5}\varphi _{22}-5.0\times 10^{5}\varphi _{32}-1.6\times 10^{5}\\ 7.5\times 10^{5}\varphi _{32}-5.0\times 10^{5}\varphi _{22}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \end {align*}
Solving gives \(\varphi _{22}=-0.731\,\) and \(\varphi _{32}=-0.476\,.\)Second eigenvector is \[ \varphi _{2}=\begin {bmatrix} 1\\ -0.731\,\\ -0.476 \end {bmatrix} \]
For \(\omega _{3}=81.3\)\[ \left ( \begin {bmatrix} 160000 & -160000 & 0\\ -160000 & 660000 & -500000\\ 0 & -500000 & 1580000 \end {bmatrix} -81.3^{2}\begin {bmatrix} 100 & 0 & 0\\ 0 & 200 & 0\\ 0 & 0 & 300 \end {bmatrix} \right ) \begin {bmatrix} 1\\ \varphi _{23}\\ \varphi _{33}\end {bmatrix} =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \]
Hence\begin {align*} \begin {bmatrix} -5.\,01\times 10^{5} & -1.6\times 10^{5} & 0\\ -1.6\times 10^{5} & -6.62\times 10^{5} & -5.0\times 10^{5}\\ 0 & -5.0\times 10^{5} & -4.03\times 10^{5}\end {bmatrix}\begin {bmatrix} 1\\ \varphi _{23}\\ \varphi _{33}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \\\begin {bmatrix} -1.6\times 10^{5}\varphi _{23}-5.01\times 10^{5}\\ -6.62\times 10^{5}\varphi _{23}-5.0\times 10^{5}\varphi _{33}-1.6\times 10^{5}\\ -5.0\times 10^{5}\varphi _{23}-4.03\times 10^{5}\varphi _{33}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \end {align*}
Solving gives \(\varphi _{23}=-3.13\,\) and \(\varphi _{32}=\varphi _{33}=3.82\,.\) Third eigenvector is \[ \varphi _{3}=\begin {bmatrix} 1\\ -3.13\,\\ 3.82 \end {bmatrix} \]
Eigenvectors are mass normalized. Mass normalization factors \(\mu _{i}\) are found for each eigenvector\begin {align*} \mu _{1} & =\varphi _{1}^{T}\left [ M\right ] \varphi _{1}\\ & =\begin {bmatrix} 1\\ 0.506\\ 0.188 \end {bmatrix} ^{T}\begin {bmatrix} 100 & 0 & 0\\ 0 & 200 & 0\\ 0 & 0 & 300 \end {bmatrix}\begin {bmatrix} 1\\ 0.506\\ 0.188 \end {bmatrix} =162. \end {align*}
and\begin {align*} \mu _{2} & =\varphi _{2}^{T}\left [ M\right ] \varphi _{2}\\ & =\begin {bmatrix} 1\\ -0.731\,\\ -0.476 \end {bmatrix} ^{T}\begin {bmatrix} 100 & 0 & 0\\ 0 & 200 & 0\\ 0 & 0 & 300 \end {bmatrix}\begin {bmatrix} 1\\ -0.731\,\\ -0.476 \end {bmatrix} =275. \end {align*}
and\begin {align*} \mu _{3} & =\varphi _{3}^{T}\left [ M\right ] \varphi _{3}\\ & =\begin {bmatrix} 1\\ -3.13\,\\ 3.82 \end {bmatrix} ^{T}\begin {bmatrix} 100 & 0 & 0\\ 0 & 200 & 0\\ 0 & 0 & 300 \end {bmatrix}\begin {bmatrix} 1\\ -3.13\,\\ 3.82 \end {bmatrix} =6.44\times 10^{3} \end {align*}
Normalized eigenvectors are \begin {align*} \Phi _{1} & =\frac {\varphi _{1}}{\sqrt {\mu _{1}}}=\frac {1}{\sqrt {162}}\begin {bmatrix} 1\\ 0.506\\ 0.188 \end {bmatrix} =\begin {bmatrix} 7.86\times 10^{-2}\\ 3.98\times 10^{-2}\\ 1.48\times 10^{-2}\end {bmatrix} \\ \Phi _{2} & =\frac {\varphi _{2}}{\sqrt {\mu _{2}}}=\frac {1}{\sqrt {275.}}\begin {bmatrix} 1\\ -0.731\,\\ -0.476 \end {bmatrix} =\begin {bmatrix} 6.03\times 10^{-2}\\ -4.41\times 10^{-2}\\ -2.87\times 10^{-2}\end {bmatrix} \\ \Phi _{3} & =\frac {\varphi _{3}}{\sqrt {\mu _{3}}}=\frac {1}{\sqrt {6.44\times 10^{3}}}\begin {bmatrix} 1\\ -3.13\,\\ 3.82 \end {bmatrix} =\begin {bmatrix} 1.25\times 10^{-2}\\ -0.039\,\\ 4.76\times 10^{-2}\end {bmatrix} \end {align*}
Verification of the above result follows
EDU>> k=[160000 -160000 0;-160000 660000 -500000;0 -500000 1580000]; EDU>> M=[100 0 0;0 200 0;0 0 300]; EDU>> [eigV,lam]=eig(k,M) eigV = 0.0786 0.0606 0.0124 0.0398 -0.0437 -0.0389 0.0148 -0.0289 0.0477 lam = 1.0e+03 * 0.7897 0 0 0 2.7528 0 0 0 6.6242 EDU>> sqrt(diag(lam)) ans = 28.1013 52.4674 81.3889
Initial conditions are \(\theta _{i}\relax (0) =0\) for \(i=1,2,3\) and \(\theta _{1}^{\prime }\relax (0) =\) \(\theta _{3}^{\prime }\relax (0) =0\) but \(\theta _{2}^{\prime }\relax (0) =2\) rad/sec.
The generalized coordinates are shown above. kinetic energy is\[ T=\frac {1}{2}I\left (\theta _{1}^{\prime }\right ) ^{2}+\frac {1}{2}I\left ( \theta _{2}^{\prime }\right ) ^{2}+\frac {1}{2}I\left (\theta _{3}^{\prime }\right ) ^{2}\]
where \(I=\frac {1}{3}mL^{2}.\) Mas matrix becomes\[ \left [ M\right ] =\frac {1}{3}mL^{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} \theta _{1}^{\prime \prime }\\ \theta _{2}^{\prime \prime }\\ \theta _{3}^{\prime \prime }\end {bmatrix} \]
Spring potential energy is\begin {align*} V_{spring} & =\frac {1}{2}k\left (L\theta _{2}-L\theta _{1}\right ) ^{2}+\frac {1}{2}k\left (L\theta _{3}-L\theta _{2}\right ) ^{2}\\ & =\frac {1}{2}kL^{2}\left (\theta _{2}^{2}+\theta _{1}^{2}-2\theta _{1}\theta _{2}\right ) +\frac {1}{2}kL^{2}\left (\theta _{3}^{2}+\theta _{2}^{2}-2\theta _{2}\theta _{3}\right ) \\ & =\theta _{1}^{2}\left (\frac {1}{2}kL^{2}\right ) +\theta _{2}^{2}\left ( \frac {1}{2}kL^{2}+\frac {1}{2}kL^{2}\right ) +\theta _{3}^{2}\left (\frac {1}{2}kL^{2}\right ) +\theta _{1}\theta _{2}\left (-kL^{2}\right ) +\theta _{1}\theta _{3}\relax (0) +\theta _{2}\theta _{3}\left (-kL^{2}\right ) \end {align*}
Hence stiffness matrix due to spring is\[ \left [ K\right ] _{spring}=kL^{2}\begin {bmatrix} 1 & -1 & 0\\ -1 & 2 & -1\\ 0 & -1 & 1 \end {bmatrix} \]
Assume the zero PE for gravity is taken as the top of the bar. Stiffness due to gravity is\[ V_{g}=-mg\frac {L}{2}\left (\cos \theta _{1}+\cos \theta _{2}+\cos \theta _{3}\right ) \]
\(V_{11}=\frac {\partial ^{2}V_{g}}{\partial \theta _{1}^{2}}=mg\frac {L}{2}\left ( \cos \theta _{1}\right ) .\) Evaluate this at static position \(\theta _{1}=0,\)hence \(V_{11}=m\frac {L}{2}\). Similarly, \(V_{22}=V_{33}=m\frac {L}{2}\,\). All other terms are zero.
Hence stiffness matrix due to gravity is\[ \left [ K\right ] _{g}=mg\frac {L}{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix} \]
Therefore, complete stiffness matrix is\[ kL^{2}\begin {bmatrix} 1 & -1 & 0\\ -1 & 2 & -1\\ 0 & -1 & 1 \end {bmatrix} +mg\frac {L}{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix} \]
There are no generalized forces. Hence EOM is\[ \frac {1}{3}mL^{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} \theta _{1}^{\prime \prime }\\ \theta _{2}^{\prime \prime }\\ \theta _{3}^{\prime \prime }\end {bmatrix} +\left (kL^{2}\begin {bmatrix} 1 & -1 & 0\\ -1 & 2 & -1\\ 0 & -1 & 1 \end {bmatrix} +mg\frac {L}{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix} \right ) \begin {bmatrix} \theta _{1}\\ \theta _{2}\\ \theta _{3}\end {bmatrix} =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \]
For case \(k=0.05\frac {mg}{L}\), Hence for \(\sigma =0.05\) then \(k=\sigma \frac {mg}{L}.\) EOM becomes
\begin {align*} \frac {1}{3}mL^{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} \theta _{1}^{\prime \prime }\\ \theta _{2}^{\prime \prime }\\ \theta _{3}^{\prime \prime }\end {bmatrix} +\left (\sigma mgL\begin {bmatrix} 1 & -1 & 0\\ -1 & 2 & -1\\ 0 & -1 & 1 \end {bmatrix} +mg\frac {L}{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix} \right ) \begin {bmatrix} \theta _{1}\\ \theta _{2}\\ \theta _{3}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \\ \frac {1}{3}mL^{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} \theta _{1}^{\prime \prime }\\ \theta _{2}^{\prime \prime }\\ \theta _{3}^{\prime \prime }\end {bmatrix} +mgL\begin {bmatrix} \frac {1}{2}+\sigma & -\sigma & 0\\ -\sigma & \frac {1}{2}+2\sigma & -\sigma \\ 0 & -\sigma & \frac {1}{2}+\sigma \end {bmatrix} \begin {bmatrix} \theta _{1}\\ \theta _{2}\\ \theta _{3}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \\\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} \theta _{1}^{\prime \prime }\\ \theta _{2}^{\prime \prime }\\ \theta _{3}^{\prime \prime }\end {bmatrix} +\frac {3g}{L}\begin {bmatrix} \frac {1}{2}+\sigma & -\sigma & 0\\ -\sigma & \frac {1}{2}+2\sigma & -\sigma \\ 0 & -\sigma & \frac {1}{2}+\sigma \end {bmatrix} \begin {bmatrix} \theta _{1}\\ \theta _{2}\\ \theta _{3}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \end {align*}
Let \(L=1,g=10\). The above becomes\[\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} \theta _{1}^{\prime \prime }\\ \theta _{2}^{\prime \prime }\\ \theta _{3}^{\prime \prime }\end {bmatrix} +\begin {bmatrix} 15+30\sigma & -30\sigma & 0\\ -30\sigma & 15+60\sigma & -30\sigma \\ 0 & -30\sigma & 15+30\sigma \end {bmatrix}\begin {bmatrix} \theta _{1}\\ \theta _{2}\\ \theta _{3}\end {bmatrix} =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \]
Natural frequencies of the system are found by solving the eigenvalue problem.\[ \det \left ( \begin {bmatrix} 15+30\sigma & -30\sigma & 0\\ -30\sigma & 15+60\sigma & -30\sigma \\ 0 & -30\sigma & 15+30\sigma \end {bmatrix} -\omega ^{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix} \right ) =0 \]
Substituting \(\sigma =0.05\) gives\begin {align*} \det \left ( \begin {bmatrix} 16.5 & -1.5 & 0\\ -1.5 & 18.0 & -1.5\\ 0 & -1.5 & 16.5 \end {bmatrix} -\omega ^{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix} \right ) & =0\\ \det \begin {bmatrix} 16.5-\omega ^{2} & -1.5 & 0\\ -1.5 & 18-\omega ^{2} & -1.5\\ 0 & -1.5 & 16.5-\omega ^{2}\end {bmatrix} & =0\\ -\omega ^{6}+51\omega ^{4}-861.75\omega ^{2}+4826.3 & =0 \end {align*}
Positive roots of this polynomial are \(\omega =3.87,\omega =4.062,\omega =4.416\).
Associated eigenvectors are found by solving for \(\varphi _{i}\) in \(\left ( \left [ K\right ] -\omega ^{2}\left [ M\right ] \right ) \varphi _{i}=0\) for each eigenvalue \(\omega _{i}\).
For \(\omega _{1}=3.87\) \begin {align*} \begin {bmatrix} 16.5-\omega _{1}^{2} & -1.5 & 0\\ -1.5 & 18-\omega _{1}^{2} & -1.5\\ 0 & -1.5 & 16.5-\omega _{1}^{2}\end {bmatrix}\begin {bmatrix} 1\\ \varphi _{21}\\ \varphi _{31}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \\\begin {bmatrix} 16.5-3.87^{2} & -1.5 & 0\\ -1.5 & 18-3.87^{2} & -1.5\\ 0 & -1.5 & 16.5-3.87^{2}\end {bmatrix}\begin {bmatrix} 1\\ \varphi _{21}\\ \varphi _{31}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \\\begin {bmatrix} 1.523\, & -1.5 & 0\\ -1.5 & 3.023\, & -1.5\\ 0 & -1.5 & 1.523\, \end {bmatrix}\begin {bmatrix} 1\\ \varphi _{21}\\ \varphi _{31}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \\\begin {bmatrix} 1.523\,-1.5\varphi _{21}\\ 3.023\,\varphi _{21}-1.5\varphi _{31}-1.5\\ 1.523\,\varphi _{31}-1.5\varphi _{21}\end {bmatrix} & =\begin {bmatrix} 0\\ 0\\ 0 \end {bmatrix} \end {align*}
Solving gives \(\varphi _{21}=1.015\,3\) and \(\varphi _{31}=1.046\,2\). First eigenvector is\[ \varphi _{1}=\begin {bmatrix} 1\\ \varphi _{21}\\ \varphi _{31}\end {bmatrix} =\begin {bmatrix} 1\\ 1.0153\\ 1.0462 \end {bmatrix} \]
Similarly, second and the third eigenvectors are found.
Eigenvectors are mass normalized. First the mass normalization factors \(\mu _{i}\) are found for each eigenvector\begin {align*} \mu _{1} & =\varphi _{1}^{T}\left [ M\right ] \varphi _{1}\\ & =\begin {bmatrix} 1\\ 1.015\,3\\ 1.046\,2 \end {bmatrix} ^{T}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} 1\\ 1.0153\\ 1.0462 \end {bmatrix} =3.125\,4 \end {align*}
Normalized eigenvector is\[ \Phi _{1}=\frac {\varphi _{1}}{\sqrt {3.125\,4}}=\frac {1}{\sqrt {3.792}}\begin {bmatrix} 1\\ 1.015\,3\\ 1.046\,2 \end {bmatrix} =\begin {bmatrix} 0.51353\\ 0.52139\\ 0.53726 \end {bmatrix} \]
Verification of the above result (Matlab result is more accurate due to more accurate method used)
EDU>> k=[0.55 -0.05 0;-0.05 0.6 -0.05;0 -0.05 0.55]; EDU>> M=eye(3); EDU>> [eigV,lam]=eig(k,M) eigV = -0.5774 -0.7071 0.4082 -0.5774 -0.0000 -0.8165 -0.5774 0.7071 0.4082 EDU>> sqrt(diag(lam)) ans = 0.7071 0.7416 0.8062
Transformation matrix (based on Matlab more accurate result) is \[ \Phi =\left [ \Phi _{1}\Phi _{2}\Phi _{3}\right ] =\begin {bmatrix} -0.577\, & -0.7073 & 0.4082\\ -0.577\, & 0 & -0.8165\\ -0.577\, & 0.706\,9 & 0.4082 \end {bmatrix} \] Mapping from physical coordinates \(\theta \) to modal coordinates \(\eta \) is\[ \mathbf {\theta }=\left [ \Phi \right ] \mathbf {\eta }\]
Bold face is used to indicate a column vector. EOM’s are written in modal coordinates resulting in\begin {align*} \begin {bmatrix} 1\, & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} \eta _{1}^{\prime \prime }\,\\ \eta _{2}^{\prime \prime }\\ \eta _{3}^{\prime \prime }\end {bmatrix} +\begin {bmatrix} \omega _{1}^{2}\, & 0 & 0\\ 0 & \omega _{2}^{2} & 0\\ 0 & 0 & \omega _{2}^{2}\end {bmatrix}\begin {bmatrix} \eta _{1}\,\\ \eta _{2}\\ \eta _{3}\end {bmatrix} & =\begin {bmatrix} 0\,\\ 0\\ 0 \end {bmatrix} \\\begin {bmatrix} 1\, & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} \eta _{1}^{\prime \prime }\,\\ \eta _{2}^{\prime \prime }\\ \eta _{3}^{\prime \prime }\end {bmatrix} +\begin {bmatrix} 0.7071^{2}\, & 0 & 0\\ 0 & 0.7416^{2} & 0\\ 0 & 0 & 0.8062^{2}\end {bmatrix}\begin {bmatrix} \eta _{1}\,\\ \eta _{2}\\ \eta _{3}\end {bmatrix} & =\begin {bmatrix} 0\,\\ 0\\ 0 \end {bmatrix} \end {align*}
Initial conditions are transformed to modal coordinates using \(\mathbf {\eta }\relax (0) =\left [ \Phi \right ] ^{T}\left [ M\right ] \mathbf {\theta }\relax (0) \) and \(\mathbf {\eta }^{\prime }\relax (0) =\left [ \Phi \right ] ^{T}\left [ M\right ] \mathbf {\theta }^{\prime }\relax (0) \), since \(\mathbf {\theta }\relax (0) =0\) then \(\mathbf {\eta }\left ( 0\right ) =0\), however \(\mathbf {\theta }^{\prime }\relax (0) \) is not all zero, hence\begin {align*} \begin {bmatrix} \eta _{1}^{\prime }\relax (0) \,\\ \eta _{2}^{\prime }\relax (0) \\ \eta _{3}^{\prime }\relax (0) \end {bmatrix} & =\begin {bmatrix} -0.577\, & -0.7073 & 0.4082\\ -0.577\, & 0 & -0.8165\\ -0.577\, & 0.7069 & 0.4082 \end {bmatrix} ^{T}\begin {bmatrix} 1\, & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} 0\,\\ 2\\ 0 \end {bmatrix} \\ & =\begin {bmatrix} -1.154\\ 0\\ -1.633 \end {bmatrix} \end {align*}
Initial conditions in modal coordinates are found. The solution can be found. The solution to \(\eta ^{\prime \prime }+\omega ^{2}\eta =0\) with initial conditions \(\eta \relax (0) \) and \(\eta ^{\prime }\relax (0) \) is \(\eta \relax (t) =\eta \relax (0) \cos \omega t+\frac {\eta ^{\prime }\relax (0) }{\omega }\sin \omega t\). Therefore modal solutions are\begin {align*} \eta _{1}\relax (t) & =\frac {-1.154}{0.7071}\sin \left ( 0.7071t\right ) =-1.632\sin \left (0.707\,t\right ) \\ \eta _{2}\relax (t) & =0\\ \eta _{3}\relax (t) & =\frac {-1.633}{0.8062}\sin \left ( 0.8062t\right ) =-2.026\sin \left (0.8062t\right ) \end {align*}
Solution in the normal coordinates is\begin {align*} \begin {bmatrix} \theta _{1}\relax (t) \,\\ \theta _{2}\relax (t) \\ \theta _{3}\relax (t) \end {bmatrix} & =\begin {bmatrix} -0.577\, & -0.7073 & 0.4082\\ -0.577\, & 0 & -0.8165\\ -0.577\, & 0.706\,9 & 0.4082 \end {bmatrix}\begin {bmatrix} -1.632\sin \left (0.707\,t\right ) \,\\ 0\\ -2.026\sin \left (0.8062t\right ) \end {bmatrix} \\ & =\begin {bmatrix} 0.94166\sin \left (0.707\,t\right ) -0.82701\sin \left (0.8062t\right ) \\ 0.94166\sin \left (0.707\,t\right ) +1.6542\sin \left (0.8062t\right ) \\ 0.94166\sin \left (0.707\,t\right ) -0.82701\sin \left (0.8062t\right ) \end {bmatrix} \end {align*}
Using part (a), but with \(\sigma =2\) results in\[ \det \left ( \begin {bmatrix} 15+30\sigma & -30\sigma & 0\\ -30\sigma & 15+60\sigma & -30\sigma \\ 0 & -30\sigma & 15+30\sigma \end {bmatrix} -\omega ^{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix} \right ) =0 \]
\begin {align*} \det \left ( \begin {bmatrix} 15+30\relax (2) & -30\relax (2) & 0\\ -30\relax (2) & 15+60\relax (2) & -30\relax (2) \\ 0 & -30\relax (2) & 15+30\relax (2) \end {bmatrix} -\omega ^{2}\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix} \right ) & =0\\ \det \begin {bmatrix} 75.0-\omega ^{2} & -60.0 & 0\\ -60.0 & 135.0-\omega ^{2} & -60.0\\ 0 & -60.0 & 75.0-\omega ^{2}\end {bmatrix} & =0 \end {align*}
Similar steps as repeated as part (a) above. The final result are shown below using Matlab
EDU>> k=[75 -60 0;-60 135 -60;0 -60 75] EDU>> M=eye(3); [eigV,lam]=eig(k,M) eigV = -0.5774 -0.7071 0.4082 -0.5774 0.0000 -0.8165 -0.5774 0.7071 0.4082 EDU>> sqrt(diag(lam)) 3.8730 8.6603 13.9642
Transformation matrix is \(\Phi =\left [ \Phi _{1}\Phi _{2}\Phi _{3}\right ] =\begin {bmatrix} -0.577\, & -0.7071 & 0.4082\\ -0.577\, & 0 & -0.816\,5\\ -0.577\, & 0.7071 & 0.4082 \end {bmatrix} .\) Mapping from \(\theta \) to modal coordinates \(\eta \) is\[ \mathbf {\theta }=\left [ \Phi \right ] \mathbf {\eta }\]
Bold face is used to indicate a column vector. EOM’s are written in modal coordinates resulting in\begin {align*} \begin {bmatrix} 1\, & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} \eta _{1}^{\prime \prime }\,\\ \eta _{2}^{\prime \prime }\\ \eta _{3}^{\prime \prime }\end {bmatrix} +\begin {bmatrix} \omega _{1}^{2}\, & 0 & 0\\ 0 & \omega _{2}^{2} & 0\\ 0 & 0 & \omega _{2}^{2}\end {bmatrix}\begin {bmatrix} \eta _{1}\,\\ \eta _{2}\\ \eta _{3}\end {bmatrix} & =\begin {bmatrix} 0\,\\ 0\\ 0 \end {bmatrix} \\\begin {bmatrix} 1\, & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {bmatrix} \eta _{1}^{\prime \prime }\,\\ \eta _{2}^{\prime \prime }\\ \eta _{3}^{\prime \prime }\end {bmatrix} +\begin {bmatrix} 3.8730^{2}\, & 0 & 0\\ 0 & 8.6603^{2} & 0\\ 0 & 0 & 13.9642^{2}\end {bmatrix}\begin {bmatrix} \eta _{1}\,\\ \eta _{2}\\ \eta _{3}\end {bmatrix} & =\begin {bmatrix} 0\,\\ 0\\ 0 \end {bmatrix} \end {align*}
Initial conditions are transformed to modal coordinates using \(\mathbf {\eta }\relax (0) =\left [ \Phi \right ] ^{T}\left [ M\right ] \mathbf {x}\relax (0) \) and \(\mathbf {\eta }^{\prime }\relax (0) =\left [ \Phi \right ] ^{T}\left [ M\right ] \theta ^{\prime }\relax (0) \), since \(\theta \relax (0) =0\) then \(\mathbf {\eta }\relax (0) =0\), however \(\theta ^{\prime }\relax (0) \) is not all zero. Similar to part (a), initial conditions are found \[\begin {bmatrix} \eta _{1}^{\prime }\relax (0) \,\\ \eta _{2}^{\prime }\relax (0) \\ \eta _{3}^{\prime }\relax (0) \end {bmatrix} =\begin {bmatrix} -1.154\\ 0\\ -1.633 \end {bmatrix} \]
The solution to \(\eta ^{\prime \prime }+\lambda ^{2}\eta =0\) with initial conditions \(\eta \relax (0) \) and \(\eta ^{\prime }\relax (0) \) is given by \(\eta \relax (t) =\eta \relax (0) \cos \lambda t+\frac {\eta ^{\prime }\relax (0) }{\lambda }\sin \lambda t\). The solutions are\begin {align*} \eta _{1}\relax (t) & =\frac {-1.154}{3.8730}\sin \left (3.873t\right ) =-0.297\,96\sin \left (3.873\,t\right ) \\ \eta _{2}\relax (t) & =0\\ \eta _{3}\relax (t) & =\frac {-1.633}{13.9642}\sin \left ( 13.9642\right ) =-0.116\,94\sin \left (13.9642\right ) \end {align*}
Solution in the physical coordinates is\begin {align*} \begin {bmatrix} \theta _{1}\relax (t) \,\\ \theta _{2}\relax (t) \\ \theta _{3}\relax (t) \end {bmatrix} & =\begin {bmatrix} -0.577\, & -0.7071 & 0.4082\\ -0.577\, & 0 & -0.816\,5\\ -0.577\, & 0.7071 & 0.4082 \end {bmatrix}\begin {bmatrix} -0.297\,96\sin \left (3.8730\,t\right ) \,\\ 0\\ -0.116\,94\sin \left (13.9642t\right ) \end {bmatrix} \\ & =\begin {bmatrix} 0.171\,92\sin \left (3.873t\right ) -4.773\,5\times 10^{-2}\sin \left ( 13.964t\right ) \\ 9.548\,2\times 10^{-2}\sin \left (13.964t\right ) +0.171\,92\sin \left ( 3.873t\right ) \\ 0.171\,92\sin \left (3.873t\right ) -4.773\,5\times 10^{-2}\sin \left ( 13.964t\right ) \end {bmatrix} \end {align*}
Summary table
| \(k\) | frequencies | \(\left [ \Phi \right ] \) | solutions in \(\theta \) |
| \(0.05\frac {mg}{L}\) | \(\begin {Bmatrix} 0.7071\\ 0.7416\\ 0.8062 \end {Bmatrix} \) | \(\begin {bmatrix} -0.5774 & -0.7071 & 0.4082\\ -0.5774 & 0 & -0.8165\\ -0.5774 & 0.7071 & 0.4082 \end {bmatrix} \) | \(\begin {bmatrix} 0.94166\sin \left (0.707\,t\right ) -0.82701\sin \left (0.8062t\right ) \\ 0.94166\sin \left (0.707\,t\right ) +1.6542\sin \left (0.8062t\right ) \\ 0.94166\sin \left (0.707\,t\right ) -0.82701\sin \left (0.8062t\right ) \end {bmatrix} \) |
| \(2\frac {mg}{L}\) | \(\begin {Bmatrix} 3.8730\\ 8.6603\\ 13.9642 \end {Bmatrix} \) | \(\begin {bmatrix} -0.5774 & -0.7071 & 0.4082\\ -0.5774 & 0 & -0.8165\\ -0.5774 & 0.7071 & 0.4082 \end {bmatrix} \) | \(\begin {bmatrix} 0.17192\sin \left (3.873t\right ) -4.773\,5\times 10^{-2}\sin \left ( 13.964t\right ) \\ 9.5482\times 10^{-2}\sin \left (13.964t\right ) +0.171\,92\sin \left ( 3.873t\right ) \\ 0.17192\sin \left (3.873t\right ) -4.773\,5\times 10^{-2}\sin \left ( 13.964t\right ) \end {bmatrix} \) |
Even though the normalized natural frequencies are different, the shape functions are the same.
Plots of the solutions of \(\theta _{i}\relax (t) \) for both cases are made. For the case of \(k=0.05\frac {mg}{L}\)
For \(k=2\frac {mg}{L}\)
In addition, a small program is written to animate both the full solution and the modal solutions. The program to animate the full solution is at http://12000.org/my_courses/univ_wisconson_madison/spring_2013/EMA_545_Mechanical_Vibrations/HWs/HW10/HW10p2.m.txt while the program that animate the modal solution is number 112 at bottom of this page http://12000.org/my_notes/my_matlab_functions/index.htm
EOM is
\[\begin {bmatrix} 600 & 400 & 200\\ 400 & 1200 & 0\\ 200 & 0 & 800 \end {bmatrix}\begin {Bmatrix} x_{1}^{\prime \prime }\\ x_{2}^{\prime \prime }\\ x_{3}^{\prime \prime }\end {Bmatrix} +\begin {bmatrix} 500 & 300 & -400\\ 300 & 900 & 600\\ -400 & 600 & 1300 \end {bmatrix}\begin {Bmatrix} x_{1}^{\prime }\\ x_{2}^{\prime }\\ x_{3}^{\prime }\end {Bmatrix} +10^{3}\begin {bmatrix} 300 & 0 & -200\\ 0 & 500 & 300\\ -200 & 300 & 700 \end {bmatrix}\begin {Bmatrix} x_{1}\\ x_{2}\\ x_{3}\end {Bmatrix} =\begin {Bmatrix} 200\cos \left (16t\right ) \\ 0\\ 0 \end {Bmatrix} \]
Initial conditions are \(\begin {Bmatrix} x_{1}\relax (0) \\ x_{2}\relax (0) \\ x_{3}\relax (0) \end {Bmatrix} =\begin {Bmatrix} 0\\ 0\\ 0 \end {Bmatrix} \) and \(\begin {Bmatrix} x_{1}^{\prime }\relax (0) \\ x_{2}^{\prime }\relax (0) \\ x_{3}^{\prime }\relax (0) \end {Bmatrix} =\begin {Bmatrix} 0\\ 0\\ 0 \end {Bmatrix} \).
Solve the eigenvalue problem to determine the natural frequencies of the system\begin {align*} \det \left (\left [ K\right ] -\omega ^{2}\left [ M\right ] \right ) & =0\\ \det \left ( \begin {bmatrix} 300\times 10^{3} & 0 & -200\times 10^{3}\\ 0 & 500\times 10^{3} & 300\times 10^{3}\\ -200\times 10^{3} & 300\times 10^{3} & 700\times 10^{3}\end {bmatrix} -\omega ^{2}\begin {bmatrix} 600 & 400 & 200\\ 400 & 1200 & 0\\ 200 & 0 & 800 \end {bmatrix} \right ) & =0\\ -4.0\times 10^{8}\omega ^{6}+1.044\times 10^{12}\omega ^{4}-4.72\times 10^{14}\omega ^{2}+5.8\times 10^{16} & =0 \end {align*}
Positive roots are \(\left \{ \omega =15.052,\omega =17.562,\omega =45.552\right \} .\) For each natural frequency the corresponding eigenvector is found. A program is now used to compute these values.
EDU>> k = [300 0 -200;0 500 300;-200 300 700]*10^3; M = [600 400 200;400 1200 0;200 0 800]; C = [500 300 -400;300 900 600;-400 600 1300]; [PHI,lam] = eig(k,M); PHI lam = sqrt(diag(lam)) CC = PHI'*C*PHI; zeta1 = CC(1,1)/(2*lam(1)) zeta2 = CC(2,2)/(2*lam(2)) zeta3 = CC(3,3)/(2*lam(3)) PHI = -0.0216 0.0232 -0.0373 0.0203 0.0168 0.0201 -0.0220 0.0023 0.0302 lam = 15.0519 17.5624 45.5522 zeta1 = 0.0018 zeta2 = 0.0219 zeta3 = 0.0376
\(\left [ \Phi \right ] =\begin {bmatrix} -0.0216 & 0.0232 & -0.0373\\ 0.0203 & 0.0168 & 0.0201\\ -0.0220 & 0.0023 & 0.0302 \end {bmatrix} .\) In modal coordinates EOM is
\begin {align*} \begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {Bmatrix} \eta _{1}^{\prime \prime }\\ \eta _{2}^{\prime \prime }\\ \eta _{3}^{\prime \prime }\end {Bmatrix} +\left [ \Phi \right ] ^{T}\begin {bmatrix} 500 & 300 & -400\\ 300 & 900 & 600\\ -400 & 600 & 1300 \end {bmatrix} \left [ \Phi \right ] \begin {Bmatrix} \eta _{1}^{\prime }\\ \eta _{2}^{\prime }\\ \eta _{3}^{\prime }\end {Bmatrix} +\begin {bmatrix} \omega _{1}^{2} & 0 & 0\\ 0 & \omega _{2}^{2} & 0\\ 0 & 0 & \omega _{3}^{2}\end {bmatrix}\begin {Bmatrix} \eta _{1}\\ \eta _{2}\\ \eta _{3}\end {Bmatrix} & =\left [ \Phi \right ] ^{T}\begin {Bmatrix} 200\cos \left (16t\right ) \\ 0\\ 0 \end {Bmatrix} \\\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {Bmatrix} \eta _{1}^{\prime \prime }\\ \eta _{2}^{\prime \prime }\\ \eta _{3}^{\prime \prime }\end {Bmatrix} +\begin {bmatrix} 5.419\times 10^{-2} & 5.331\times 10^{-2} & -0.416\\ 5.331\times 10^{-2} & 0.768 & -3.52\times 10^{-4}\\ -0.4156 & -3.52\times 10^{-4} & 3.428 \end {bmatrix}\begin {Bmatrix} \eta _{1}^{\prime }\\ \eta _{2}^{\prime }\\ \eta _{3}^{\prime }\end {Bmatrix} +\begin {bmatrix} \omega _{1}^{2} & 0 & 0\\ 0 & \omega _{2}^{2} & 0\\ 0 & 0 & \omega _{3}^{2}\end {bmatrix}\begin {Bmatrix} \eta _{1}\\ \eta _{2}\\ \eta _{3}\end {Bmatrix} & =\left [ \Phi \right ] ^{T}\begin {Bmatrix} 200\cos \left (16t\right ) \\ 0\\ 0 \end {Bmatrix} \\\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {Bmatrix} \eta _{1}^{\prime \prime }\\ \eta _{2}^{\prime \prime }\\ \eta _{3}^{\prime \prime }\end {Bmatrix} +\begin {bmatrix} 2\zeta _{1}\omega _{1} & 0 & 0\\ 0 & 2\zeta _{2}\omega _{2} & 0\\ 0 & 0 & 2\zeta _{3}\omega _{3}\end {bmatrix}\begin {Bmatrix} \eta _{1}^{\prime }\\ \eta _{2}^{\prime }\\ \eta _{3}^{\prime }\end {Bmatrix} +\begin {bmatrix} \omega _{1}^{2} & 0 & 0\\ 0 & \omega _{2}^{2} & 0\\ 0 & 0 & \omega _{3}^{2}\end {bmatrix}\begin {Bmatrix} \eta _{1}\\ \eta _{2}\\ \eta _{3}\end {Bmatrix} & =\left [ \Phi \right ] ^{T}\begin {Bmatrix} 200\cos \left (16t\right ) \\ 0\\ 0 \end {Bmatrix} \end {align*}
In the above \(2\zeta _{1}\omega _{1}=0.0542\), \(2\zeta _{2}\omega _{2}=0.7676\) and \(2\zeta _{3}\omega _{3}=3.424\,7\). Hence \(\zeta _{1}=\frac {5.419\,3\times 10^{-2}}{2\left (15.0519\right ) }=0.0018\) and \(\zeta _{2}=\frac {0.76755}{2\left (17.5624\right ) }=0.0219\) and \(\zeta _{3}=\frac {3.424\,7}{2\left (45.5522\right ) }=0.0376\)
Final EOM in modal coordinates is\[\begin {bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end {bmatrix}\begin {Bmatrix} \eta _{1}^{\prime \prime }\\ \eta _{2}^{\prime \prime }\\ \eta _{3}^{\prime \prime }\end {Bmatrix} +\begin {bmatrix} 0.0542 & 0 & 0\\ 0 & 0.768 & 0\\ 0 & 0 & 3.425 \end {bmatrix}\begin {Bmatrix} \eta _{1}^{\prime }\\ \eta _{2}^{\prime }\\ \eta _{3}^{\prime }\end {Bmatrix} +\begin {bmatrix} 226.56 & 0 & 0\\ 0 & 308.44 & 0\\ 0 & 0 & 2075 \end {bmatrix}\begin {Bmatrix} \eta _{1}\\ \eta _{2}\\ \eta _{3}\end {Bmatrix} =\begin {Bmatrix} -4.32\cos \left (16.0t\right ) \\ 4.64\cos \left (16.0t\right ) \\ -7.46\cos \left (16.0t\right ) \end {Bmatrix} \]
EOM’s to solve are\begin {align*} \eta _{1}^{\prime \prime }+2\zeta _{1}\omega _{1}\eta _{1}^{\prime }+\omega _{1}^{2}\eta _{1} & =-4.32\cos \left (16.0t\right ) \\ \eta _{2}^{\prime \prime }+2\zeta _{2}\omega _{2}\eta _{2}^{\prime }+\omega _{2}^{2}\eta _{2} & =4.64\cos \left (16.0t\right ) \\ \eta _{3}^{\prime \prime }+2\zeta _{3}\omega _{3}\eta _{3}^{\prime }+\omega _{3}^{2}\eta _{3} & =-7.46\cos \left (16.0t\right ) \end {align*}
Initial conditions are zero. The solution in modal coordinates is given in appendix B for underdamped case. Complete solution for the case of underdamped is given in appendix B as\[ \eta \relax (t) =\frac {F_{0}}{\beta ^{2}+4\zeta ^{2}\omega ^{2}\varpi ^{2}}\left \{ \beta \cos \left (\varpi t\right ) +2\zeta \omega \varpi \sin \left ( \varpi t\right ) -e^{-\zeta \omega t}\left [ \beta \cos \left (\omega _{d}t\right ) +\frac {\zeta \omega \beta }{\omega _{d}}\sin \left (\omega _{d}t\right ) \right ] \right \} h\relax (t) \]
\(\beta =\left (\omega ^{2}-\varpi ^{2}\right ) \), \(\omega _{d}=\omega \sqrt {1-\zeta ^{2}}.\)
The solutions in modal coordinates are now found. Recall that \(\omega _{1}=\) \(15.0519,\omega _{2}=17.5624,\omega _{3}=45.5522\) and \(\zeta _{1}=0.0018\),\(\zeta _{2}=\) \(0.0219\) and \(\zeta _{3}=0.0376\)
Next step is to transform the solution to the physical coordinates using \(q_{j}={\displaystyle \sum \limits _{m=1}^{3}} \Phi \left (j,m\right ) \eta \relax (m) \), or\[ \mathbf {q}=\left [ \Phi \right ] \mathbf {\eta }\]
In component form\begin {align*} q_{1}\relax (t) & =\Phi \left (1,1\right ) \eta _{1}\relax (t) +\Phi \left (1,2\right ) \eta _{2}\relax (t) +\Phi \left (1,3\right ) \eta _{3}\relax (t) \\ q_{2}\relax (t) & =\Phi \left (2,1\right ) \eta _{1}\relax (t) +\Phi \left (2,2\right ) \eta _{2}\relax (t) +\Phi \left (2,3\right ) \eta _{3}\relax (t) \\ q_{3}\relax (t) & =\Phi \left (3,1\right ) \eta _{1}\relax (t) +\Phi \left (3,2\right ) \eta _{2}\relax (t) +\Phi \left (3,3\right ) \eta _{3}\relax (t) \end {align*}
Program was written to complete the computation and make plots. Here is the result showing plots of each of the above \(q_{i}\relax (t) \) vs. time
function nma_HW10_problem_3_EMA_545()
%solve for q(t) using modal analysis, by Nasser M. Abbasi
close all;
syms t;
N = 3;
k = [300 0 -200;0 500 300;-200 300 700]*10^3;
M = [600 400 200;400 1200 0;200 0 800];
C = [500 300 -400;300 900 600;-400 600 1300];
wF = 16;
F = [200*cos(wF*t); 0; 0];
[PHI,lam] = eig(k,M);
lam = sqrt(diag(lam));
CC = PHI'*C*PHI
F = PHI.'*F;
eta = sym(zeros(N, 1));
time_constant = zeros(3,1);
for i=1:N
w = lam(i);
b = w^2-wF^2;
zeta = CC(i,i)/(2*w);
wd = w*sqrt(1-zeta^2);
eta(i) = F(i)/(b^2+4*zeta^2*w^2*wF^2) * ...
( b*cos(wF*t)+2*zeta*w*wF*sin(wF*t)- ...
exp(-zeta*w*t)* ( b*cos(wd*t)+ zeta*w*b/wd * sin(wd*t) ) ...
);
time_constant(i) = 1/(zeta*w);
end
q=PHI*eta;
time_constant
time_constant = sum(time_constant);
% plot the generalized solutions
lims= [-0.004 0.003;
-0.002 0.007;
-0.006 0.002
];
for i=1:N
subplot(3,1,i);
ezplot(q(i),[0,100]);
ylim(lims(i,:));
title(sprintf('q(%d) solution, time constant = %f',i,time_constant));
xlabel('time (sec)');
ylabel('q(t) Newton');
end
end
From above, the time to reach steady state is about \(90\) seconds based on looking at \(q_{1}\relax (t) \) since that takes the longest time to each steady state out of the three coordinates.
The time constant for each \(\eta _{i}\relax (t) \) solution was calculated giving \(\tau _{1}=\frac {1}{\zeta _{1}\omega _{1}}=37.4471\) and \(\tau _{2}=2.602\) and \(\tau _{3}=0.58\). The first time constant \(\tau _{1}=37.4\) seconds dominated the result in the response in the physical coordinates.
This means the dominant time constant found in modal analysis is one to use to estimate how long it will take for the response in physical coordinates to reach steady state. Each modal solution contributes to each physical solution. The one with the longest time constant affects more than any other mode how long the physical solution takes to reach steady state.
\(\left [ M\right ] =\begin {bmatrix} 5 & -3\\ -3 & 4 \end {bmatrix} \)kg,\(\omega _{1}=15.68\) rad/sec,\(\omega _{2}=40.78\) rad/sec. \(\mathbf {\varphi }_{1}=\begin {Bmatrix} 1\\ 1.366 \end {Bmatrix} ,\mathbf {\varphi }_{2}=\begin {Bmatrix} 1\\ -0.366 \end {Bmatrix} \)\begin {align*} \mu _{1} & =\begin {Bmatrix} 1\\ 1.366 \end {Bmatrix} ^{T}\begin {bmatrix} 5 & -3\\ -3 & 4 \end {bmatrix}\begin {Bmatrix} 1\\ 1.366 \end {Bmatrix} =4.267\,8\\ \mu _{2} & =\begin {Bmatrix} 1\\ -0.366 \end {Bmatrix} ^{T}\begin {bmatrix} 5 & -3\\ -3 & 4 \end {bmatrix}\begin {Bmatrix} 1\\ -0.366 \end {Bmatrix} =7.731\,8 \end {align*}
Normalized eigenvectors are \begin {align*} \mathbf {\Phi }_{1} & =\frac {1}{\sqrt {\mu _{1}}}\mathbf {\varphi }_{1}=\frac {1}{\sqrt {4.267\,8}}\begin {Bmatrix} 1\\ -0.366 \end {Bmatrix} =\begin {Bmatrix} 0.484\,06\\ -0.177\,17 \end {Bmatrix} \\ \mathbf {\Phi }_{2} & =\frac {1}{\sqrt {\mu _{2}}}\mathbf {\varphi }_{2}=\frac {1}{\sqrt {7.731\,8}}\begin {Bmatrix} 1\\ -0.366 \end {Bmatrix} =\begin {Bmatrix} 0.359\,63\\ -0.131\,63 \end {Bmatrix} \end {align*}
Hence\[ \left [ \Phi \right ] =\left [ \mathbf {\Phi }_{1}\mathbf {\Phi }_{2}\right ] =\begin {bmatrix} 0.484\,06 & 0.359\,63\\ -0.177\,17 & -0.131\,63 \end {bmatrix} \]
EOM in modal coordinates is\begin {align*} \begin {bmatrix} 1 & 0\\ 0 & 1 \end {bmatrix}\begin {Bmatrix} \eta _{1}^{\prime \prime }\\ \eta _{2}^{\prime \prime }\end {Bmatrix} +\begin {bmatrix} 2\left (0.08\right ) \left (15.68\right ) & 0\\ 0 & 2\left (0.08\right ) \left (40.78\right ) \end {bmatrix}\begin {Bmatrix} \eta _{1}^{\prime }\\ \eta _{2}^{\prime }\end {Bmatrix} +\begin {bmatrix} 15.68^{2} & 0\\ 0 & 40.78^{2}\end {bmatrix}\begin {Bmatrix} \eta _{1}\\ \eta _{2}\end {Bmatrix} & =\left [ \Phi \right ] ^{T}\begin {Bmatrix} 50\sin \left (20t\right ) \\ 100\cos \left (20t\right ) \end {Bmatrix} \\\begin {bmatrix} 1 & 0\\ 0 & 1 \end {bmatrix}\begin {Bmatrix} \eta _{1}^{\prime \prime }\\ \eta _{2}^{\prime \prime }\end {Bmatrix} +\begin {bmatrix} 2.509 & 0\\ 0 & 6.524\,8 \end {bmatrix}\begin {Bmatrix} \eta _{1}^{\prime }\\ \eta _{2}^{\prime }\end {Bmatrix} +\begin {bmatrix} 245.86 & 0\\ 0 & 1663 \end {bmatrix}\begin {Bmatrix} \eta _{1}\\ \eta _{2}\end {Bmatrix} & =\begin {bmatrix} 24.203\sin \left (20.0t\right ) -17.717\cos \left (20.0t\right ) \\ 17.982\sin \left (20.0t\right ) -13.163\cos \left (20.0t\right ) \end {bmatrix} \end {align*}
The two EOMs to solve are\begin {align*} \eta _{1}^{\prime \prime }\relax (t) +2.509\eta _{1}^{\prime }\left ( t\right ) +245.86\eta _{1}\relax (t) & =24.203\sin \left (20t\right ) -17.717\cos \left (20t\right ) =\operatorname {Re}\left \{ \frac {24.203}{i}e^{i20t}\right \} +\operatorname {Re}\left \{ -17.717e^{i20t}\right \} \\ \eta _{2}^{\prime \prime }\relax (t) +6.525\eta _{2}^{\prime }\left ( t\right ) +1663\eta _{2}\relax (t) & =17.982\sin \left (20t\right ) -13.163\cos \left (20t\right ) =\operatorname {Re}\left \{ \frac {17.982}{i}e^{i20t}\right \} +\operatorname {Re}\left \{ -13.163e^{i20t}\right \} \end {align*}
Hence\begin {align*} \eta _{1}^{\prime \prime }\relax (t) +2.509\eta _{1}^{\prime }\left ( t\right ) +245.86\eta _{1}\relax (t) & =24.203\sin \left (20t\right ) -17.717\cos \left (20t\right ) =\operatorname {Re}\left \{ \left ( -24.203i-17.717\right ) e^{i20t}\right \} \\ \eta _{2}^{\prime \prime }\relax (t) +6.525\eta _{2}^{\prime }\left ( t\right ) +1663\eta _{2}\relax (t) & =17.982\sin \left (20t\right ) -13.163\cos \left (20t\right ) =\operatorname {Re}\left \{ \left ( -17.982i-13.163\right ) e^{i20t}\right \} \end {align*}
In matrix form\[ \left [ I\right ] \mathbf {\eta }^{\prime \prime }+\left [ C\right ] \mathbf {\eta }^{\prime }+\left [ K\right ] \mathbf {\eta =}\operatorname {Re}\left \{ \mathbf {F}e^{i\varpi t}\right \} \]
Where \(\varpi =20\) rad/sec. \(\mathbf {F}\) is the complex amplitude of the input \[ \mathbf {F}=\begin {Bmatrix} -24.203i-17.717\\ -17.982i-13.163 \end {Bmatrix} \]
Using method of transfer functions (since steady state response is needed), response is\[ \mathbf {\eta =}\operatorname {Re}\left \{ \mathbf {X}e^{i20t}\right \} \]
Where \[ X_{j}=\frac {F_{j}}{-\varpi ^{2}+2i\zeta _{j}\omega _{j}\varpi +\omega _{j}^{2}}\]
Steady state solutions in modal coordinates is\begin {align*} \eta _{1}\relax (t) & =\operatorname {Re}\left \{ \frac {-24.203i-17.717}{-\varpi ^{2}+2.5088i\varpi +245.86}e^{i\varpi t}\right \} \\ & =\operatorname {Re}\left \{ \frac {-24.203i-17.717}{-400+50.176i+245.86}e^{i\varpi t}\right \} \\ & =\operatorname {Re}\left \{ \left (5.77\times 10^{-2}+0.176i\right ) e^{i\varpi t}\right \} \\ \eta _{2}\relax (t) & =\operatorname {Re}\left \{ \frac {-17.982i-13.163}{-\varpi ^{2}+6.525i\varpi +1663}e^{i\varpi t}\right \} \\ & =\operatorname {Re}\left \{ \frac {-17.982i-13.163}{-400+130.5i+1663}e^{i\varpi t}\right \} \\ & =\operatorname {Re}\left \{ \left (-1.178\times 10^{-2}-1.302\times 10^{-2}i\right ) e^{i\varpi t}\right \} \end {align*}
Solutions are transformed back to normal coordinates \[ \mathbf {q}=\left [ \Phi \right ] \mathbf {\eta }\]
Hence\begin {align*} q_{j}\relax (t) & =\sum _{n}^{2}\Phi \left (j,n\right ) \eta \left ( n\right ) \\ & =\sum _{n}^{2}\Phi \left (j,n\right ) \operatorname {Re}\left \{ X\left ( n\right ) e^{i\varpi t}\right \} \\ & =\operatorname {Re}\sum _{n}^{2}\Phi \left (j,n\right ) X\relax (n) e^{i\varpi t} \end {align*}
Since\(\left [ \Phi \right ] =\) \(\begin {bmatrix} 0.484\,06 & 0.359\,63\\ -0.177\,17 & -0.131\,63 \end {bmatrix} \) then\begin {align*} q_{1}\relax (t) & =\operatorname {Re}\left (\left \{ 0.484\,06\left ( 5.77\times 10^{-2}+0.176i\right ) +0.359\,63\left (-1.178\times 10^{-2}-1.302\times 10^{-2}i\right ) \right \} e^{i20t}\right ) \\ q_{2}\relax (t) & =\operatorname {Re}\left (\left \{ -0.177\,17\left (5.77\times 10^{-2}+0.176i\right ) -0.131\,63\left ( -1.178\times 10^{-2}-1.302\times 10^{-2}i\right ) \right \} e^{i20t}\right ) \end {align*}
or\begin {align*} q_{1}\relax (t) & =\operatorname {Re}\left (\left \{ 2.369\times 10^{-2}+8.051\times 10^{-2}i\right \} e^{i20t}\right ) \\ q_{2}\relax (t) & =\operatorname {Re}\left (\left \{ -8.672\times 10^{-3}-2.947\times 10^{-2}i\right \} e^{i20t}\right ) \end {align*}
Therefore \begin {align*} Y_{1} & =2.369\times 10^{-2}+8.051\times 10^{-2}i\\ Y_{2} & =-8.672\times 10^{-3}-2.947\times 10^{-2}i \end {align*}
sectionProblem 5
First step is to determine EOM. The kinetic energy \(T\) is\[ T=\frac {1}{2}I\left (\theta ^{\prime }\right ) ^{2}+\frac {1}{2}m_{2}\left ( x^{\prime }\right ) ^{2}\]
\(I=\frac {1}{3}m_{1}L^{2}.\)Assuming small angle, stiff spring approximation and zero gravity datum at the level where pendulum is hinged, spring potential energy \(V\) is\begin {align*} V & =\frac {1}{2}k\left (x-\frac {L}{2}\theta \right ) ^{2}\\ & =\frac {1}{2}k\left (x^{2}+\frac {L^{2}}{4}\theta ^{2}-xL\theta \right ) \\ & =\theta ^{2}\left (\frac {L^{2}}{8}k\right ) +x^{2}\left (\frac {1}{2}k\right ) +x\theta \left (-\frac {kL}{2}\right ) \end {align*}
Stiffness matrix due to spring is\[ K_{spring}=\begin {bmatrix} \frac {L^{2}}{4}k & -\frac {kL}{2}\\ -\frac {kL}{2} & k \end {bmatrix} \]
Potential energy due to gravity is \(V_{g}=-mg\frac {L}{2}\cos \theta \). Hence \(V_{g_{11}}=\frac {\partial ^{2}V_{g}}{\partial ^{2}\theta }=\left (mg\frac {L}{2}\cos \theta \right ) _{\theta =0}=mg\frac {L}{2}.\) All other terms are zero. The stiffness matrix due to gravity is\[ K_{spring}=\begin {bmatrix} mg\frac {L}{2} & 0\\ 0 & 0 \end {bmatrix} \]
Combined stiffness matrix is\[ K=\begin {bmatrix} \frac {L^{2}}{4}k+mg\frac {L}{2} & -\frac {kL}{2}\\ -\frac {kL}{2} & k \end {bmatrix} \]
EOM is\[\begin {bmatrix} I & 0\\ 0 & m_{2}\end {bmatrix}\begin {Bmatrix} \theta ^{\prime \prime }\\ x^{\prime \prime }\end {Bmatrix} +\begin {bmatrix} \frac {L^{2}}{4}k+mg\frac {L}{2} & -\frac {kL}{2}\\ -\frac {kL}{2} & k \end {bmatrix}\begin {Bmatrix} \theta \\ x \end {Bmatrix} =\begin {Bmatrix} Q_{\theta }\\ Q_{x}\end {Bmatrix} \]
Generalized forces are now found. \(Q_{\theta }=FL\) since \(F\) is only external forces acting on the first d.o.f. \(\theta \) and the work done by this force is \(FL\delta \theta \) for small virtual angle. For \(Q_{x}\) work is done only by damper and acts to remove energy, hence negative in sign. \(Q_{x}\) \(=-cx^{\prime }\). The above becomes\[\begin {bmatrix} I & 0\\ 0 & m_{2}\end {bmatrix}\begin {Bmatrix} \theta ^{\prime \prime }\\ x^{\prime \prime }\end {Bmatrix} +\begin {bmatrix} \frac {L^{2}}{4}k+mg\frac {L}{2} & -\frac {kL}{2}\\ -\frac {kL}{2} & k \end {bmatrix}\begin {Bmatrix} \theta \\ x \end {Bmatrix} =\begin {Bmatrix} FL\\ -cx^{\prime }\end {Bmatrix} \]
Rearranging\[\begin {bmatrix} I & 0\\ 0 & m_{2}\end {bmatrix}\begin {Bmatrix} \theta ^{\prime \prime }\\ x^{\prime \prime }\end {Bmatrix} +\begin {bmatrix} 0 & 0\\ 0 & c \end {bmatrix}\begin {Bmatrix} \theta ^{\prime }\\ x^{\prime }\end {Bmatrix} +\begin {bmatrix} \frac {L^{2}}{4}k+m_{2}g\frac {L}{2} & -\frac {kL}{2}\\ -\frac {kL}{2} & k \end {bmatrix}\begin {Bmatrix} \theta \\ x \end {Bmatrix} =\begin {Bmatrix} FL\\ 0 \end {Bmatrix} \]
Each EOM is\begin {align*} I\theta ^{\prime \prime }+\left (\frac {L^{2}}{4}k+m_{1}g\frac {L}{2}\right ) \theta -\frac {kL}{2}x & =FL\\ m_{2}x^{\prime \prime }+cx^{\prime }-\frac {kL}{2}\theta +kx & =0 \end {align*}
Units checking: First EOM. each term must have units of torque. \(\frac {L^{2}}{4}k\theta \) have units of torque OK. \(mg\frac {L}{2}\theta \) have units of torque OK. \(kLx\) have units of torque OK.
second EOM Each term must have units of force. \(cx^{\prime }\) have units of force OK. \(kL\theta \) have units of force, OK. \(kx\) have units of force, OK.
Transfer function is now found Let \(x=\operatorname {Re}\left \{ Xe^{i\varpi t}\right \} ,\theta =\operatorname {Re}\left \{ Ye^{i\varpi t}\right \} ,F=\operatorname {Re}\left \{ \hat {F}e^{i\varpi t}\right \} .\)Substitute in the above EOM\begin {align*} \operatorname {Re}\left \{ \left [ \left (-I\varpi ^{2}Y\right ) +\left ( \frac {L^{2}}{4}k+m_{2}g\frac {L}{2}\right ) Y-\frac {kL}{2}X\right ] e^{i\varpi t}\right \} & =\operatorname {Re}\left \{ \hat {F}Le^{i\varpi t}\right \} \\ \operatorname {Re}\left \{ \left [ -m_{2}\varpi ^{2}X+ic\varpi X-\frac {kL}{2}Y+kX\right ] e^{i\varpi t}\right \} & =0 \end {align*}
Simplify\begin {align} \left (-I\varpi ^{2}+\frac {L^{2}}{4}k+m_{2}g\frac {L}{2}\right ) Y-\frac {kL}{2}X & =\hat {F}L\label {eq:3.1}\\ \left (-m_{2}\varpi ^{2}+ic\varpi +k\right ) X & =\frac {kL}{2}Y\label {eq:3.2} \end {align}
The above two equations are solved to obtain the required transfer functions \(X/F\) and \(Y/F\,\). To obtain \(Y/F,~\)the second equation solved for \(X\) in terms of \(Y\) \[ X=\frac {\frac {kL}{2}}{-m_{2}\varpi ^{2}+ic\varpi +k}Y \]
\(X\) in first equation is replaced by the giving\begin {align*} \left (-I\varpi ^{2}+\frac {L^{2}}{4}k+m_{2}g\frac {L}{2}\right ) Y-\frac {kL}{2}\frac {\frac {kL}{2}}{-m_{2}\varpi ^{2}+ic\varpi +k}Y & =\hat {F}L\\ \left (-I\varpi ^{2}+\frac {L^{2}}{4}k+m_{2}g\frac {L}{2}-\frac {\frac {k^{2}L^{2}}{4}}{-m_{2}\varpi ^{2}+ic\varpi +k}\right ) Y & =\hat {F}L \end {align*}
Hence\[ Y=\frac {1}{\left (-\frac {1}{3}m_1L\varpi ^2+\frac {L}{4}k+\frac {m_2g}{2}-\frac {k^2L/4}{-m_2\varpi ^2+ic\varpi +k}\right ) }\hat {F} \]
To obtain the transfer function \(X/F\), the second equation is solved for \(Y\) in terms of \(X\)\[ Y=\frac {\left (-m_{2}\varpi ^{2}+ic\varpi +k\right ) }{kL/2}X \]
This is substituted in the first equation giving\begin {align*} \left (-I\varpi ^{2}+\frac {L^{2}}{4}k+m_{2}g\frac {L}{2}\right ) \frac {\left ( -m_{2}\varpi ^{2}+ic\varpi +k\right ) }{kL/2}X-\frac {kL}{2}X & =\hat {F}L\\ \left [ \frac {\left (-\frac {1}{3}m_{1}L\varpi ^{2}+\frac {L}{4}k+\frac {m_{2}g}{2}\right ) \left (-m_{2}\varpi ^{2}+ic\varpi +k\right ) }{k/2}-\frac {kL}{2}\right ] X & =\hat {F}L \end {align*}
Hence\[ X=\frac {kL}{\left (-\frac {1}{3}m_1L\varpi ^2+\frac {L}{4}k+\frac {m_2g}{2}\right ) \left (-m_2\varpi ^2+ic\varpi +k\right ) -k^2L}\hat {F} \]
This complete part(a). These are the analytical expressions for the transfer functions.
Let \(m_{1}=m_{2}=1\) kg, \(k=3\) N/m, \(L=1\) m, \(g=9.81\) m/s\(^{2},\) \(c=0.1\) N-s/m\(.\)
A program was written to plot the magnitude and phase spectrums of \(x\left ( t\right ) \) and \(\theta \relax (t) \) using the above numerical values. This was done for a range of forcing frequencies to cover both natural frequencies and beyond. Natural frequencies are found by solving the eigenvalue problem \(\det \left (\left [ K\right ] -\omega ^{2}\left [ M\right ] \right ) =0\) \begin {align*} \omega _{1} & =1.1308\text { rad/sec}\\ \omega _{2} & =4.3228\text { rad/sec} \end {align*}
The magnitude and phase of each transfer function are evaluated when \(\varpi =\omega _{1}\) and when \(\varpi =\omega _{2}.F=1\) was assumed since its numerical value was not given. Result is shown below. From these plots, magnitude and phase values are determined at the natural frequencies.\begin {align*} x\relax (t) & =\operatorname {Re}\left \{ Xe^{i\varpi t}\right \} \\ \theta \relax (t) & =\operatorname {Re}\left \{ Ye^{i\varpi t}\right \} \end {align*}
Table of results
| response | magnitude at \(\omega _{1}\) | phase at \(\omega _{1}\) | magnitude at \(\omega _{2}\) | phase at \(\omega _{2}\) |
| \(x\relax (t) \) | \(4.25\) | \(-83^{0}\) | \(2.62\) | \(131.7^{0}\) |
| \(\theta \relax (t) \) | \(2.55\) | \(-80^{0}\) | \(11.5\) | \(-50^{0}\) |
| ratio | \(4.25/2.55=\) \(1.666\,7\) | \(11.5/2.62=\) \(4.\,389\,3\) | ||
Plots used to obtain these results
The function used to generate the plots
function nma_HW10_problem_5_EMA_545_spectrum()
%plots the spectrums of problem 5, HW10, by Nasser M. Abbasi
close all;
c = 0.1;
g = 9.81;
L = 1;
k = 3;
m1 = 1;
m2 = 1;
F = 1;
M = [1/3*m1*L^2 0;0 m2];
K = [L^2/4*k+m2*g*L/2 -k*L/2;-k*L/2 k];
C = [0 0;0 c];
[PHI,w] = eig(K,M);
lam = sqrt(diag(w))
I = sqrt(-1);
X = @(wf) ((k*L)./((-1/3*m1*L*wf.^2+L/4*k+m2*g/2).*(-m2*wf.^2+I*c*wf+k)- (k^2*L)))*F;
Y = @(wf) (1./((-1/3*m1*L*wf.^2+L/4*k+m2*g/2-( k^2*L./(-m2*wf.^2+I*c*wf+k)))))*F;
N = 2;
for i=1:N
figure(i);
wf = 0:0.1:6.5;
if i==1
name_='X';
tf_ = X(wf);
else
name_='Y';
tf_ = Y(wf);
end
subplot(2,1,1);
plot(wf,abs(tf_));
hold on;
line([lam(1) lam(1)],[0 5],'LineStyle','-.');
line([lam(2) lam(2)],[0 5],'LineStyle','-.');
title(sprintf('magnitude spectrum of %c, $\\omega_1=%f$, $\\omega_2=%f$',name_,lam(1),lam(2)),'interpreter','latex','FontSize',12);
xlabel('forcing frequency (rad/sec)');
ylabel(sprintf('$|%c|$',name_),'interpreter','latex','FontSize',12);
grid;
subplot(2,1,2);
plot(wf,angle(tf_));
line([lam(1) lam(1)],[-5 5],'LineStyle','-.');
line([lam(2) lam(2)],[-5 5],'LineStyle','-.');
title(sprintf('phase spectrum of X, $\\omega_1=%f$, $\\omega_2=%f$',lam(1),lam(2)),'interpreter','latex','FontSize',12);
xlabel('forcing frequency (rad/sec)');
ylabel(sprintf('$arg(%c)$',name_),'interpreter','latex','FontSize',12);
grid;
end
end
Eigenvectors \(\Phi _{1}\) and \(\Phi _{1}\) are now found, using modal analysis, which de-couples the EOM. The ratio of one component of the same eigenvector to its other component is found and compared with the result found above. The eigenvectors found are\begin {align*} \Phi _{1} & =\begin {Bmatrix} -0.5446\\ -0.9493 \end {Bmatrix} \\ \Phi _{2} & =\begin {Bmatrix} -1.6442\\ 0.3145 \end {Bmatrix} \end {align*}
The ratios are \(0.9493/0.5446=1.743\,1\) and \(1.6442/0.3145=5.228\,0\). Compare these to the ratios found
| response | magnitude at \(\omega _{1}\) | phase at \(\omega _{1}\) | magnitude at \(\omega _{2}\) | phase at \(\omega _{2}\) |
| \(x\relax (t) \) | \(4.25\) | \(-83^{0}\) | \(2.62\) | \(131.7^{0}\) |
| \(\theta \relax (t) \) | \(2.55\) | \(-80^{0}\) | \(11.5\) | \(-50^{0}\) |
| ratio | \(4.25/2.55=1.666\,7\) | \(11.5/2.62=\) \(4.389\,3\) | ||
These ratios are close to each others. Ratio \(\Phi _{1j}/\Phi _{2j}\) shows how much one dof (1) will change relative to dof (2) in mode \(j\)
Transfer functions are plotted in part(a). From magnitude spectrum of \(Y\) it is seen that \(\left \vert Y\right \vert =0\) when \(\varpi \) between \(1.5\) and \(2.0\) rad/sec and also when \(\varpi >6\) rad/sec. \(\omega _{cart}=\sqrt {\frac {k}{m_{2}}}=\sqrt {\frac {3}{1}}=1.732\,1\) rad/sec. This agrees with range found in plots. When \(\sqrt {\frac {k}{m_{2}}}=1.73\), top mass acts as vibration absorber, and rod will not oscillate when \(F\relax (t) \) is at this specific frequency.
From HW6, problem 3\[ f\relax (t) =\begin {array} [c]{lll}\frac {P}{\tau }t & & 0<t<\tau \\ 0 & & \tau <t<2\tau \end {array} \]
Let \(y_{ss}\relax (t) \) be the solution from problem 3 found using FFT technique. Let the full solution for deflection of the above pillar be \[ \chi \left (y,t\right ) =y\relax (t) \psi \relax (y) \]
\(y\relax (t) \) is the time dependent (dynamic) part of the solution. This solution is \(y_{ss}\relax (t) \) found in problem 3. \(\psi \left ( y\right ) \) is solution due to static loading. Also called the shape function. For cantilever beam with static force \(P\) at its end, deflection curve due to static loading \(P\) at end is\[ \psi \relax (x) =\frac {P}{6EI}\left (3Lx^{2}-x^{3}\right ) \]
Internal bending moment \(M\left (x,t\right ) =EI\frac {d^{2}\chi \left ( x,t\right ) }{dx^{2}}\) and direct stress \(\sigma =\) \(\frac {M\left (x,t\right ) c}{I}\) where \(c\) is the section modulus. Assume \(c=\frac {h}{2}\). For yield, let \(\sigma =40MPa\), then\begin {align*} M\left (x,t\right ) & =\frac {\sigma I}{c}\\ EI\frac {d^{2}\chi \left (x,t\right ) }{dx^{2}} & =\frac {\sigma I}{c} \end {align*}
\(I=\frac {1}{12}bh^{3}\).\begin {align*} \frac {d^{2}\chi \left (x,t\right ) }{dx^{2}} & =y_{ss}\relax (t) \frac {d^{2}}{dx^{2}}\frac {P}{6EI}\left (3Lx^{2}-x^{3}\right ) \\ & =y_{ss}\relax (t) \frac {PL}{EI} \end {align*}
Solve for \(P\) at yield\begin {align*} y_{ss}\relax (t) \frac {PL}{EI} & =\frac {\sigma _{yield}I}{c}\\ P & =\frac {\sigma _{yield}I}{y_{ss}\relax (t) \frac {h}{2}L}EI \end {align*}
\(y_{ss}\relax (t) \) from problem 3 has maximum value of \(1.8\) at \(t=10\) sec. Given numerical values in the problem and using this maximum value of \(y_{ss}\relax (t) \) then \(P\) can be found from above.
I am not sure this is the correct approach to solve this.We did not have any practice or examples on solving this type of vibration problem before. Need more time to study this subject.