The APDL was modified to use solid pipe288 and put the mass element at the location as shown in the problem statement. The following are the four modes generated by ANSYS
| set number | mode | frequency (Hz) |
| \(1\) | Torsion | \(10.437\) |
| \(2\) | First transverse (bending) | \(14.1815\) |
| \(3\) | Second transverse (bending) | \(35.384\) |
| \(4\) | First longitudinal (axial) | \(447.98\) |
The following are the four plots showing the mode shapes for each of the above modes
The above result was next compared to the analytical result that was done in HW 3, by using the numerical value given in this problem. The numerical values for this problem are listed here
| variable name | numerical value |
| \(L\) (length of pipe) | \(6\) ft |
| \(a\) | \(2\) ft |
| \(b\) | \(4\) ft |
| \(d\) (diameter of pipe) | \(1.2\) in \(=\frac{1.2}{12}=\) \(0.1\) ft |
| \(W\) (weight of flywheel) | \(100\) lb |
| \(r\) (outer radius of flywheel) | \(16\) in \(=\frac{16}{12}=\) \(1.3333\) ft |
| \(r_{f}\) (radius of gyration) | \(\sqrt{\frac{r^{2}}{2}}=\sqrt{\frac{1.3333^{2}}{2}}=\) \(0.94279\) ft |
| \(E\) (Elastic modulus of pipe material, steel) | \(29007547.546\times 144\) psf (\(200\) GPa) |
| \(G\) (shear modulus for pipe material, steel) | \(11196913.353\times 144\) psf (\(7.2\) GPa) |
| Poisson’s ratio for steel | \(0.295\) |
| \(I\) area moment of inertia for pipe section | \(\frac{\pi }{4}\left ( \frac{d}{2}\right ) ^{4}=4.90874\times 10^{-6}\) ft\(^{4}\) |
| \(I_{flywheel}\) mass moment of inertial of flywheel | \(\frac{W}{g}r_{f}^{2}=5.52105\) slug-ft\(^{2}\) |
The above values were now used in the derivations from HW3 to obtain numerical values for the natural frequencies. The following are the results obtained (using analytical result from HW3 derivation)
| mode | Analytical result | Numerical calculation \(\omega _{n}\) in rad/sec | Hz |
| Torsion | \(\omega _{n}=\sqrt{\frac{gG\pi d^{4}}{32Wr_{f}^{2}}\left ( \frac{1}{a}+\frac{1}{b}\right ) }\) | \(\sqrt{\frac{\left ( 32.2\right ) \left ( 11196913.353\times 144\right ) \pi \left ( 0.1\right ) ^{4}}{32\left ( 100\right ) \left ( 0.94279\right ) ^{2}}\left ( \frac{1}{2}+\frac{1}{4}\right ) }=\) \(65.58\) | \(\allowbreak 10.437\) |
| bending (1) | \(\omega _{n}=\sqrt{\frac{3gEI}{W}\left ( \frac{L}{ab}\right ) ^{3}}\) | \(\sqrt{\frac{3\left ( 32.2\right ) \left ( 29007547.546\times 144\right ) \left ( 4.90874\times 10^{-6}\right ) }{100}\left ( \frac{6}{\left ( 2\right ) \left ( 4\right ) }\right ) ^{3}}=\) \(91.412\) | \(14.549\) |
| axial | \(\omega _{n}=\sqrt{\frac{gAE}{W}\left ( \frac{1}{a}+\frac{1}{b}\right ) }\) | \(\sqrt{\frac{\left ( 32.2\right ) \left ( \pi \left ( \frac{0.1}{2}\right ) ^{2}\right ) \left ( 29007547.546\times 144\right ) }{100}\left ( \frac{1}{2}+\frac{1}{4}\right ) }=\) \(2814.8\) | \(447.99\) |
The following table compares the above analytical result with the ANSYS result shown earlier with the percentage error
| mode | ANSYS result (Hz) | Analytical result (Hz) | error percentage |
| Torsion | \(10.437\) | \(\allowbreak 10.437\) | \(0\%\) |
| First bending | \(14.1815\) | \(14.549\) | \(\left ( \frac{14.1815-14.549}{14.1815}\right ) \times 100=\) \(2.59\%\) |
| First axial | \(447.98\) | \(447.99\) | \(\frac{447.98-447.99}{447.98}\times 100=0.002\%\) |
The mode that has most error is the first bending (transverse) mode. This was the case also in HW7 ANSYS problem. ANSYS result is the more accurate one. The analytical result for this mode was derived The transverse case uses stiffness \(3EI\left ( \frac{L}{ab}\right ) ^{3}\) due to load at \(a\) distance from one end of fixed-free beam and \(b\) distance from the other end of the fixed beam. But this derivation does not account for any bending rotation in the beam as the ANSYS result would do.
The first step is to determine the natural frequency of the system. Since the springs are in parallel then \[ k_{eq}=6k \] And the equivalent mass is \(m_{eq}=\frac{W}{g}\) where \(W=700\) N. Hence\[ \omega _{n}=\sqrt{\frac{k_{eq}}{m_{eq}}}=\sqrt{\frac{6k}{\frac{W}{g}}}=\sqrt{\frac{6\left ( 6000\right ) }{\frac{700}{9.81}}}=22.461\text{ rad/sec}\] Since this is undamped system, then the steady state solution (particular solution) is given by\begin{equation} y_{p}\left ( t\right ) =\frac{x_{st}}{\sqrt{\left ( 1-r^{2}\right ) ^{2}}}\cos \omega t \tag{1} \end{equation} Where \(r=\frac{\omega }{\omega _{n}}\) and \(\omega \) is the driving frequency, which is \[ \omega =1000\left ( \frac{2\pi }{rev}\right ) \left ( \frac{\min }{60}\right ) =1000\left ( \frac{2\pi }{60}\right ) =104.72\text{ rad/sec}\] From (1), we see that the maximum steady state response is \begin{equation} y_{ss}=\frac{x_{st}}{\sqrt{\left ( 1-r^{2}\right ) ^{2}}} \tag{2} \end{equation} We now just need to determine \(x_{st}\) which is the static deflection. Let \(m_{0}\) be the unbalanced mass which is spinning inside, and let \(e\) be the radius around the spin axis. Therefore, and assuming \(\omega \) is constant, this mass will have only radial acceleration towards the center of \(e\omega ^{2}\) and therefore it will induce a centripetal force \(m_{0}e\omega ^{2}\).
From the above we see that the vertical force is\[ F\left ( t\right ) =\overset{F_{0}}{\overbrace{m_{0}e\omega ^{2}}}\sin \theta \left ( t\right ) \] Hence the static deflection is\[ x_{st}=\frac{F_{0}}{k_{eq}}=\frac{m_{0}e\omega ^{2}}{6k}\] Substituting this into (2) gives\begin{equation} y_{ss}=\frac{\frac{m_{0}e\omega ^{2}}{6k}}{\sqrt{\left ( 1-r^{2}\right ) ^{2}}}=\frac{m_{0}e\omega ^{2}}{6k\sqrt{\left ( 1-r^{2}\right ) ^{2}}} \tag{3} \end{equation} But \(r\) is\[ r=\frac{\omega }{\omega _{n}}=\frac{104.72}{22.461}=4.662\,3 \] Since \(r>1\) then we now can simplify \(\sqrt{\left ( 1-r^{2}\right ) ^{2}}=r^{2}-1\) and (3) becomes\[ y_{ss}=\frac{m_{0}e\omega ^{2}}{6k\left ( r^{2}-1\right ) }\] Since \(\ \)we want to limit deflection to \(5\) mm peak to peak, then we want to limit \(y_{ss}=2.5\) mm (which is half of the peak-to-peak). The above equation becomes\begin{align*} 2.5\times 10^{-3} & =\frac{m_{0}e\left ( 104.72\right ) ^{2}}{6\left ( 6000\right ) \left ( 4.662\,5^{2}-1\right ) }\\ & =\frac{m_{0}e\left ( 104.72\right ) ^{2}}{36000\left ( 20.739\right ) }\\ & =\frac{m_{0}e\left ( 104.72\right ) ^{2}}{7.466\times 10^{5}} \end{align*}
Solving for unbalance \(m_{0}e\) gives \[ m_{0}e=\frac{\left ( 2.5\times 10^{-3}\right ) \left ( 7.466\times 10^{5}\right ) }{\left ( 104.72\right ) ^{2}}\] Or \[ \fbox{$m_0e=0.1702$ kg-meter}\] This means to limit \(m_{0}e\) below this value in order to limit vibration to \(5\) mm, peak-to-peak.
The first step is to make a FBD and corresponding inertia diagram Where it is assumed the left spring is in tension and the right side spring is in compression.
Taking moments around the pivot \(o\) where the bar is rotating around, and using anti-clockwise as positive gives (this assumes small angle approximation) \begin{align} \sum M & =I_{o}\ddot{\theta }\nonumber \\ -k\left ( \frac{L}{4}\theta \right ) \frac{L}{4}-k\left ( \frac{3L}{4}\theta \right ) \frac{3L}{4}+kx\left ( t\right ) \left ( \frac{3L}{4}\right ) & =I_{o}\ddot{\theta }\tag{1} \end{align}
But \(I_{o}\) is the mass moment of inertia around \(o\), which is\begin{align*} I_{o} & =\overset{I_{cg}}{\overbrace{\frac{1}{12}mL^{2}}}+\overset{\text{parallel axis}}{\overbrace{m\left ( \frac{1}{4}L\right ) ^{2}}}\\ & =\frac{7}{48}L^{2}m \end{align*}
Therefore the equation of motion (1) becomes\begin{align} \frac{7}{48}L^{2}m\ddot{\theta } & =-k\left ( \frac{L^{2}}{16}\theta +\frac{9L^{2}}{16}\theta \right ) +k\frac{3L}{4}x\left ( t\right ) \nonumber \\ L^{2}m\ddot{\theta }+\theta \left ( k\frac{10}{16}L^{2}\right ) \frac{48}{7} & =k\frac{48}{7}\left ( \frac{3L}{4}\right ) x\left ( t\right ) \nonumber \\ m\ddot{\theta }+\theta \left ( \frac{30}{7}k\right ) & =k\frac{36}{7}\frac{1}{L}x\left ( t\right ) \tag{2} \end{align}
Therefore \[ \fbox{$\omega _n=\sqrt{\frac{30}{7}\frac{k}{m}}$}\] We now need to expand \(x\left ( t\right ) \) in Fourier series. \(x\left ( t\right ) \) has period of \(\tau \). This is not even and not odd function. \[ x\left ( t\right ) =\frac{X}{\tau }t \] Hence\begin{align*} a_{0} & =\frac{1}{\frac{\tau }{2}}\int _{0}^{\tau }\frac{X}{\tau }tdt=\frac{2}{\tau }\frac{X}{\tau }\left ( \frac{t^{2}}{2}\right ) _{0}^{\tau }=\frac{X}{\tau ^{2}}\tau ^{2}=X\\ a_{n} & =\frac{1}{\frac{\tau }{2}}\int _{0}^{\tau }\frac{X}{\tau }t\cos \left ( \frac{2\pi }{\tau }nt\right ) dt\\ & =\frac{2}{\tau }\frac{X}{\tau }\int _{0}^{\tau }t\cos \left ( \frac{2\pi }{\tau }nt\right ) dt\\ & =\frac{2}{\tau }\frac{X}{\tau }\left ( 0\right ) \\ & =0 \end{align*}
And\begin{align*} b_{n} & =\frac{1}{\frac{\tau }{2}}\int _{0}^{\tau }\frac{X}{\tau }t\sin \left ( \frac{2\pi }{\tau }nt\right ) dt\\ & =\frac{2}{\tau }\frac{X}{\tau }\int _{0}^{\tau }t\sin \left ( \frac{2\pi }{\tau }nt\right ) dt\\ & =\frac{2}{\tau }\frac{X}{\tau }\left ( -\frac{\tau ^{2}}{2n\pi }\right ) \\ & =-\frac{X}{n\pi } \end{align*}
Hence \begin{align*} x\left ( t\right ) & \approx \frac{a_{0}}{2}+\sum _{n=1}^{\infty }b_{n}\sin \left ( \frac{2\pi }{\tau }nt\right ) \\ & \approx \frac{X}{2}-\frac{X}{\pi }\sum _{n=1}^{\infty }\frac{1}{n}\sin \left ( \frac{2\pi }{\tau }nt\right ) \\ & \approx \frac{X}{2}-\frac{X}{\pi }\sum _{n=1}^{\infty }\frac{1}{n}\sin \left ( \frac{2\pi }{\tau }nt\right ) \end{align*}
To verify this solution, the above is plotted for number of terms to see if it will approximate the original \(x\left ( t\right ) .\)
Now we go back to the original equation of motion (2), and replace \(x\left ( t\right ) \) by its Fourier series expansion\begin{align} m\ddot{\theta }+\theta \left ( \frac{30}{7}k\right ) & =k\frac{36}{7}\frac{1}{L}\left ( \frac{X}{2}-\frac{X}{\pi }\sum _{n=1}^{\infty }\frac{1}{n}\sin \left ( \frac{2\pi }{\tau }nt\right ) \right ) \nonumber \\ & =k\frac{18}{7}\frac{X}{L}-k\frac{1}{\pi }\frac{36}{7}\frac{X}{L}\left ( \sin \left ( \frac{2\pi }{\tau }t\right ) +\frac{1}{2}\sin \left ( \frac{2\pi }{\tau }2t\right ) +\frac{1}{3}\sin \left ( \frac{2\pi }{\tau }3t\right ) +\cdots \right ) \nonumber \\ & =k\frac{18}{7}\frac{X}{L}-k\frac{1}{\pi }\frac{36}{7}\frac{X}{L}\left ( \sin \left ( \omega t\right ) +\frac{1}{2}\sin \left ( 2\omega t\right ) +\frac{1}{3}\sin \left ( 3\omega t\right ) +\frac{1}{4}\sin \left ( 4\omega t\right ) +\cdots \right ) \tag{3} \end{align}
Linearity is now used to find the solution to the above by adding the the steady state response to each of the terms. The steady state response to the first term above, which is \(\frac{18}{7}k\frac{X}{mL}\) is the steady state response to the ODE\[ m\ddot{\theta }+\theta \left ( \frac{30}{7}k\right ) =\left ( \frac{18}{7}k\frac{X}{L}\right ) \] Which Is given by\[ y_{ss}=\left ( k\frac{18}{7}\frac{X}{L}\right ) \frac{1}{k_{eq}}\] But \(k_{eq}=\frac{30}{7}k\), therefore\begin{align*} y_{ss} & =\left ( \frac{18}{7}k\frac{X}{L}\right ) \frac{7}{30k}\\ & =\frac{9}{15}\frac{X}{L} \end{align*}
This is the response to only the first term in (3). Now we do the same for each of the trig terms. But we only need to consider one general term. The ODE we will look at now is\begin{align*} m\ddot{\theta }+\theta \left ( \frac{30}{7}k\right ) & =k\frac{1}{\pi }\frac{36}{7}\frac{X}{L}\sum _{n=1}^{\infty }\frac{1}{n}\sin \left ( \frac{2\pi }{\tau }nt\right ) \\ & =k\frac{1}{\pi }\frac{36}{7}\frac{X}{L}\left ( \sin \left ( \frac{2\pi }{\tau }t\right ) +\frac{1}{2}\sin \left ( \frac{2\pi }{\tau }2t\right ) +\frac{1}{3}\sin \left ( \frac{2\pi }{\tau }3t\right ) +\cdots \right ) \\ & =k\frac{1}{\pi }\frac{36}{7}\frac{X}{L}\left ( \sin \left ( \omega t\right ) +\frac{1}{2}\sin \left ( 2\omega t\right ) +\frac{1}{3}\sin \left ( 2\omega t\right ) +\cdots \right ) \end{align*}
Considering one general term\begin{align} m\ddot{\theta }+\theta \left ( \frac{30}{7}k\right ) & =k\left ( \frac{1}{\pi }\frac{36}{7}\frac{X}{L}\frac{1}{n}\right ) \sin \left ( n\omega t\right ) \nonumber \\ & =F_{0}\sin \left ( n\omega t\right ) \tag{4} \end{align}
Where \begin{align} F_{0} & =\left ( k\frac{1}{\pi }\frac{36}{7}\frac{X}{L}\frac{1}{n}\right ) \nonumber \\ x_{st} & =\frac{F_{0}}{k_{eq}}\nonumber \\ & =\frac{k\frac{1}{\pi }\frac{36}{7}\frac{X}{L}\frac{1}{n}}{\frac{30}{7}k}\nonumber \\ & =\frac{6}{5\pi L}\frac{X}{n}\tag{5} \end{align}
We know the steady state (particular) solution for (4) is \begin{equation} \theta _{ss}\left ( t\right ) =\frac{x_{st}}{\left ( 1-\left ( nr\right ) ^{2}\right ) }\sin \left ( n\omega t\right ) \tag{6} \end{equation} Where \(r\) is\begin{equation} r=\frac{\omega }{\omega _{n}}=\frac{\frac{2\pi }{\tau }}{\sqrt{\frac{30}{7}\frac{k}{m}}}=\frac{2\pi }{\tau \sqrt{\frac{30}{7}\frac{k}{m}}}\tag{7} \end{equation} The above is the steady state response for the \(n^{th}\) term. So the total response is the sum of all these responses. Putting all this together, we now obtain the steady state solution as\begin{equation} \theta _{ss}\left ( t\right ) =k\frac{9X}{15kL}-\sum _{n=1}^{\infty }\frac{x_{st}}{\left ( 1-\left ( nr\right ) ^{2}\right ) }\sin \left ( n\omega t\right ) \tag{8} \end{equation} Where \(x_{st}\) is given (5) and \(r\) is given by (7) and \(\omega =\frac{2\pi }{\tau }\). To try verify the above, it is plotted using the following values \(X=1,L=10\) meter,\(k=100\) N/m,\(\tau =3\sec \) and \(m=5\) kg. This is the result (for \(30\) terms in Fourier sum)
