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 \] Where \(\gamma \) is the angle between horizontal plane and the velocity vector \(V\). This is illustrated in figure 2.25 from our text book on page 19.
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} \] Taking the \(\lim _{m\rightarrow \infty }\) results in\[\begin{pmatrix} \dot{w}\\ \dot{q}\end{pmatrix} =\begin{pmatrix} 0 & u_{o}\\ \frac{1}{I_{y}}M_{w} & \frac{1}{I_{y}}\left ( M_{q}+M_{\dot{w}}u_{o}\right ) \end{pmatrix}\begin{pmatrix} w\\ q \end{pmatrix} \] The characteristic equation is found from\begin{align} \det \left ( A-\lambda I\right ) & =0\nonumber \\\begin{vmatrix} -\lambda & u_{o}\\ \frac{1}{I_{y}}M_{w} & \frac{1}{I_{y}}\left ( M_{q}+M_{\dot{w}}u_{o}\right ) -\lambda \end{vmatrix} & =0\nonumber \\ \lambda ^{2}-\lambda \frac{1}{I_{y}}\left ( M_{q}+M_{\dot{w}}u_{o}\right ) -\frac{1}{I_{y}}M_{w}u_{o} & =0 \tag{1} \end{align}
Let the natural frequency\[ \boxed{ \omega _{n}=\sqrt{-M_{w}\frac{u_{o}}{I_{y}}} } \] And the damping ratio\[ \boxed{ \zeta =\frac{-\left ( M_{q}+M_{\dot{w}}u_{o}\right ) }{I_{y}\sqrt{-M_{w}\frac{u_{o}}{I_{y}}}} } \] Equation (1) can be written in standard form as\[ \lambda ^{2}+2\zeta \omega _{n}\lambda +\omega _{n}^{2}=0 \] The eigenvalues are\begin{align*} \lambda & =\frac{-b\pm \sqrt{b^{2}-4ac}}{2a}\\ & =\frac{1}{2I_{y}}\left ( M_{q}+M_{\dot{w}}u_{o}\right ) \pm \frac{1}{2I_{y}}\sqrt{\left ( M_{q}+M_{\dot{w}}u_{o}\right ) ^{2}+4I_{y}M_{w}u_{o}}\\ & =\frac{1}{2I_{y}}\left ( M_{q}+M_{\dot{w}}u_{o}\right ) \pm \frac{i}{2I_{y}}\sqrt{-\left ( M_{q}+M_{\dot{w}}u_{o}\right ) ^{2}-4I_{y}M_{w}u_{o}}\\ & =n\pm i\omega \end{align*}
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}}} } \] The time to double is\begin{align*} t & =\frac{0.693}{\left \vert n\right \vert }\\ & =\frac{\left ( 2I_{y}\right ) 0.693}{\left \vert M_{q}+M_{\dot{w}}u_{o}\right \vert }\\ & = I_{y}\frac{1.386}{\left \vert M_{q}+M_{\dot{w}}u_{o}\right \vert } \end{align*}
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 ] \] Using table 7.2, each of the above expressions are evaluated. Some are functions of \(C_{L}\). \begin{align*} C_{l_{\beta }} & =-0.0689-0.0917C_{L}\\ C_{n_{r}} & =-0.048-0.0238C_{L}^{2}\\ C_{l_{r}} & =-0.0144+0.271C_{L}\\ C_{n_{\beta }} & =0.01326+0.017C_{L}^{2}\\ C_{l_{p}} & =-0.441\\ C_{n_{p}} & =-0.00109-0.0966C_{L} \end{align*}
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}\] At trim, \(L=mg\), hence\[ C_{L}=\frac{mg}{\frac{1}{2}\rho V^{2}S}\] Or\[ V=\sqrt{\frac{2mg}{\rho SC_{L}}}\] Substituting the numerical values given in the problem above (SI) gives\[ V=\sqrt{\frac{2\left ( 10675\right ) }{\rho \left ( 14.9\right ) C_{L}}}\] The air density \(\rho \) is found from appendix D since the airplane is at see level. From appendix D \begin{align*} \rho & =2.3769\times 10^{-3}\text{ lb }\frac{\text{sec}^{2}}{ft^{4}}\\ & = \boxed{1.225\text{ kg/m}^{3} } \end{align*}
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} The meaning of these coefficients is explained more in this table
| \(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 \] Or \[ \left ( -0.2797\times -0.2737-0.304\times 0.1946\right ) \cos \theta _{0}+ \left ( -0.3295\times 0.1946-\left ( -0.2797\right ) \times -0.04073\right ) \sin \theta _{0}>0 \] Hence\begin{align*} 1.739\,5\times 10^{-2}\cos \theta _{0}-7.551\,3\times 10^{-2}\sin \theta _{0} & >0\\ 1.739\,5\times 10^{-2}-7.551\,3\times 10^{-2}\tan \theta _{0} & >0\\ \tan \theta _{0} & >\frac{1.739\,5\times 10^{-2}}{7.551\,3\times 10^{-2}}\\ \tan \theta _{0} & >0.230\,36 \end{align*}
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} The appendix is now used to replace the expressions in the above. From appendix B.9 and for \(A>1\) \begin{equation} C_{l_{\beta }}=C_{L}\left [ \left ( \frac{C_{l_{\beta }}}{C_{L}}\right ) _{\Lambda _{c/2}}K_{M_{A}}+\left ( \frac{C_{l_{\beta }}}{C_{L}}\right ) _{A}\right ] +\Gamma \left ( \frac{C_{l_{\beta }}}{\Gamma }K_{M_{\Gamma }}\right ) \tag{3} \end{equation} The rest of the expression in B.9,1 was not used, since \(\theta =0\). For \(C_{n_{\beta }}\), looking at equation B.9.3 in the appendix B.9, shows it does not depend on \(\Gamma \), therefore its current numerical value from the above table is used \begin{equation} C_{n_{\beta }}=0.1946 \tag{4} \end{equation} Looking at \(C_{l_{r}}\), and from appendix B.11 \begin{equation} C_{l_{r}}=C_{L}\left ( \frac{C_{L_{r}}}{C_{L}}\right ) _{C_{L_{M}}=0}+\left ( \frac{\Delta C_{l_{r}}}{\Gamma }\right ) \Gamma \tag{5} \end{equation} Where the last term was not used since \(\theta =0\). Finally from B.11.4 one sees that \(C_{n_{r}}\) does not depend on \(\Gamma \), therefore its current numerical value is used \begin{equation} C_{n_{r}}=-0.2737 \tag{6} \end{equation} Substituting (3,4,5,6) into (2) results in \begin{multline} -0.2737\overbrace{\left ( C_{L}\left [ \left ( \frac{C_{l_{\beta }}}{C_{L}}\right ) _{\Lambda _{c/2}}K_{M_{A}}+\left ( \frac{C_{l_{\beta }}}{C_{L}}\right ) _{A}\right ] +\Gamma \left ( \frac{C_{l_{\beta }}}{\Gamma }K_{M_{\Gamma }}\right ) \right ) }^{C_{l_{\beta }}} \tag{7}\\ -0.1946\overbrace{\left ( C_{L}\left ( \frac{C_{L_{r}}}{C_{L}}\right ) _{C_{L_{M}}=0}+\left ( \frac{\Delta C_{l_{r}}}{\Gamma }\right ) \Gamma \right ) }^{C_{l_{r}}} >0\nonumber \end{multline} In the above, \(\frac{\Delta C_{l_{r}}}{\Gamma }\) is found from (B.11,3) as\[ \frac{\Delta C_{l_{r}}}{\Gamma }=\frac{1}{12}\frac{\pi A\sin \Lambda _{c/4}}{A+4\cos \Lambda _{c/4}}\] Hence (7) becomes\begin{multline} -0.2737\overbrace{\left ( C_{L}\left [ \left ( \frac{C_{l_{\beta }}}{C_{L}}\right ) _{\Lambda _{c/2}}K_{M_{A}}+\left ( \frac{C_{l_{\beta }}}{C_{L}}\right ) _{A}\right ] +\Gamma \left ( \frac{C_{l_{\beta }}}{\Gamma }K_{M_{\Gamma }}\right ) \right ) }^{C_{l_{\beta }}} \tag{8}\\ -0.1946\overbrace{\left ( C_{L}\left ( \frac{C_{L_{r}}}{C_{L}}\right ) _{C_{L_{M}}=0}+\left ( \frac{1}{12}\frac{\pi A\sin \Lambda _{c/4}}{A+4\cos \Lambda _{c/4}}\right ) \Gamma \right ) }^{C_{l_{r}}} >0\nonumber \end{multline} In order to determine what happens as \(\Gamma \) changes, it is assumed that all terms that do not involve \(\Gamma \) above are fixed for the time being and can be renamed to constants \(z_{1}\) and \(z_{2}\) for the purpose of finding the effect of changing \(\Gamma \). Writing (8) as\begin{equation} -z_{1}\Gamma \left ( \frac{C_{l_{\beta }}}{\Gamma }K_{M_{\Gamma }}\right ) -z_{2}\frac{1}{12}\frac{\pi A\sin \Lambda _{c/4}}{A+4\cos \Lambda _{c/4}}\Gamma >0 \tag{9} \end{equation} This above can be simplified more since terms such as \(\pi \) and \(A\) are fixed and do not change with changing \(\Gamma \). In addition, \(\Lambda _{c/4}\) do not change with \(\Gamma \). The above is simplified more resulting in \begin{equation} -\overbrace{z_{1}\Gamma \left ( \frac{C_{l_{\beta }}}{\Gamma }K_{M_{\Gamma }}\right ) }^{C_{l_{\beta }}} -\overbrace{z_{2}\Gamma }^{C_{l_{r}}}>0 \tag{10} \end{equation} The LHS above is required to remain positive for stability. The more positive it is, the more stable the system is. As \(\Gamma \) increases, then \(C_{l_{r}}=z_{2}\Gamma \) will become more positive. But what happens to \(\Gamma \left ( \frac{C_{l_{\beta }}}{\Gamma }K_{M_{\Gamma }}\right ) \) as \(\Gamma \) increases? To find what happens to \(\frac{C_{l_{\beta }}}{\Gamma }\) and what happens to \(K_{M_{\Gamma }}\), the hints given are used. From figure B.9,4 it is seen that \(\frac{C_{l_{\beta }}}{\Gamma }\) is negative curve, as the y-axis is negative for positive aspect ratio (the x-axis). This was the case for all values of the taper ratios \(\lambda \).
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} or\begin{equation} \overbrace{\Gamma \left ( \operatorname{positive}\right ) }^{C_{l_{\beta }}} +\overbrace{\Gamma \left ( \operatorname{negative}\right ) }^{C_{l_{r}}} >0 \tag{12} \end{equation} \(E\) has to be positive for spiral stability. Between the two terms above, \(C_{l_{\beta }}\) will cause \(E\) to be become more positive as \(\Gamma \) increases. While the second term \(C_{l_{r}}\) will have the opposite effect, it will reduce \(E\) and cause it to become negative.
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} Where \begin{align*} I_{x}^{\prime } & =\frac{\left ( I_{x}I_{z}-I_{zx}^{2}\right ) }{I_{z}}=\frac{\left ( 230\times 1778-0\right ) }{1778}=230\text{ kg m}^{2}\\ I_{z}^{\prime } & =\frac{\left ( I_{x}I_{z}-I_{zx}^{2}\right ) }{I_{x}}=\frac{\left ( 230\times 1778-I_{zx}^{2}\right ) }{230}=1778\text{ kg m}^{2}\\ I_{zx}^{\prime } & =\frac{I_{zx}}{\left ( I_{x}I_{z}-I_{xz}^{2}\right ) }=0 \end{align*}
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}\] But at trim, \(L=mg\), hence\[ C_{L}=\frac{mg}{\frac{1}{2}\rho V^{2}S}=\frac{10675}{\frac{1}{2}\left ( 1.225\right ) \left ( 77.3\right ) ^{2}\left ( 14.9\right ) }= \boxed{0.196} \] Now the numerical values of the dimensional derivatives are calculated
| 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} \] And since \(I_{zx}^{\prime }=0\) the above reduces to\[\begin{pmatrix} \dot{v}\\ \dot{p}\\ \dot{r}\\ \dot{\phi }\end{pmatrix} =\begin{pmatrix} -0.09076\,1 & -0.115\,54 & -76.811 & 9.81\\ -2.434\,8 & -56.5 & 4.951\,7 & 0\\ 0.05045\, & -0.331\,47 & -0.810\,63 & 0\\ 0 & 1 & 0 & 0 \end{pmatrix}\begin{pmatrix} v\\ p\\ r\\ \phi \end{pmatrix} \] The above is in the form of \(\dot{x}=Ax\) . The eigenvalues of \(A\) are found using Matlab. They are
\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}\] Where in the above, the vector \(\begin{pmatrix} x_{1i}\\ x_{2i}\\ x_{3i}\\ x_{4i}\end{pmatrix} \) is the \(i^{th}\) eigenvector that corresponds to the \(i^{th}\) eigenvalue. The above can also be written as\[ \boxed{ \mathbf{x}=\mathbf{x}_{1}e^{\lambda _{1}t}+\mathbf{x}_{2}e^{\lambda _{2}t}+\mathbf{x}_{3}e^{\lambda _{3}t}+\mathbf{x}_{4}e^{\lambda _{4}t} } \] Where the eigenvalues \(\lambda _{i}\) are known, but not the eigenvectors \(\mathbf{x}_{i}\). By definition, an eigenvector \(x_{i}\) corresponding to eigenvalue \(\lambda _{i}\) can be found from \begin{align*} Ax_{i} & =\lambda _{i}x_{i}\,\\ \left ( A-\lambda _{i}I\right ) x_{i} & =0 \end{align*}
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 above is in the form of \(\dot{x}=Ax\) . The eigenvalues of \(A\) (using eig command) are \begin{align*} \lambda _{dutch} & =-0.44044\pm 1.9015i\\ \lambda _{rolling} & =-56.505 \end{align*}
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}}} \] This has 4 roots. They all should have negative real part for stability. Using the hint given, we set this as quadratic polynomial in say \(r=s^{2}\) and solve for \(r\) first. Letting \(s^{2}=r\) in the above gives\[ p\left ( r\right ) =r^{2}+\frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}r+\frac{LM}{I_{x}I_{y}}\] This can be solved using the quadratic equation\begin{align*} r & =s^{2}=\frac{-b}{2a}\pm \frac{\sqrt{b^{2}-4ac}}{2a}\\ s^{2} & =-\frac{1}{2}\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) \pm \frac{1}{2}\sqrt{\left ( \frac{H^{2}-LI_{y}-MI_{x}}{I_{x}I_{y}}\right ) ^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) } \end{align*}
\(\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 ) }\] There are two sub-cases to consider under this case. If the expression \(A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) \) under the square root is positive or negative. Let \(A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) =B\) therefore
If \(B<0\) then\[ s^{2}=-\frac{1}{2}A+\frac{1}{2}i\sqrt{\left \vert B\right \vert }\] This is a complex number, say \(z\). To make it easier to proceed, it is written in polar form as \(\left \vert z\right \vert e^{i\arg \left ( z\right ) }\). Hence in polar form\[ s^{2}=z=\left \vert z\right \vert e^{i\arg z}\] Where \(|z|\) is the magnitude of the complex root, which is always positive. Taking the square root of the complex root gives \begin{align*} s_{1,2} & =\pm \left ( \left \vert z\right \vert e^{i\arg z}\right ) ^{\frac{1}{2}}\\ & =\pm \sqrt{\left \vert z\right \vert }e^{\frac{i}{2}\arg z}\\ & =\{+\sqrt{\left \vert z\right \vert }e^{\frac{i}{2}\arg z},-\sqrt{z}e^{\frac{i}{2}\arg z}\} \end{align*}
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 ) \] Now we can take the square root\begin{align*} s & =\pm \sqrt{\frac{1}{2}}\sqrt{\left \vert D\right \vert -A}\\ & =\left \{ \sqrt{\frac{1}{2}}\sqrt{\left \vert D\right \vert -A},-\sqrt{\frac{1}{2}}\sqrt{\left \vert D\right \vert -A}\right \} \end{align*}
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 ) }\] There is now 2 sub cases to consider. If the expression \(\sqrt{A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) }\)under the square root is positive or negative. Let \(A^{2}-4\left ( \frac{LM}{I_{x}I_{y}}\right ) =B\) therefore
If \(B<0\) then\[ s^{2}=-\frac{1}{2}A-\frac{1}{2}i\sqrt{\left \vert B\right \vert }\] This is a complex number. To make it easier to proceed, it is written in polar form as \(\left \vert s^{2}\right \vert e^{i\arg \left ( s^{2}\right ) }\). Hence in polar form\[ s^{2}=\left \vert s^{2}\right \vert e^{i\arg \left ( s^{2}\right ) }\] This is the same as case 1. Will just use that result. hence this is unstable condition\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\), 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.