Vibration isolation was based on reducing absolute acceleration of passenger under turbulent external forces. This was done by isolating the passenger from the base motion subjected to external absolute acceleration. Hence the model is based on the following diagram
Hence EQM of motion is\begin {align} my^{\prime \prime }+c\left (y^{\prime }-z^{\prime }\right ) +k\left (y-z\right ) & =0\nonumber \\ my^{\prime \prime }+cy^{\prime }+ky & =cz^{\prime }+kz \label {eq:1} \end {align}
We are given the time history of the turbulent acceleration. Hence in frequency domain we can write\[ z^{\prime \prime }=\operatorname {Re}\left \{ Z_{n}^{acc}e^{i\left (\omega _{1}n\right ) t}\right \} \]
Where \(Z_{n}^{acc}\) is the complex amplitude of the \(n^{th}\) harmonic component in the acceleration data. Let \(\omega _{1}n\equiv \varpi _{n}\) then using the above, In frequency domain Eq ?? becomes\begin {align*} \operatorname {Re}\left \{ \left (-m\varpi _{n}^{2}+i\varpi _{n}c+k\right ) Y_{n}e^{i\varpi _{n}t}\right \} & =\operatorname {Re}\left \{ \left ( c\frac {Z_{n}^{acc}}{i\varpi _{n}}+k\frac {Z_{n}^{acc}}{-\varpi _{n}^{2}}\right ) e^{i\varpi _{n}t}\right \} \\ Y_{n} & =\left (\frac {\frac {c}{i\varpi _{n}}-\frac {k}{\varpi _{n}^{2}}}{-m\varpi _{n}^{2}+i\varpi _{n}c+k}\right ) Z_{n}^{acc} \end {align*}
The above gives the transfer function between the displacement of the passenger and the external acceleration. In otherwords\[ y\relax (t) =\operatorname {Re}\left \{ \left (\frac {\frac {c}{i\varpi _{n}}-\frac {k}{\varpi _{n}^{2}}}{-m\varpi _{n}^{2}+i\varpi _{n}c+k}\right ) Z_{n}^{acc}e^{i\left (\omega _{1}n\right ) t}\right \} \]
Let \[ Y_{n}=\left (\frac {\frac {c}{i\varpi _{n}}-\frac {k}{\varpi _{n}^{2}}}{-m\varpi _{n}^{2}+i\varpi _{n}c+k}\right ) Z_{n}^{acc}\] then the transfer function is\begin {align*} H\left (\varpi _{n}\right ) & =\frac {Y_{n}}{Z_{n}^{acc}}=\frac {\frac {-ic}{\varpi _{n}}-\frac {k}{\varpi _{n}^{2}}}{-m\varpi _{n}^{2}+i\varpi _{n}c+k}\\ & =-\frac {1}{\varpi _{n}^{2}}\frac {\left (k+ic\varpi _{n}\right ) }{\left ( k-m\varpi _{n}^{2}\right ) +i\varpi _{n}c} \end {align*}
Hence phase is\[ \arg \left (H\left (\varpi _{n}\right ) \right ) =\tan ^{-1}\left ( \frac {c\varpi _{n}}{k}\right ) -\tan ^{-1}\left (\frac {\varpi _{n}c}{k-m\varpi _{n}^{2}}\right ) \]
and magnitude is\[ \left \vert H\left (\varpi _{n}\right ) \right \vert =\left \vert \frac {Y_{n}}{Z_{n}^{acc}}\right \vert =\frac {1}{\varpi _{n}}\frac {\sqrt {k^{2}+c^{2}\varpi _{n}^{2}}}{\sqrt {\left (k-m\varpi _{n}^{2}\right ) ^{2}+\left ( \varpi _{n}c\right ) ^{2}}}\]
These can be written in terms of \(\zeta \) and \(\omega _{nat}\) as follows. From \(H\left (\varpi _{n}\right ) =-\frac {1}{\varpi _{n}^{2}}\frac {\left ( k+ic\varpi _{n}\right ) }{\left (k-m\varpi _{n}^{2}\right ) +i\varpi _{n}c}\), dividing numerator and denominator by \(k=m\omega _{nat}^{2}\) and using \(c=2\zeta m\omega _{nat}\) then\[ H\left (\varpi _{n}\right ) =-\frac {1}{\varpi _{n}^{2}}\frac {\left ( 1+\frac {i2\zeta m\omega _{nat}\varpi _{n}}{m\omega _{nat}^{2}}\right ) }{\left ( 1-\frac {m\varpi _{n}^{2}}{m\omega _{nat}^{2}}\right ) +\frac {i\varpi _{n}2\zeta m\omega _{nat}}{m\omega _{nat}^{2}}}=-\frac {1}{\varpi _{n}^{2}}\frac {\left ( 1+\frac {i2\zeta \varpi _{n}}{\omega _{nat}}\right ) }{\left (1-\frac {\varpi _{n}^{2}}{\omega _{nat}^{2}}\right ) +\frac {i\varpi _{n}2\zeta }{\omega _{nat}}}\]
Let \(r_{n}=\frac {\varpi _{n}}{\omega _{nat}}\) then the above becomes\[ H\left (\varpi _{n}\right ) =-\frac {1}{\varpi _{n}^{2}}\frac {\left (1+i2\zeta r_{n}\right ) }{\left (1-r_{n}^{2}\right ) +i2r_{n}\zeta }\]
Hence\begin {align*} \left \vert H\left (\varpi _{n}\right ) \right \vert & =\frac {1}{\varpi _{n}}\frac {\sqrt {1+\left (2\zeta r_{n}\right ) ^{2}}}{\sqrt {\left (1-r_{n}^{2}\right ) ^{2}+\left (2r_{n}\zeta \right ) ^{2}}}\\ \arg \left (H\left (\varpi _{n}\right ) \right ) & =\tan ^{-1}\left (2\zeta r_{n}\right ) -\tan ^{-1}\frac {2r_{n}\zeta }{1-r_{n}^{2}} \end {align*}
The following is a plot showing the passenger absolute acceleration \(y^{\prime \prime }\relax (t) \) over the period of \(80\) seconds against the turbulent acceleration \(z^{\prime \prime }\relax (t) \). We now see that passenger absolute acceleration is close to the nominal acceleration. This was done using the following values for the vibration isolation
| \(M\) | \(100000\) kg |
| \(\zeta \) | \(0.72\) |
| \(k\) | \(38924\) N/m |
| \(c\) | \(57746\) Ns/m |
The plot on the right side is the absolute acceleration of the passenger during flight in the turbulent case.
The length of first class cabinet was estimated to be \(L=15\) meters from looking at Boeing web page.
Using Steel, Structural ASTM-A36 \(I\) beam as a cantilever beam for the implementation, then using \(k=\frac {3EI}{L^{3}}\) results in\begin {align*} 38924 & =\frac {3\left (200\times 10^{9}\right ) I}{15^{3}}\\ I & =2.1895\times 10^{-4}\text { m}^{4} \end {align*}
Using rectangle cross section \(I=\frac {bh^{3}}{12}\). Letting \(h=20\) cm, then \(b=\frac {\left (2.1895\times 10^{-4}\right ) 12}{0.2^{3}}=0.32843\) meter or \(32\) cm.
\[ Q=2000t\frac {\left (T-t\right ) }{T^{2}}\left [ h\relax (t) -h\left ( t-T\right ) \right ] \]
\begin {align*} m & =0.5\text { kg}\\ \omega _{n} & =2\pi f_{n}\\ f_{n} & =100\text { Hz} \end {align*}
Hence pulse duration is \(\frac {1}{f}=0.01\) sec\(.\)\[ my^{\prime \prime }+cy^{\prime }+ky=Q\relax (t) \]
In the frequency domain assuming that the force \(Q\relax (t) \) can be represented in its Fourier series as \[ Q\relax (t) =\operatorname {Re}\left (\sum _{n}Q_{n}e^{i\omega _{1}nt}\right ) \] where \(\omega _{1}\) is the fundamental frequency for \(Q\relax (t) \) which depends on the period we choose to select to sample over. In this example, I selected \(3T\) as the overall period to sample over so that it covers the pulse duration and an additional time to show the free vibration part as well and to compare to the analytical solution. Hence the EQM becomes\[ Y_{n}=\frac {Q_{n}}{-m\left (n\omega _{1}\right ) ^{2}+ic\left (n\omega _{1}\right ) +k}\]
\(k=\omega _{n}^{2}m\) hence dividing the numerator and denominator by \(k\) we obtain\begin {align*} Y_{n} & =\frac {\frac {Q_{n}}{K}}{\left (1-\frac {m\left (n\omega _{1}\right ) ^{2}}{\omega _{n}^{2}m}\right ) +\frac {ic\left (n\omega _{1}\right ) }{\omega _{n}^{2}m}}\\ & =\frac {1}{k}\frac {1}{\left (1-r_{n}^{2}\right ) +i2\zeta r_{n}}Q_{n} \end {align*}
where \(r_{n}=\frac {n\omega _{1}}{\omega _{n}}\).Hence response is \begin {align*} y\relax (t) & =\operatorname {Re}\left (\sum _{n}Y_{n}e^{i\omega _{1}nt}\right ) \\ & =\operatorname {Re}\left (\sum _{n}\frac {1}{k}\frac {1}{\left (1-r_{n}^{2}\right ) +i2\zeta r_{n}}Q_{n}e^{i\omega _{1}nt}\right ) \end {align*}
\(y\relax (t) \) is found by taking the IFFT of \(\sum _{n}\frac {1}{k}\frac {1}{\left (1-r_{n}^{2}\right ) +i2\zeta r_{n}}Q_{n}\).
\(Q_{n}\) values are found by taking the FFT of \(Q\relax (t) \). We start by sampling \(Q\relax (t) \). To obtain the solution for say \(t=0\cdots 3T\), then we have to assume that the period of the signal is actually \(3T\) and sample over this whole time from \(0\cdots 3T-delt\). Then we use FFT on the result. Then find the response by doing IFFT. Using \(N=128\) over \(t=0\cdots 0.03\) seconds, the following solution was obtained
%by Nasser M. Abbasi, HW 7, EMA 545
close all;
T = 0.01; %sec
duration = 3*T; %duration to find solution over
N = 128;
delT = duration/(N-1);
w1 = 2*pi/duration; %fundamental freq rad/sec
t = linspace(0,(duration-delT),N);
Qt = @(t) (2000*t.*(T-t))/T^2.*(t<=T)+0*(t>T)
subplot(2,1,1)
plot(t,Qt(t),'r-o');
hold on;
plot(0:delT:duration,Qt(0:delT:duration),'r');
title(sprintf('force Q(t) and its reponse. 16 samples, delT=%f',delT));
xlabel('time sec');
grid;
m = 0.5; %mass kg
wn = 2*pi*100; %natural freq
k = wn^2*m; %stiffness N/meter
[Q,ws] = fft_easy(Qt(t),delT);
zeta = 0.002;
I = sqrt(-1);
y = ifft_easy( (Q/k)./( (1-(ws/wn).^2) + 2*I*zeta*ws/wn),ws);
subplot(2,1,2);
plot(t,y,'r');
title(sprintf('reponse at zeta=%f',zeta));
xlabel('time sec');
grid;
For \(\zeta =0.002\) the above Matlab script was modified and the following solution resulted.
Now we compare the above with the analytical solution.
The pulse can be written as \begin {align*} F & =Q\relax (t) \left [ h\relax (t) -h\left (t-T\right ) \right ] \\ & =Q\relax (t) h\relax (t) -Q\relax (t) h\left ( t-T\right ) \end {align*}
Let \(t^{\prime }=t-T\), hence \(t=t^{\prime }+T\), therefore the above becomes\[ F=Q\relax (t) h\relax (t) -Q\left (t^{\prime }+T\right ) h\left (t^{\prime }\right ) \]
But \(Q\relax (t) =\frac {2000t\left (T-t\right ) }{T^{2}}\). Let \(\frac {2000}{T^{2}}=\beta \) since it is a constant. Hence \(Q\relax (t) =\beta t\left (T-t\right ) \). Now we write the above \(F\) as\begin {align} F & =\beta t\left (T-t\right ) h\relax (t) -\beta \left (t^{\prime }+T\right ) \left (T-\left (t^{\prime }+T\right ) \right ) h\left ( t^{\prime }\right ) \nonumber \\ & =\left (\beta Tt-\beta t^{2}\right ) h\relax (t) -\beta \left ( t^{\prime }+T\right ) \left (-t^{\prime }\right ) h\left (t^{\prime }\right ) \nonumber \\ & =\left (\beta Tt-\beta t^{2}\right ) h\relax (t) +\beta \left ( \left (t^{\prime }\right ) ^{2}+Tt^{\prime }\right ) h\left (t^{\prime }\right ) \nonumber \\ & =\beta Tth\relax (t) -\beta t^{2}h\relax (t) +\beta T\left ( t^{\prime }\right ) ^{2}+\beta Tt^{\prime }h\left (t^{\prime }\right ) \label {eq:2.1} \end {align}
So we see that the response to \(F\) will be the response to a unit impulse \(h\relax (t) \) with forcing basis functions that are \(1,t,t^{2}\). Now we can use the solution from back of the book appendix \(B\) to sum the responses in order to find the final response and compare to the FFT method.
From appendix B, the response to unit ramp \(th\relax (t) \)is\[ r\left (th\relax (t) \right ) =\frac {1}{m\omega _{n}^{3}}\left ( \omega _{n}t-2\zeta +e^{-\zeta \omega _{n}t}\left [ 2\zeta \cos \omega _{d}t-\left ( 1-2\zeta ^{2}\right ) \frac {\omega _{n}}{\omega _{d}}\sin \omega _{d}t\right ] \right ) h\relax (t) \]
and the response to quadratic \(t^{2}h\relax (t) \) is\[ s\left (t^{2}h\relax (t) \right ) =\frac {1}{m\omega _{n}^{4}}\left ( \left (\omega _{n}t\right ) ^{2}-4\zeta \omega _{n}t-2\left (1-4\zeta ^{2}\right ) +e^{-\zeta \omega _{n}t}\left [ 2\left (1-4\zeta ^{2}\right ) \cos \omega _{d}t+\left (6\zeta -8\zeta ^{3}\right ) \frac {\omega _{n}}{\omega _{d}}\sin \omega _{d}t\right ] \right ) h\relax (t) \]
Now that we have the basis solutions, we can apply them to EQ ??\begin {align*} F & =\beta T\left (r\relax (t) +r\left (t^{\prime }\right ) \right ) -\beta T\left (s\relax (t) -s\left (t^{\prime }\right ) \right ) \\ & =\beta \left (r\relax (t) +r\left (t-T\right ) \right ) -\beta T\left (s\relax (t) -s\left (t-T\right ) \right ) \\ & =\left (\beta T\right ) \frac {1}{m\omega _{n}^{3}}\left (\omega _{n}t-2\zeta +e^{-\zeta \omega _{n}t}\left [ 2\zeta \cos \omega _{d}t-\left ( 1-2\zeta ^{2}\right ) \frac {\omega _{n}}{\omega _{d}}\sin \omega _{d}t\right ] \right ) h\relax (t) \\ & +\left (\beta T\right ) \frac {1}{m\omega _{n}^{3}}\left (\omega _{n}t^{\prime }-2\zeta +e^{-\zeta \omega _{n}t^{\prime }}\left [ 2\zeta \cos \omega _{d}t^{\prime }-\left (1-2\zeta ^{2}\right ) \frac {\omega _{n}}{\omega _{d}}\sin \omega _{d}t^{\prime }\right ] \right ) h\left (t^{\prime }\right ) \\ & -\left (\beta \right ) \frac {1}{m\omega _{n}^{4}}\left (\left (\omega _{n}t\right ) ^{2}-4\zeta \omega _{n}t-2\left (1-4\zeta ^{2}\right ) +e^{-\zeta \omega _{n}t}\left [ 2\left (1-4\zeta ^{2}\right ) \cos \omega _{d}t+\left (6\zeta -8\zeta ^{3}\right ) \frac {\omega _{n}}{\omega _{d}}\sin \omega _{d}t\right ] \right ) h\relax (t) \\ & +\left (\beta \right ) \frac {1}{m\omega _{n}^{4}}\left (\left (\omega _{n}t^{\prime }\right ) ^{2}-4\zeta \omega _{n}t^{\prime }-2\left (1-4\zeta ^{2}\right ) +e^{-\zeta \omega _{n}t^{\prime }}\left [ 2\left (1-4\zeta ^{2}\right ) \cos \omega _{d}t^{\prime }+\left (6\zeta -8\zeta ^{3}\right ) \frac {\omega _{n}}{\omega _{d}}\sin \omega _{d}t^{\prime }\right ] \right ) h\left (t^{\prime }\right ) \end {align*}
In the above, \(\omega _{d}=\omega _{n}\sqrt {1-\zeta ^{2}}\). To plot this solution, the following small script was used and was run for both \(\zeta =0.2\) and \(\zeta =0.002\)
For \(\zeta =0.2\)
For \(\zeta =0.002\)
The analytical solution, using superposition agreed with the FFT solution for \(\zeta =0.2\). However, for some reason which I am not able to determine why yet, the FFT solution when \(\zeta =0.002\) did not agree with the analytical solution. The analytical solution was verified to be correct using another numerical ODE solver. So the FFT method for some reason is not giving accurate result for \(\zeta =0.002\). The same Matlab script was used for both cases. I tried increasing the sampling rate but that did not change the result. Please see Appendix for verification and the code used to plot the analytical solutions.
Let \(T\) be the kinetic energy and \(V\) be the potential energy. Then equation of motion for a generalized coordinate \(q_{i}\) is given by
\[ \frac {d}{dt}\left (\frac {\partial L}{\partial \dot {q}_{i}}\right ) -\frac {\partial L}{\partial q_{i}}=Q_{i}\]
Where \(L\) is the Lagrangian \(L=T-V\) and \(Q_{i}\) is the generalized force in the \(q_{i}\) direction.
Assuming \(x_{2}>x_{1}\) and masses are moving to the right. For \(x_{1}\) we obtain\begin {align*} T & =\frac {1}{2}m_{1}\dot {x}_{1}^{2}+\frac {1}{2}m_{1}\dot {x}_{2}^{2}\\ V & =\frac {1}{2}k_{1}x_{1}^{2}+\frac {1}{2}k_{2}\left (x_{2}-x_{1}\right ) ^{2}+\frac {1}{2}k_{4}x_{1}^{2}+\frac {1}{2}k_{3}x_{2}^{2}\\ Q_{1} & =F\\ Q_{2} & =0 \end {align*}
Hence \begin {align*} L & =T-V\\ & =\frac {1}{2}m_{1}\dot {x}_{1}^{2}+\frac {1}{2}m_{1}\dot {x}_{2}^{2}-\left ( \frac {1}{2}k_{1}x_{1}^{2}+\frac {1}{2}k_{2}\left (x_{2}-x_{1}\right ) ^{2}+\frac {1}{2}k_{4}x_{1}^{2}+\frac {1}{2}k_{3}x_{2}^{2}\right ) \end {align*}
\begin {align*} \frac {\partial L}{\partial \dot {x}_{1}} & =m_{1}\dot {x}_{1}\\ \frac {d}{dt}\left (\frac {\partial L}{\partial \dot {x}_{1}}\right ) & =m_{1}\ddot {x}_{1}\\ \frac {\partial L}{\partial x_{1}} & =-k_{1}x_{1}-k_{2}\left (x_{2}-x_{1}\right ) \left (-1\right ) -k_{4}x_{1} \end {align*}
and
\begin {align*} \frac {\partial L}{\partial \dot {x}_{2}} & =m_{1}\dot {x}_{2}\\ \frac {d}{dt}\left (\frac {\partial L}{\partial \dot {x}_{2}}\right ) & =m_{1}\ddot {x}_{2}\\ \frac {\partial L}{\partial x_{2}} & =-k_{2}\left (x_{2}-x_{1}\right ) \relax (1) -k_{3}x_{2} \end {align*}
Hence the 2 EOM are for \(x_{1}\)\begin {align*} \frac {d}{dt}\left (\frac {\partial L}{\partial \dot {x}_{1}}\right ) -\frac {\partial L}{\partial x_{1}} & =F\\ m_{1}\ddot {x}_{1}-\left (-k_{1}x_{1}+k_{2}\left (x_{2}-x_{1}\right ) -k_{4}x_{1}\right ) & =F\\ m_{1}\ddot {x}_{1}+k_{1}x_{1}-k_{2}\left (x_{2}-x_{1}\right ) +k_{4}x_{1} & =F \end {align*}
Therefore EOM 1\[ m_{1}\ddot {x}_{1}+\left (k_{1}+k_{2}+k_{4}\right ) x_{1}-k_{2}x_{2}=F \]
and for \(x_{2}\)\begin {align*} \frac {d}{dt}\left (\frac {\partial L}{\partial \dot {x}_{2}}\right ) -\frac {\partial L}{\partial x_{2}} & =0\\ m_{1}\ddot {x}_{2}-\left (-k_{2}\left (x_{2}-x_{1}\right ) -k_{3}x_{2}\right ) & =0\\ m_{1}\ddot {x}_{2}+k_{2}\left (x_{2}-x_{1}\right ) +k_{3}x_{2} & =0 \end {align*}
Hence EOM 2\[ m_{1}\ddot {x}_{2}+\left (k_{2}+k_{3}\right ) x_{2}-k_{2}x_{1}=0 \]
Hence in Matrix form EOM are\begin {align*} MX^{\prime \prime }+KX & =Q\\\begin {pmatrix} m_{1} & 0\\ 0 & m_{2}\end {pmatrix}\begin {pmatrix} x_{1}^{\prime \prime }\\ x_{2}^{\prime \prime }\end {pmatrix} +\begin {pmatrix} \left (k_{1}+k_{2}+k_{4}\right ) & -k_{2}\\ -k_{2} & \left (k_{2}+k_{3}\right ) \end {pmatrix}\begin {pmatrix} x_{1}\\ x_{2}\end {pmatrix} & =\begin {pmatrix} F\\ 0 \end {pmatrix} \end {align*}
If \(m_{2}\) do not exist, then this means the springs \(k_{2}\) and \(k_{3}\) do not have a mass between them and so these need to be replaced by single spring, say \(k_{5}\) found by finding equivalent spring in series
\begin {align*} \frac {1}{k_{5}} & =\frac {1}{k_{2}}+\frac {1}{k_{3}}\\ k_{5} & =\frac {k_{3}+k_{2}}{k_{2}k_{3}} \end {align*}
From above, EQM for \(m_{1}\) becomes\[ m_{1}\ddot {x}_{1}+\left (k_{1}+\overset {k_{5}}{\overbrace {\left (\frac {k_{3}+k_{2}}{k_{2}k_{3}}\right ) }}+k_{4}\right ) x_{1}=F \]
So now \(k_{4}\) and \(k_{4}\) are in parallel, hence we replace \(k_{5}+k_{4}\) by \(k_{6}\) found from\begin {align*} k_{6} & =k_{5}+k_{4}\\ & =\left (\frac {k_{3}+k_{2}}{k_{2}k_{3}}\right ) +k_{4}\\ k_{6} & =\frac {k_{3}+k_{2}+k_{2}k_{3}k_{4}}{k_{2}k_{3}} \end {align*}
Hence EQM for \(m_{1}\) now becomes\[ m_{1}\ddot {x}_{1}+\left (k_{1}+\overset {k_{6}}{\overbrace {\frac {k_{3}+k_{2}+k_{2}k_{3}k_{4}}{k_{2}k_{3}}}}\right ) x_{1}=F \]
and finally\[ m_{1}\ddot {x}_{1}+\frac {k_{3}+k_{2}+k_{2}k_{3}k_{4}+k_{1}k_{2}k_{3}}{k_{2}k_{3}}x_{1}=F \]