The angle \theta is the kinematic Euler angle that measures the pitch of the airplane relative to the ground. The angle of attack \alpha is an aerodynamic angle that measures the angle between the zero lift line and the velocity direction V of the airplane at any moment. The relation between \theta and \alpha is \theta =\gamma +\alpha
Equation 4.9,18 in text book represents the longitudinal equation of motion using the four longitudinal variables \left \{ u,w,q,\theta \right \} . In pure pitching mode, the two degrees of freedom are q and w (similar to short-period mode). The pure pitching mode is the short-period mode in the limit as m\longrightarrow \infty . From Mechanics of flight, Warren Phillips, page 873 :
Thus we see that as the dimensionless mass becomes large compared to the dimensionless moment of inertia, the short-period motion associated with free flight becomes pure pitching motion.
Therefore, to obtain the pure pitching mode, The short-period mode is used and then the limit is taken to obtain the pure pitching mode equations. The short-period mode is given by the following 2\times 2 system. (Equation 4.19,18 in the text book)\begin{pmatrix} \dot{w}\\ \dot{q}\end{pmatrix} =\begin{pmatrix} \frac{Z_{w}}{m-Z_{\dot{w}}} & \frac{Z_{q}+mu_{o}}{m-Z_{\dot{w}}}\\ \frac{1}{I_{y}}\left ( M_{w}+\frac{M_{\dot{w}}Z_{w}}{m-Z_{\dot{w}}}\right ) & \frac{1}{I_{y}}\left ( M_{q}+\frac{M_{\dot{w}}\left ( Z_{q}+mu_{o}\right ) }{m-Z_{\dot{w}}}\right ) \end{pmatrix}\begin{pmatrix} w\\ q \end{pmatrix}
Let the natural frequency \boxed{ \omega _{n}=\sqrt{-M_{w}\frac{u_{o}}{I_{y}}} }
Where \begin{align*} n & =\frac{1}{2I_{y}}\left ( M_{q}+M_{\dot{w}}u_{o}\right ) \\ \omega & =\frac{1}{2I_{y}}\sqrt{-\left ( M_{q}+M_{\dot{w}}u_{o}\right ) ^{2}-4I_{y}M_{w}u_{o}} \end{align*}
Hence\begin{align*} \lambda _{1} & =n+i\omega \\ \lambda _{2} & =n-i\omega \end{align*}
The period T is given by \boxed{ T=\frac{2\pi }{\omega }=\frac{4\pi I_{y}}{\sqrt{-\left ( M_{q}+M_{\dot{w}}u_{o}\right ) ^{2}-4I_{y}M_{w}u_{o}}} }
Hint: refer to Eq. (6.8,6) in the book.
Solution
Table 7.2 in the text book is
From equation 6.8,6 in the book E=g\left [ \left ( C_{l_{\beta }}C_{n_{r}}-C_{l_{r}}C_{n_{\beta }}\right ) \cos \theta _{0}+\left ( C_{l_{p}}C_{n_{\beta }}-C_{l_{\beta }}C_{n_{p}}\right ) \sin \theta _{0}\right ]
For each different value of \theta _{0}, C_{L} is changed and new value for E is obtained. This was done three times, for \theta _{0}=-10^{0},0,+10^{0}. The standard lift coefficient equation is used to obtain the corresponding value of speed for each C_{L} C_{L}=\frac{L}{\frac{1}{2}\rho V^{2}S}
Hence\begin{align*} V & =\sqrt{\frac{2\left ( 10675\right ) }{\left ( 1.225\right ) \left ( 14.9\right ) }}\sqrt{\frac{1}{C_{L}}}\\ & =34.201\sqrt{\frac{1}{C_{L}}}\text{ m/sec}\\ & =\boxed{112.18\text{ ft/sec}} \end{align*}
These are the equations to plot E vs. V. For each C_{L}, V and E are calculated using the above. A small program was written to do this. The result is in figure 2.26
The plot in figure 2.26 shows that spiral mode can become unstable depending on the speed of the airplane.
The spiral mode characteristic equation was simplified to D\lambda +E=0 as discussed in the textbook, page 193. E represents the static stability (as the case with the full, un-simplified characteristic equation). The condition for static stability is that E should remain positive. If E changes from positive to negative, which means a root has switched from negative to positive, then the response becomes unstable (diverges).
From the above plot, when the reference angle of climb \theta _{0} was larger than zero, this mode became unstable when the airplane slowed down to below critical value. The larger \theta _{0} became, the larger this critical value became.
For example, when \theta _{0}=0, this critical speed was about 27 m s−1, but when \theta _{0}=\SI{10}{\degree }, the critical speed was above 39 m s−1 . When the reference angle of climb is negative, this mode remained stable for all the speed range shown since E was positive throughout.
Solution:
The first step is to evaluate E for the given jet using the stability derivatives in table 6.6.
Where E>0 implies the following (per equation 6.8,6 on page 194 of textbook)\begin{equation} \left ( C_{l_{\beta }}C_{n_{r}}-C_{l_{r}}C_{n_{\beta }}\right ) \cos \theta _{0}+\left ( C_{l_{p}}C_{n_{\beta }}-C_{l_{\beta }}C_{n_{p}}\right ) \sin \theta _{0}>0 \tag{1} \end{equation}
| C_{l} | The rolling moment coefficient C_{l}=\frac{L}{\frac{1}{2}\rho V^{2}Sb} where L here is rolling moment and not lift |
| C_{l_{\beta }} | \frac{\partial C_{l}}{\partial \beta } where \beta is the side slip angle |
| C_{l_{r}} | \frac{\partial C_{l}}{\partial r} where r is yaw rate |
| C_{l_{p}} | \frac{\partial C_{l}}{\partial p} rolling moment coefficient of propulsion units |
| C_{n} | Yawning moment coefficient C_{n}=\frac{N}{\frac{1}{2}\rho V^{2}Sb}where N is the yawning moment |
| C_{n_{r}} | \frac{\partial C_{n}}{\partial r} where r is yaw rate |
| C_{n_{\beta }} | \frac{\partial C_{n}}{\partial \beta } where \beta is the side slip angle |
| C_{n_{p}} | \frac{\partial C_{n}}{\partial p}Yawning moment coefficient of propulsion units |
Using table 6.6 gives \begin{align*} C_{l_{\beta }} & =-0.2797\\ C_{n_{r}} & =-0.2737\\ C_{l_{r}} & =0.304\\ C_{n_{\beta }} & =0.1946\\ C_{l_{p}} & =-0.3295\\ C_{n_{p}} & =-0.04073 \end{align*}
Substituting all these values in (6.8,6) gives \left ( C_{l_{\beta }}C_{n_{r}}-C_{l_{r}}C_{n_{\beta }}\right ) \cos \theta _{0}+\left ( C_{l_{p}}C_{n_{\beta }}-C_{l_{\beta }}C_{n_{p}}\right ) \sin \theta _{0}>0
Therefore\begin{align*} \theta _{0} & >\tan ^{-1}\left ( 0.230\,36\right ) \\ & > \SI{12.97}{\degree } \end{align*}
The above implies that the climb angle has to be larger than 12.97° to insure that E remains positive and the jet remain statically stable in spiral mode at any speed.
E is now examined to see how it is depends on \Gamma (dihedral angle). The expression for E above does not have \Gamma in its as shown, but \Gamma comes into play when the coefficients are replaced by their expressions in appendix B.9 and B.11. Before doing this, the term multiplied by \sin \theta _{0} are cancelled inthe above, since the jet is in a horizontal flight or \theta _{0}=0. E>0 now implies the following \begin{equation} \left ( C_{l_{\beta }}C_{n_{r}}-C_{l_{r}}C_{n_{\beta }}\right ) >0 \tag{2} \end{equation}
Looking at figure B.9,5 shows that K_{M_{\Gamma }}>0 for all ranges defined over M\cos \Lambda _{c/2}.
This means \left ( \frac{C_{l_{\beta }}}{\Gamma }K_{M_{\Gamma }}\right ) =\left ( \operatorname{negative}\times \operatorname{positive}\right ) =\operatorname{negative}. Therefore, as \Gamma increases, C_{l_{\beta }} becomes more negative (since \Gamma is positive). To see this more clearly, (10) is written as\begin{equation} -\overbrace{\Gamma \left ( \operatorname{negative}\right )}^{C_{l_{\beta }}} -\overbrace{\Gamma \left ( \operatorname{positive}\right ) }^{C_{l_{r}}} >0 \tag{11} \end{equation}
It is assumed in all of this that \Gamma is a positive angle and remain positive.
For horizontal flight \theta =0, as the \Gamma angle increases, its effect on C_{l_{r}} is to make the airplane become unstable in spiral mode, while \Gamma effect on C_{l_{\beta }} is to make the airplane become stable.
Solution method summary: It is required to write \dot{x}=Ax for the lateral mode. This is equation 4.9,19 on page 113 of the textbook. These are expressed using dimensional derivatives Y_{v},Y_{m}\,\cdots . Table 7.2 is used, and using non-dimensional derivatives C_{y_{\beta }},C_{l_{\beta }},\cdots , with table 4.5 on page 118, the numerical values in the dimensional matrix A are found. Having obtained \dot{x}=Ax in numerical form, Matlab was used to obtain the eigenvalues and eigenvectors. The above is done for both C_{L}=0 and C_{L}\neq 0. For the case of C_{L}=0 or g=0, the A matrix becomes 3\times 3 while for g\neq 0 the A matrix remained 4\times 4.
Solution:
Table 7.2 in the text book is
Table 4.5, page 118 of the textbook is used to convert from dimensional to non-dimensional
Equation 4.9.19 for lateral mode, in dimensional form is\begin{equation} \begin{pmatrix} \dot{v}\\ \dot{p}\\ \dot{r}\\ \dot{\phi }\end{pmatrix} =\begin{pmatrix} \frac{Y_{v}}{m} & \frac{Y_{p}}{m} & \left ( \frac{Y_{r}}{m}-u_{0}\right ) & g\cos \theta _{0}\\ \left ( \frac{L_{v}}{I_{x}^{\prime }}+I_{zx}^{\prime }N_{v}\right ) & \left ( \frac{L_{p}}{I_{x}^{\prime }}+I_{zx}^{\prime }N_{p}\right ) & \left ( \frac{L_{r}}{I_{x}^{\prime }}+I_{zx}^{\prime }N_{r}\right ) & 0\\ \left ( I_{zx}^{\prime }L_{v}+\frac{N_{v}}{I_{z}^{\prime }}\right ) & \left ( I_{zx}^{\prime }L_{p}+\frac{N_{p}}{I_{z}^{\prime }}\right ) & \left ( I_{zx}^{\prime }L_{r}+\frac{N_{r}}{I_{z}^{\prime }}\right ) & 0\\ 0 & 1 & \tan \theta _{0} & 0 \end{pmatrix}\begin{pmatrix} v\\ p\\ r\\ \phi \end{pmatrix} \tag{4.9,19} \end{equation}
To find air density \rho , we are told the airplane is at see level. Hence using table in appendix D then we find the corresponding air density at that altitude\begin{align*} \rho & =2.3769\times 10^{-3}\text{ lb }\frac{\text{sec}^{2}}{ft^{4}}\\ & =1.225\text{ kg/m}^{3} \end{align*}
In this case C_{L}=\frac{L}{\frac{1}{2}\rho V^{2}S}
| non-dim. | value |
| C_{y_{\beta }} | -0.14 |
| C_{y_{p}} | -0.039 |
| C_{y_{r}} | 0.165 |
| C_{l_{\beta }} | -0.0689-0.0917C_{L}=-0.0689-0.0917\times 0.195\,76= -0.08685\, |
| C_{l_{p}} | -0.441 |
| C_{l_{r}} | -0.0144+0.271C_{L}=-0.0144+0.271\times 0.195\,76= 0.03865 |
| C_{n_{\beta }} | 0.01326+0.017C_{L}^{2}=0.01326+0.017\times 0.195\,76^{2}= 0.0139 |
| C_{n_{p}} | -0.00109-0.0966C_{L}=-0.00109-0.0966\times 0.195\,76= -0.02 |
| C_{n_{r}} | -0.048-0.0238C_{L}^{2}=-0.048-0.0238\times 0.195\,76^{2}= -0.0489 |
Hence
| dim. | non-dim. | numerical equation | value |
| Y_{v} | \frac{1}{2}\rho u_{0}SC_{y_{\beta }} | \frac{1}{2}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 14.9\right ) \left ( -0.14\right ) | -98.764 |
| Y_{p} | \frac{1}{4}\rho u_{0}bSC_{y_{p}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14\right ) \left ( 14.9\right ) \left ( -0.039\right ) | -125.73 |
| Y_{r} | \frac{1}{4}\rho u_{0}bSC_{y_{r}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14\right ) \left ( 14.9\right ) \left ( 0.165\right ) | 531.95 |
| L_{v} | \frac{1}{2}\rho u_{0}SC_{l_{\beta }} | \frac{1}{2}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14\right ) \left ( 14.9\right ) \left ( -8.685\,1\times 10^{-2}\right ) | -560.01 |
| L_{p} | \frac{1}{4}\rho u_{0}b^{2}SC_{l_{p}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14^{2}\right ) \left ( 14.9\right ) \left ( -0.441\right ) | -12995 |
| L_{r} | \frac{1}{4}\rho u_{0}b^{2}SC_{l_{r}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14^{2}\right ) \left ( 14.9\right ) \left ( 3.865\,1\times 10^{-2}\right ) | 1138.9 |
| N_{v} | \frac{1}{2}\rho u_{0}bSC_{n_{\beta }} | \frac{1}{2}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14\right ) \left ( 14.9\right ) \left ( 1.391\,1\times 10^{-2}\right ) | 89.700 |
| N_{p} | \frac{1}{4}\rho u_{0}b^{2}SC_{n_{p}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14^{2}\right ) \left ( 14.9\right ) \left ( -0.02\right ) | -589.35 |
| N_{r} | \frac{1}{4}\rho u_{0}b^{2}SC_{n_{r}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14^{2}\right ) \left ( 14.9\right ) \left ( -4.891\,2\times 10^{-2}\right ) | -1441.3 |
Now that all the numerical values are calculated, 4.19,9 is written again in order to find a numerical A matrix in order to determine its eigenvalues.
Since \theta =0, equation 4.19,9 becomes\begin{pmatrix} \dot{v}\\ \dot{p}\\ \dot{r}\\ \dot{\phi }\end{pmatrix} =\begin{pmatrix} \frac{-98.764}{10675/9.81} & \frac{-125.73}{10675/9.81} & \left ( \frac{531.95}{10675/9.81}-77.3\right ) & 9.81\\ \left ( \frac{-560.01}{230}+I_{zx}^{\prime }\left ( 89.7\right ) \right ) & \left ( \frac{-12995}{230}+I_{zx}^{\prime }\left ( -589.35\right ) \right ) & \left ( \frac{1138.9}{230}+I_{zx}^{\prime }\left ( -1441.3\right ) \right ) & 0\\ \left ( I_{zx}^{\prime }\left ( -560.01\right ) +\frac{89.700}{1778}\right ) & \left ( I_{zx}^{\prime }\left ( -12995\right ) +\frac{-589.35}{1778}\right ) & \left ( I_{zx}^{\prime }\left ( 1138.9\right ) +\frac{-1441.3}{1778}\right ) & 0\\ 0 & 1 & 0 & 0 \end{pmatrix}\begin{pmatrix} v\\ p\\ r\\ \phi \end{pmatrix}
\begin{align*} \lambda _{dutch} & =-0.4218\pm 2.2873i\\ \lambda _{spiral} & =-0.0551\\ \lambda _{rolling} & =-56.5025 \end{align*}
The characteristic polynomial is
Hence\begin{align*} p\left ( \lambda \right ) & =\lambda ^{4}+57.4014\lambda ^{3}+56.2385\lambda ^{2}+308.576\lambda +16.832\\ & =A\lambda ^{4}+B\lambda ^{3}+C\lambda ^{2}+D\lambda +E \end{align*}
Therefore\begin{align*} E & =16.832>0\\ R & =D\left ( BC-AD\right ) -B^{2}E\\ & =308.576\left ( 57.4014\times 56.2385-1\times 308.576\right ) -\left ( 57.4014^{2}\right ) \left ( 16.832\right ) \\ & =8.454\,6\times 10^{5}>0 \end{align*}
Hence all modes are stable since both E and R are positive. Now that the eigenvalues are known, the system characteristic timing table is generated, as was done in the textbook on page 188. Here is the table for this airplane
| Mode | name | \lambda =n\pm \omega i | period (sec) | t_{half}\left ( s\right ) | N_{half} (cycles) |
| \frac{2\pi }{\omega } | \frac{0.693}{\left \vert n\right \vert } | 0.11\frac{\omega }{\left \vert n\right \vert } | |||
| 1 | spiral | -0.0551 | — | \frac{0.693}{0.0551}=12.577 | — |
| 2 | rolling | -56.5025 | — | \frac{0.693}{56.5025}=0.012265 | — |
| 3 | dutch roll | -0.4218\pm 2.2873i | \frac{2\pi }{2.2873}= 2.747\,0 | \frac{0.693}{0.4218}=1.643\,0 | 0.11\frac{2.2873}{0.4218}=0.596\,50 |
The dutch roll is oscillatory, its characteristic transients is plotted below
The corresponding eigenvectors are now found in order to generate the eigenvector phase diagrams similar to figure 6.3, page 169 in the textbook. The solution will appear as follows\begin{pmatrix} v\\ p\\ r\\ \phi \end{pmatrix} =\begin{pmatrix} x_{11}\\ x_{21}\\ x_{31}\\ x_{41}\end{pmatrix} e^{\lambda _{1}t}+\begin{pmatrix} x_{12}\\ x_{22}\\ x_{32}\\ x_{42}\end{pmatrix} e^{\lambda _{2}t}+\begin{pmatrix} x_{13}\\ x_{23}\\ x_{33}\\ x_{43}\end{pmatrix} e^{\lambda _{3}t}+\begin{pmatrix} x_{14}\\ x_{24}\\ x_{34}\\ x_{44}\end{pmatrix} e^{\lambda _{4}t}
In this problem, Matlab was used to obtain the eigenvectors.
The problem asks to normalize the eigenvector using the third entry, which is r. Therefore, after finding the eigenvectors using the eig command, each entry in the eigenvector was divided by the third entry in the same eigenvector. A small function was written to automate this process for both C_{L}=0 and C_{L}\neq 0. The function and its full output are listed in the appendix of this problem. Here is the result found. The entries in the eigenvector were made to be non-dimensional. This seems to be what was done in the textbook, page 189. Non-dimensional eigenvector was generated from the eigenvector returned by Matlab as follows:
To convert eigenvector \begin{pmatrix} v\\ p\\ r\\ \Delta \phi \end{pmatrix} to non-dimensional form, we multiply elements as follows \begin{pmatrix} \hat{v}\\ \hat{p}\\ \hat{r}\\ \Delta \hat{\phi }\end{pmatrix} =\begin{pmatrix} \frac{v}{u_{0}}\\ p\frac{b}{2u_{0}}\\ r\frac{b}{2u_{0}}\\ \Delta \phi \end{pmatrix} . This was done for each eigenvector of each mode. The dutch mode has two complex conjugate eigenvalues and counts as one mode.
After the eigenvectors are found, polar form table for the eigenvectors is made, similar to table 6.4, page 168 of the textbook, then the vector phasor diagrams is drawn for the dutch mode.
The polar form of each eigenvector is summarized below
| dutch \lambda _{3,4}=-0.4218\pm 2.2873i
| spiral \lambda _{1}=-0.055268
| rolling \lambda _{2}=-56.5025
| ||||
| magnitude | phase (deg) | magnitude | phase (deg) | magnitude | phase (deg) | |
| v | 35.874 | 80.158 | 12.115 | 0 | 2.2246 | 0 |
| p | 1.5436 | -98.95 | 0.43488 | 180 | 168.36 | 0 |
| r | 1 | 0 | 1 | 0 | 1 | 0 |
| \Delta \phi | 0.66332 | -199.39 | 7.868 | 0 | 2.9796 | -180 |
non-dimensional result v,p,r are first made non-dimensional. Next, the ratio variable to r is found. The Matlab function below shows the the implementation details.
| dutch \lambda _{3,4}=-0.4218\pm 2.2873i
| spiral \lambda _{1}=-0.055268
| rolling \lambda _{2}=-56.5025
| ||||
| amplitude | phase (deg) | amplitude | phase (deg) | amplitude | phase (deg) | |
| v/r | 7.849 | 80.158 | 2.651 | 0 | 0.4867 | 0 |
| p/r | 1.543 | -98.95 | 0.4348 | 180 | 168.36 | 0 |
| r/r | 1 | 0 | 1 | 0 | 1 | 0 |
| \frac{\Delta \phi }{r} | 11.22 | -199.39 | 133.09 | 0 | 50.4 | -180 |
| norm | 13.33 | 133.09 | 176.43 | |||
Figure 2.28 is the eigenvector diagram for the dutch mode
The above calculations are now repeated for C_{L}=0 case. Calculation of the numerical values of the dimensional derivatives gives
| non-dim. | value |
| C_{y_{\beta }} | -0.14 |
| C_{y_{p}} | -0.039 |
| C_{y_{r}} | 0.165 |
| C_{l_{\beta }} | -0.0689-0.0917C_{L}=-0.0689 |
| C_{l_{p}} | -0.441 |
| C_{l_{r}} | -0.0144+0.271C_{L}=-0.0144 |
| C_{n_{\beta }} | 0.01326+0.017C_{L}^{2}=0.01326 |
| C_{n_{p}} | -0.00109-0.0966C_{L}=-0.00109 |
| C_{n_{r}} | -0.048-0.0238C_{L}^{2}=-0.048 |
Hence
| dim. | non-dim. | numerical equation | value |
| Y_{v} | \frac{1}{2}\rho u_{0}SC_{y_{\beta }} | \frac{1}{2}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 14.9\right ) \left ( -0.14\right ) | -98.764 |
| Y_{p} | \frac{1}{4}\rho u_{0}bSC_{y_{p}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14\right ) \left ( 14.9\right ) \left ( -0.039\right ) | -125.73 |
| Y_{r} | \frac{1}{4}\rho u_{0}bSC_{y_{r}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14\right ) \left ( 14.9\right ) \left ( 0.165\right ) | 531.95 |
| L_{v} | \frac{1}{2}\rho u_{0}SC_{l_{\beta }} | \frac{1}{2}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14\right ) \left ( 14.9\right ) \left ( -0.0689\right ) | -444.26 |
| L_{p} | \frac{1}{4}\rho u_{0}b^{2}SC_{l_{p}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14^{2}\right ) \left ( 14.9\right ) \left ( -0.441\right ) | -12995 |
| L_{r} | \frac{1}{4}\rho u_{0}b^{2}SC_{l_{r}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14^{2}\right ) \left ( 14.9\right ) \left ( -0.0144\right ) | -424.32 |
| N_{v} | \frac{1}{2}\rho u_{0}bSC_{n_{\beta }} | \frac{1}{2}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14\right ) \left ( 14.9\right ) \left ( 0.01326\right ) | 85.499 |
| N_{p} | \frac{1}{4}\rho u_{0}b^{2}SC_{n_{p}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14^{2}\right ) \left ( 14.9\right ) \left ( -0.00109\right ) | -32.119 |
| N_{r} | \frac{1}{4}\rho u_{0}b^{2}SC_{n_{r}} | \frac{1}{4}\left ( 1.225\right ) \left ( 77.3\right ) \left ( 9.14^{2}\right ) \left ( 14.9\right ) \left ( -0.048\right ) | -1441.3 |
With the above, equation 4.19,9 is written again in order to find a numerical A matrix to use to find its eigenvalues, and since \theta =0 then 4.19,9 becomes (but remember to use a 3\times 3 matrix in this case, since the 4^{th} column is now zero column)\begin{align*} \begin{pmatrix} \dot{v}\\ \dot{p}\\ \dot{r}\end{pmatrix} & =\begin{pmatrix} \frac{Y_{v}}{m} & \frac{Y_{p}}{m} & \left ( \frac{Y_{r}}{m}-u_{0}\right ) \\ \left ( \frac{L_{v}}{I_{x}^{\prime }}+I_{zx}^{\prime }N_{v}\right ) & \left ( \frac{L_{p}}{I_{x}^{\prime }}+I_{zx}^{\prime }N_{p}\right ) & \left ( \frac{L_{r}}{I_{x}^{\prime }}+I_{zx}^{\prime }N_{r}\right ) \\ \left ( I_{zx}^{\prime }L_{v}+\frac{N_{v}}{I_{z}^{\prime }}\right ) & \left ( I_{zx}^{\prime }L_{p}+\frac{N_{p}}{I_{z}^{\prime }}\right ) & \left ( I_{zx}^{\prime }L_{r}+\frac{N_{r}}{I_{z}^{\prime }}\right ) \end{pmatrix} \\ & =\begin{pmatrix} \frac{-98.764}{10675/9.81} & \frac{-125.73}{10675/9.81} & \left ( \frac{531.95}{10675/9.81}-77.3\right ) \\ \left ( \frac{-444.26}{230}+I_{zx}^{\prime }\left ( 85.499\right ) \right ) & \left ( \frac{-12995}{230}+I_{zx}^{\prime }\left ( -32.119\right ) \right ) & \left ( \frac{-424.32}{230}+I_{zx}^{\prime }\left ( -1441.3\right ) \right ) \\ \left ( I_{zx}^{\prime }\left ( -444.26\right ) +\frac{85.499}{1778}\right ) & \left ( I_{zx}^{\prime }\left ( -12995\right ) +\frac{-32.119}{1778}\right ) & \left ( I_{zx}^{\prime }\left ( -424.32\right ) +\frac{-1441.3}{1778}\right ) \end{pmatrix}\begin{pmatrix} v\\ p\\ r \end{pmatrix} \end{align*}
And since I_{zx}^{\prime }=0 the above reduces to\begin{pmatrix} \dot{v}\\ \dot{p}\\ \dot{r}\end{pmatrix}\begin{pmatrix} -0.09076\,1 & -0.115\,54 & -76.811\\ -1.931\,6 & -56.5 & -1.844\,9\\ 0.04808\,7 & -0.01806\,5 & -0.79551 \end{pmatrix}\begin{pmatrix} v\\ p\\ r \end{pmatrix}
The characteristic polynomial is
Hence\begin{align*} p\left ( \lambda \right ) & =0x^{4}+x^{3}+57.3858x^{2}+53.5831x+215.257\\ & =A\lambda ^{4}+B\lambda ^{3}+C\lambda ^{2}+D\lambda +E \end{align*}
Therefore\begin{align*} E & =215.257>0\\ R & =D\left ( BC-AD\right ) -B^{2}E\\ & =53.5831\left ( 1\times 57.3858\right ) -\left ( 1^{2}\right ) \left ( 215.257\right ) \\ & =2859.7>0 \end{align*}
All modes are stable since both E and R are positive. Now that the eigenvalues are known, the characteristic times table is generated, as was done in the textbook on page 188.
| Mode | name | \lambda =n\pm \omega i | period (sec) | t_{half}\left ( s\right ) | N_{half} (cycles) |
| \frac{2\pi }{\omega } | \frac{0.693}{\left \vert n\right \vert } | 0.11\frac{\omega }{\left \vert n\right \vert } | |||
| 1 | rolling | -56.505 | — | \frac{0.693}{56.505}=0.01226 | — |
| 2 | dutch roll | -0.44044\pm 1.9015i | \frac{2\pi }{1.9015}=3.304\,3 | \frac{0.693}{0.44044}=1.573\,4 | 0.11\frac{1.9015}{0.44044}=0.474\,9 |
The dutch roll is oscillatory. The characteristic transient is plotted below
Comparing the above to g\neq 0, it is seen that there is very little change in the dutch mode when g=0.
Similar to first part, now that all the eigenvectors are found, the polar form table for the eigenvectors is made which is similar to table 6.4, page 168 of the textbook. This is followed by making the vector phasor diagrams for the dutch mode.
The polar form of each eigenvector is summarized below
| dutch \lambda _{1,2}=-0.44044\pm 1.9015i
| rolling \lambda _{3}=-56.505
| |||
| amplitude | phase (deg) | amplitude | phase (deg) | |
| v | 39.71 | 79.465 | 7.71 | 0 |
| p | 1.374 | 256.17 | 3104.4 | 0 |
| r | 1 | 0 | 1 | 0 |
| dutch \lambda _{1,2}=-0.44044\pm 1.9015i
| rolling \lambda _{2}=-56.505
| |||
| amplitude | phase (deg) | amplitude | phase (deg) | |
| \frac{v}{r} | 8.6893 | 79.465 | 1.69 | 0 |
| \frac{p}{r} | 1.3739 | 256.17 | 3104.4 | 0 |
| \frac{r}{r} | 1 | 0 | 1 | 0 |
| norm | 9.2688 | 3105.4 | ||
Figure 2.30 shows the eigenvector diagram for the dutch mode for g=0
The first observation is that, when gravity is absent, the spiral mode vanishes. There are only three modes when g=0, and these are the dutch (complex conjugate eigenvalues) and rolling (real eigenvalue) and yaw r. Only the dutch and roll modes can be compared to each others for both cases. Spiral mode can not be compared since this mode does not exist when g=0.
For the dutch mode (the oscillatory mode), There is no significant change that can be seen. So one can conclude that gravity has little effect on the dutch mode.
The final solution is now written for both cases, and plotted against each others to compare graphically. (this uses the ratios in the eigenvectors components taken from the non-dimensional results obtained above). The solution for g\neq 0 is
\begin{align*} \begin{pmatrix} \frac{v}{r}\\ \frac{p}{r}\\ \frac{r}{r}\\ \phi _{r}\end{pmatrix} & =\begin{pmatrix} x_{11}\\ x_{21}\\ x_{31}\\ x_{41}\end{pmatrix} _{spiral}e^{\lambda _{1}t}+\begin{pmatrix} x_{12}\\ x_{22}\\ x_{32}\\ x_{42}\end{pmatrix} _{rolling}e^{\lambda _{2}t}+\begin{pmatrix} x_{13}\\ x_{23}\\ x_{33}\\ x_{43}\end{pmatrix} _{dutch}e^{\lambda _{3}t}+\begin{pmatrix} x_{14}\\ x_{24}\\ x_{34}\\ x_{44}\end{pmatrix} _{dutch}e^{\lambda _{4}t}\\ & =\begin{pmatrix} 2.651\\ -0.4348\\ 1\\ 133.09 \end{pmatrix} e^{-0.055268t}+\begin{pmatrix} 0.4867\\ 168.36\\ 1\\ -50.4 \end{pmatrix} e^{-56.5025t}+\\ & e^{-0.4218t}\begin{pmatrix} 7.849\left ( \cos \left ( 2.2873t+80.15^{o}\right ) +\sin \left ( 2.2873t+80.15^{o}\right ) \right ) \\ 1.543\left ( \cos \left ( 2.2873t-98.95^{o}\right ) +\sin \left ( 2.2873t-98.95^{o}\right ) \right ) \\ \left ( \cos 2.2873t+\sin 2.2873t\right ) \\ 11.22\left ( \cos \left ( 2.2873t-199.39^{o}\right ) +\sin \left ( 2.2873t-199.39^{o}\right ) \right ) \end{pmatrix} \end{align*}
Hence
\begin{align*} \frac{v}{r} & =2.651e^{-0.055268t}+0.4867e^{-56.5025t}+7.849e^{-0.4218t}\left ( \cos \left ( 2.2873t+80.15^{o}\right ) +\sin \left ( 2.2873t+80.15^{o}\right ) \right ) \\ \frac{p}{r} & =-0.4348e^{-0.055268t}+\overbrace{168.36e^{-56.5025t}}^{\text{compare to g=0 below}} +1.543e^{-0.4218t}\left ( \cos \left ( 2.2873t-98.95^{o}\right ) +\sin \left ( 2.2873t-98.95^{o}\right ) \right ) \\ \frac{r}{r} & =e^{-0.055268t}+e^{-56.5025t}+e^{-0.4218t}\left ( \cos 2.2873t+\sin 2.2873t\right ) \\ \frac{\phi }{r} & =133.09e^{-0.055268t}-50.4e^{-56.5025t}+11.22e^{-0.4218t}\left ( \cos \left ( 2.2873t-199.39^{o}\right ) +\sin \left ( 2.2873t-199.39^{o}\right ) \right ) \end{align*}
The solution for g=0 is
\begin{align*} \begin{pmatrix} \frac{v}{r}\\ \frac{p}{r}\\ \frac{r}{r}\end{pmatrix} & =\begin{pmatrix} x_{12}\\ x_{22}\\ x_{32}\end{pmatrix} _{rolling}e^{\lambda _{2}t}+\begin{pmatrix} x_{13}\\ x_{23}\\ x_{33}\end{pmatrix} _{dutch}e^{\lambda _{3}t}+\begin{pmatrix} x_{14}\\ x_{24}\\ x_{34}\end{pmatrix} _{dutch}e^{\lambda _{4}t}\\ & =\begin{pmatrix} 1.69\\ 3104.4\\ 1 \end{pmatrix} e^{-56.505t}+e^{-0.44044t}\begin{pmatrix} 8.6893\left ( \cos \left ( 1.9015t+79.465^{o}\right ) +\sin \left ( 1.9015t+79.465^{o}\right ) \right ) \\ 1.3739\left ( \cos \left ( 1.9015t+256.17^{o}\right ) +\sin \left ( 1.9015t+256.17^{o}\right ) \right ) \\ \left ( \cos 1.9015t+\sin 1.9015t\right ) \end{pmatrix} \end{align*}
Hence
\begin{align*} \frac{v}{r} & =1.69e^{-56.505t}+8.68939e^{-0.44044t}\left ( \cos \left ( 1.9015t+79.465^{o}\right ) +\sin \left ( 1.9015t+79.465^{o}\right ) \right ) \\ \frac{p}{r} & =\overbrace{3104.4e^{-56.505t}}^{\text{larger but quickly damps}} +1.3739e^{-0.44044t}\left ( \cos \left ( 1.9015t\right ) +\sin \left ( 1.9015t\right ) \right ) \\ \frac{r}{r} & =e^{-56.505t}+e^{-0.44044t}\left ( \cos 1.9015t+\sin 1.9015t\right ) \end{align*}
From above, we see when g=0 the contribution from rolling mode has much larger amplitude (3104 vs. 168.36). But this damps out very fast thanks to the large negative exponent on the exponential.
To compare the different modes for gravity and no gravity, each component from each solution is plotted against the similar component from the other solution. For example, \left ( \frac{v}{r}\right ) _{g=0} and \left ( \frac{v}{r}\right ) _{g\neq 0} are plotted on same plot. The same for \left ( \frac{p}{r}\right ) _{g=0} and \left ( \frac{p}{r}\right ) _{g\neq 0}. Figure 2.31 shows the the result.
We see that there is little change in this part of the solution. Figure 2.32 shows the p/r result
Looking Figure 2.32, there does not seem to be difference as well. But that is because \frac{p}{r} damped very quickly, even though it was much larger in the case when there is no gravity. We can see this more clearly if we zoom to a smaller time span, say for t=0.1 seconds only. Figure 2.33 shows the result
This means that the rolling eigen mode will contribute much more to the p(t) component of the total solution with no gravity when compared to with gravity. In other words, rolling becomes the more dominant motion with no gravity. But this only affect the solution for short amount of time. It will depend on other conditions and on the system being considered if this is important or not.
This also seems to imply that rolling control becomes more important in outer space since gravity becomes very weak. We do not need to worry about spiral mode control since it does not exist. Dutch mode is not affected by gravity.
solution
We are told in hint to use x= \begin{pmatrix} \theta \\ q\\ \phi \\ p \end{pmatrix} as the state vector. Therefore we need to set up a state equation \dot{x}=Ax that has this form
\begin{equation} \begin{pmatrix} \dot{\theta }\\ \dot{q}\\ \dot{\phi }\\ \dot{p}\end{pmatrix} =A\begin{pmatrix} \theta \\ q\\ \phi \\ p \end{pmatrix} \tag{1} \end{equation}
So the question is, how to build the A matrix above? We have solved this same problem in HW 2, and in there we found expressions for \left \{ \dot{\theta },\dot{q},\dot{\phi },\dot{p}\right \} and these will be used as is in this problem. These are the equations found in HW2 for this problem
\begin{align*} \dot{\theta } & =q\\ \dot{q} & =\frac{M-pH}{I_{y}}\\ \dot{\phi } & =p\\ \dot{p} & =\frac{L+qH}{I_{x}} \end{align*}
Where, from problem 4.10 statement, M=M_{\theta }\theta and L=L_{\phi }\phi , hence the above becomes
\begin{align*} \dot{\theta } & =q\\ \dot{q} & =\frac{M_{\theta }}{I_{y}}\theta -\frac{H}{I_{y}}p\\ \dot{\phi } & =p\\ \dot{p} & =\frac{L_{\phi }}{I_{x}}\phi +q\frac{H}{I_{x}} \end{align*}
Substituting the above in (1) results in
\begin{equation} \begin{pmatrix} \dot{\theta }\\ \dot{q}\\ \dot{\phi }\\ \dot{p}\end{pmatrix} =\begin{pmatrix} 0 & 1 & 0 & 0\\ \frac{M_{\theta }}{I_{y}} & 0 & 0 & -\frac{H}{I_{y}}\\ 0 & 0 & 0 & 1\\ 0 & \frac{H}{I_{x}} & \frac{L_{\phi }}{I_{x}} & 0 \end{pmatrix}\begin{pmatrix} \theta \\ q\\ \phi \\ p \end{pmatrix} \tag{2} \end{equation}
The characteristic equation can now be easily found from the definition
\begin{align*} \left \vert A-\lambda I\right \vert & =0\\\begin{vmatrix} -\lambda & 1 & 0 & 0\\ \frac{M_{\theta }}{I_{y}} & -\lambda & 0 & -\frac{H}{I_{y}}\\ 0 & 0 & -\lambda & 1\\ 0 & \frac{H}{I_{x}} & \frac{L_{\phi }}{I_{x}} & -\lambda \end{vmatrix} & =0 \end{align*}
The determinant is found using CAS and this is the result
The characteristic equation is found
Therefore
\boxed{ p\left ( \lambda \right ) =\lambda ^{4}+\frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\lambda ^{2}+\frac{LM}{I_{x}I_{y}}}
\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2} is always positive. Let A=\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) , hence the above can be written as\begin{align*} s^{2} & =-\frac{1}{2}A\pm \frac{1}{2}\sqrt{A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }\\ & =\left \{ \overbrace{-\frac{1}{2}A+\frac{1}{2}\sqrt{A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }}^{\text{case 1}} , \overbrace{-\frac{1}{2}A-\frac{1}{2}\sqrt{A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }}^{\text{case 2}} \right \} \end{align*}
We want to find out if the real part of s is negative or not to decide on stability. Considering each case at a time.
s^{2}=-\frac{1}{2}A+\frac{1}{2}\sqrt{A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }
If B<0 then s^{2}=-\frac{1}{2}A+\frac{1}{2}i\sqrt{\left \vert B\right \vert }
The real part of +\sqrt{z}e^{\frac{i}{2}\arg z} is \operatorname{Re}\left ( \sqrt{z}\left [ \cos \left ( \frac{\arg z}{2}\right ) +i\sin \left ( \frac{\arg z}{2}\right ) \right ] \right ) =\sqrt{z}\cos \left ( \frac{\arg z}{2}\right ) and the real part of the second root is \operatorname{Re}\left ( -\sqrt{z}\left [ \cos \left ( \frac{\arg z}{2}\right ) +i\sin \left ( \frac{\arg z}{2}\right ) \right ] \right ) =-\sqrt{\left \vert z\right \vert }\cos \left ( \frac{\arg z}{2}\right )
These two complex roots have the following real parts\begin{align*} \operatorname{Re}\left ( s_{1}\right ) & =-\sqrt{\left \vert z\right \vert }\cos \left ( \frac{\arg z}{2}\right ) \\ \operatorname{Re}\left ( s_{2}\right ) & =\sqrt{\left \vert z\right \vert }\cos \left ( \frac{\arg z}{2}\right ) \end{align*}
These real parts can not be both negative or both positive at the same time. If we assume \cos \left ( \frac{\arg z}{2}\right ) >0 for some \frac{\arg z}{2}, then \operatorname{Re}\left ( s_{1}\right ) <0 (since \sqrt{\left \vert z\right \vert }>0 all the time, since \left \vert z\right \vert is the magnitude of the complex root). But \operatorname{Re}\left ( s_{2}\right ) >0 or unstable.
On the other hand, if we assume \cos \left ( \frac{\arg z}{2}\right ) <0 then \operatorname{Re} \left ( s_{2}\right ) <0 but now \operatorname{Re}\left ( s_{1}\right ) >0, hence unstable. Therefore if B<0 then the system is not stable as one of the roots must have positive real part. Since B=A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) , this implies that 4\left ( \frac{LM}{I_{x}I_{y}}\right ) >A^{2}. Since A^{2}=\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2} then this condition means \begin{align*} 4\left ( \frac{LM}{I_{x}I_{y}}\right ) & >\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}\\ 4LM & >\frac{\left ( H^{2}-LI_{y}-MI_{x}\right ) ^{2}}{I_{x}I_{y}}\text{\ \ (unstable)} \end{align*}
Therefore, if B<0, the system is always unstable as one of the roots is unstable.
If B>0 then\begin{align*} s^{2} & =-\frac{1}{2}A+\frac{1}{2}\sqrt{\left \vert B\right \vert }\\ & =\frac{1}{2}\left ( \sqrt{\left \vert B\right \vert }-A\right ) \end{align*}
Now \sqrt{\left \vert B\right \vert } is positive number. Let \sqrt{\left \vert B\right \vert }=\left \vert D\right \vert then s^{2}=\frac{1}{2}\left ( \left \vert D\right \vert -A\right )
If \left ( \left \vert D\right \vert -A\right ) >0 then s_{1} is real and positive, hence unstable. This means \left \vert D\right \vert >A or \sqrt{A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }>A and since A=\frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}} then the condition is\begin{align*} \sqrt{\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) } & >\frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\\ \left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) & >\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}\\ -4\left ( \frac{LM}{I_{x}I_{y}}\right ) & >0\\ \frac{LM}{I_{x}I_{y}} & <0\\ LM & <0\text{ \ \ \ (unstable)} \end{align*}
But if \left ( \left \vert D\right \vert -A\right ) <0 then s_{1} and s_{2} are pure imaginary numbers. Hence the system is marginally stable (real part is zero. Some books call this marginally unstable). This means \left \vert D\right \vert <A or \sqrt{A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }<A then the condition is \begin{align*} \sqrt{\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) } & <\frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\\ \left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) & <\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}\\ -4\left ( \frac{LM}{I_{x}I_{y}}\right ) & <0\\ LM & >0\text{ \ (marginally stable)} \end{align*}
We need repeat all the above for case 2
s^{2}=-\frac{1}{2}A-\frac{1}{2}\sqrt{A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }
If B<0 then s^{2}=-\frac{1}{2}A-\frac{1}{2}i\sqrt{\left \vert B\right \vert }
Therefore, if B<0, system unstable.
If B>0 then\begin{align*} s^{2} & =-\frac{1}{2}A-\frac{1}{2}\sqrt{\left \vert B\right \vert }\\ & =-\frac{1}{2}\left ( \sqrt{\left \vert B\right \vert }+A\right ) \end{align*}
Now \sqrt{\left \vert B\right \vert } is positive number, since B>0. Then\begin{align*} s & =\pm i\sqrt{\frac{1}{2}}\sqrt{\sqrt{\left \vert B\right \vert }+A}\\ & =\left \{ i\sqrt{\frac{1}{2}}\sqrt{\sqrt{\left \vert B\right \vert }+A},-i\sqrt{\frac{1}{2}}\sqrt{\sqrt{\left \vert B\right \vert }+A}\right \} \end{align*}
If \sqrt{\left \vert B\right \vert }+A>0 then both s_{1} and s_{2} are pure imaginary and the system is marginally stable (real part is zero). This condition is \sqrt{\left \vert B\right \vert }+\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) >0 or \sqrt{\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }+\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) >0\text{\ \ \ \ \ (marginally stable)}
If \sqrt{\left \vert B\right \vert }+A<0 then s_{1}=i\sqrt{\frac{1}{2}}i\left \vert \sqrt{\sqrt{\left \vert B\right \vert }+A}\right \vert =-\sqrt{\frac{1}{2}}\left \vert \sqrt{\sqrt{\left \vert B\right \vert }+A}\right \vert which is negative, hence s_{1} is stable. and s_{2}=-i\sqrt{\frac{1}{2}}i\left \vert \sqrt{\sqrt{\left \vert B\right \vert }+A}\right \vert =\sqrt{\frac{1}{2}}\left \vert \sqrt{\sqrt{\left \vert B\right \vert }+A}\right \vert which is positive, hence s_{2} is not stable.
This condition is \sqrt{\left \vert B\right \vert }+\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) <0 or \sqrt{\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }+\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) <0\text{\ \ \ \ (not stable)}
This covers all 4 cases. We now summarize the findings
| condition | meaning |
| LM>\frac{1}{4}\frac{\left ( H^{2}-LI_{y}-MI_{x}\right ) ^{2}}{I_{x}I_{y}}\text{\ \ } | unstable |
| LM<0\text{ \ \ \ } | \text{unstable} |
| LM>0 | marginally stable |
| \sqrt{\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }+\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) >0 | marginally stable |
| \sqrt{\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }+\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) <0 | not stable |
I do not think I did this correctly. LM>\frac{1}{4}\frac{\left ( H^{2}-LI_{y}-MI_{x}\right ) ^{2}}{I_{x}I_{y}} comes out as unstable. But \frac{1}{4}\frac{\left ( H^{2}-LI_{y}-MI_{x}\right ) ^{2}}{I_{x}I_{y}} is positive quantity, while LM>0 came out as marginally stable. We can’t have both cases be true.