Modeling and control of hypersonic vehicle dynamic under centroid shift

Abstract

Due to the huge flight scale and the fast speed of hypersonic vehicle, the system must be of strong nonlinearity, coupling, and fast time variability, which give rise to the huge challenge for the design of controller. The good control performance must be based on the elaborately designed controller, which is established in the carefully designed center of mass. Once the center of mass moves unexpectedly, it is bound to affect the ability of existing flight controller to maintain the stability of hypersonic vehicle, resulting in serious consequences, even loss of control. Based on Newton’s laws and Varignon’s theory, a mathematical model for hypersonic vehicle with centroid shift is built up to research the influence of centroid on the motion of hypersonic vehicle. The zero-input response tests are conducted from the different aircraft body axes of the coordinate. Simulation results show that such influence is coupling, abrupt, irregular, and time-variant. In order to inhibit the bad influence of unexpected centroid shift, terminal sliding mode controller combined with radial basis function neural networks and just terminal sliding mode controller alone are adopted to handle such problems in view of robust control itself and auxiliary compensation. Simulation results show that such influence can be inhibited and compensated in a certain region, and the further research is still needed.

Introduction

Hypersonic vehicle (HSV) is an aerospace vehicle flying at Mach 5 or above. ¹ Although HSV possesses many advantages, such as a large flight envelope, high maneuverability, and good penetrability, it is very sensitive to the variation of flight conditions owing to the ultra-high altitude and Mach numbers. Besides, its dynamic model also reflects fast time varying, high nonlinear, strong coupling, and large parametric uncertainties. ² Therefore, HSV is likely to suffer a variety of uncertainties, disturbances, and faults that come from the internal and external, such as un-modeled dynamics, external disturbances, fuel sloshing, some damages, or ablation to fuselage, and some faults of control actuators or sensors, which undoubtedly pose a great challenge to the control design of HSV.

As for the anti-disturbance, anti-uncertainties, and fault-tolerant control, a large number of controllers have been designed to handle those problems, involving the state feedback control, fuzzy control, sliding mode control, adaptive control, neural networks, based on data driven, and so on, seen in previous studies,^3–7 achieving the good control effects. However, the ideal control effect is based on the carefully designed center of gravity and the high reliable controller. Center of gravity is crucial for the controller design. Once the centroid shifts unexpectedly, it must arise HSV out of the trim that can adversely affect the ability of the existing flight control system to maintain the stability of HSV, even causing the out of control or instability. As usual, the centroid shift of airplane is considered as the harmful factor, coming from consumption of fuel, airdropping freight, fuselage damage, and others, ⁸ resulting in the disastrous consequences for the control of the aircraft. However, not all centroid shifts will cause the loss of control for the aircraft. The World War II aviation history was filled with the stories of aircrafts coming back home safely despite suffering serious centroid shift. ⁹ In recent years, the mass moving control (controlled centroid shift) has been widely spread.^10,11 For one reason, the actuator of mass moving control is located in the inside of the aircraft, thus avoiding the problem of the ablation of control surface existing in aerodynamic control scheme and maintaining the good aerodynamic shape of the aircraft. For another reason, the mass moving control can make up the lack of aerodynamic force under the rarefied air and improve the maneuvering performance of the aircraft, which has become one of the hot topics in the field of HSV control.

Flight control of HSV with centroid shift in off-nominal flight conditions poses enormous challenges in many fields, including aerodynamics, structural dynamics, flight dynamics, control, and so on. Thus, a comprehensive investigation from the aircraft-integrated system perspective is needed to research, aiming to provide an integrated approach to centroid shift effect physics-based modeling and simulation, safety-of-flight assessment, flight control, and recovery. ⁸ The general equations of motion for damaged asymmetric aircraft are established in the work by Bacon and Irene, ¹³ based on Newton’s laws and Varignon’s theory, which describes a set of flight dynamics equations of motion for a rigid body not necessarily referenced to the body’s center of mass. Then, the aircraft described by new model is desired to track the motion of the body’s previous center of mass reference frame when the body loses a portion of its mass. The stability recovery of damaged asymmetric aircraft is conducted by Nguyen et al. ⁸ Based on a neural network parameter estimation blended with a direct adaptive law, a hybrid adaptive control method is proposed for the control of aircraft with uncertain structural change. The stability and convergence of this control strategy are presented based on Lyapunov theory. However, the model of aircraft, in the works by Bacon and Irene ¹³ and Nguyen et al., ⁸ with centroid shift caused by damage is established in the aircraft body coordinate frame, which makes it hard to analyze the effect of aerodynamic force and torque on the motion of HSV due to centroid shift. The aircraft models with centroid shift of Wang and colleagues^10,11 are built up for designing the controller instead of aircraft itself, which has fewer descriptions for the aerodynamic parameters. Based on the above research, it is a difficult task that how to obtain the proper transform matrix to establish the dynamic model of HSV with centroid shift in view of wind coordinate frame. Besides, since the center of mass is hard to measure, the variations of aerodynamic parameters and torques due to the centroid shift are also hard to obtain. In this article, the detailed deductions are provided including transformation matrix of coordinates, the process of modeling deduction, and the variation of aerodynamic force and torque.

This article focuses on modeling and control of HSV with centroid shift and the sections are organized as follows. Section “Modeling of HSV with centroid shift” provides the detailed process of modeling HSV with centroid shift in view of wind coordinate frame. According to Newton’s laws of motion and Varignon’s theory, the dynamic equations of HSV with centroid shift are established form two parts: linear acceleration and angular acceleration. Through a large number of coordinate transformations, the force and torque are transformed into the proper coordinate frame with the huge calculations. Then, the dynamic equations of HSV in view of wind coordinate frame are obtained. In section “Analysis of the aerodynamic, propulsive forces, and torques of HSV with centroid shift,” based on the parallel axis theorem, perpendicular axis theorem, and stretch rule, the influence of centroid shift on the aerodynamic force and torque and the inertia matrix are established. In section “Analysis of the effect of centroid shift on the motion of HSV,” taking the angular acceleration of HSV as an example, based on MATLAB, a zero-input response test is applied to analyze the characteristic of the motion of HSV with centroid shift along with single axis. In section “Example: dynamic response of HSV with centroid shift,” terminal sliding mode controller (TSMC) combined with radial basis function (RBF) neural networks observer and TSMC itself alone are applied to handle such problems in view of robust control itself and auxiliary compensation. The corresponding simulation results show that the effect of centroid shift on HSV can be handled in a certain range but more research works are also needed.

Modeling of HSV with centroid shift

Geometric characteristics of HSV

The vertical view of the HSV is shown in Figure 1, while the details of geometric characteristics and parameters can be seen in the work. ¹²

Figure 1.

Vertical view of the aerodynamic configuration of HSV.

The three-dimensional and two-dimensional views of centroid shift are shown in Figure 2. For establishing the dynamic model of HSV with centroid shift, Newton’s law is adopted in this article. ¹⁴ Then, the dynamic equations based on the flat non-rotating earth frame, shown in Figure 2, are developed, which is applied to describe the conditions that the body loses parts of its mass, fuel consumption, or other reasons giving rise to the centroid shifting to new place $O{\prime }_b$ . The variation of $O_b\to O{\prime }_b$ is shown in Figure 2(a). Z is an arbitrary point of the aircraft body in Figure 2(b), while $O{\prime }_b$ is a new place of the center of the mass in Figure 2(a).

Figure 2.

Arbitrary body referenced in inertial and body-centric frames: (a) three-dimensional view of centroid shift, $O_b\to O{\prime }_b$ and (b) two-dimensional view of centroid shift. ¹³

Linear acceleration

Considering the parts of HSV in Figure 2(b), the $x_g-y_g-z_g$ frame is selected as the inertial reference at original point $O_b$ , and $x_b-y_b-z_b$ frame is fixed on the body of HSV at any point Z. The body of HSV can move freely and possesses the relative angular velocity $\omega$ to $x_g-y_g-z_g$ frame. Order $m_i$ represents a part of mass of HSV defined in the $x_g-y_g-z_g$ frame as

$${\mathcal{R}}_i={\mathcal{R}}_Z+{\rho }_i$$

(1)

where ${\mathcal{R}}_Z$ represents the distance of Z and O. ${\rho }_i$ is the distance between $m_i$ and Z. Thus, the centroid G is described by $\overline{\rho }$ with respect to Z and $\overline{\mathcal{R}}$ with respect to O. According to Varignon’s theory, $\overline{\rho }$ and $\overline{\mathcal{R}}$ can be formulated as $m\overline{\rho }=\sum m_i{\rho }_i$ and $m\overline{\mathcal{R}}=\sum m_i{\mathcal{R}}_i$ , respectively, where m denotes the whole mass of HSV, satisfying $m=\sum m_i$ . Based on Newton’s second law, the linear motion of HSV is charged by

$$\sum F=\sum m_i({\stackrel{··}{\overline{\mathcal{R}}}}_i{)}_{x_g{y}_g{z}_g}=m(\stackrel{··}{\overline{\mathcal{R}}}{)}_{x_g{y}_g{z}_g}$$

(2)

where $\sum F$ denotes the entire external forces applied on HSV. $(\stackrel{··}{\overline{\mathcal{R}}}{)}_{x_g{y}_g{z}_g}$ represents the acceleration of centroid with respect to $x_g-y_g-z_g$ frame.

In order to obtain $(\stackrel{··}{\overline{\mathcal{R}}}{)}_{x_g{y}_g{z}_g}$ , take two time derivations of equation (1) under the $x_g-y_g-z_g$ frame, complying with the rules ¹⁵

$${\left(\frac{d[·]}{\mathrm{dt}}\right)}_{x_g{y}_g{z}_g}={\left(\frac{d[·]}{\mathrm{dt}}\right)}_{x_b{y}_b{z}_b}+\omega \times [·]$$

(3)

So $({\stackrel{·}{\mathcal{R}}}_i{)}_{x_g{y}_g{z}_g}=({\stackrel{·}{\mathcal{R}}}_Z{)}_{x_g{y}_g{z}_g}+({\stackrel{·}{\rho }}_i{)}_{x_g{y}_g{z}_g}$ becomes

$$({\stackrel{·}{\mathcal{R}}}_i{)}_{x_g{y}_g{z}_g}=v_Z+{\stackrel{·}{\rho }}_i+\omega \times {\rho }_i$$

(4)

where $v_Z$ is $({\stackrel{·}{\mathcal{R}}}_Z{)}_{x_g{y}_g{z}_g}$ described in $x_b-y_b-z_b$ frame. ${\rho }_i$ , $\omega$ , and ${\stackrel{·}{\rho }}_i$ are all defined in $x_b-y_b-z_b$ frame. In order to obtain the concise explanation, $(·{)}_{x_b{y}_b{z}_b}$ represents the derivation under the $x_b-y_b-z_b$ frame. Under the $x_g-y_g-z_g$ frame, making the derivation of equation (4) yields

$${({\stackrel{··}{\mathcal{R}}}_i)}_{x_g{y}_g{z}_g}={\left(\frac{d(v_{{\rm Z}}+{\stackrel{.}{\rho }}_i+\omega \times {\rho }_i)}{dt}\right)}_{x_b{y}_b{z}_b}+\omega \times (v_{{\rm Z}}+{\stackrel{.}{\rho }}_i+\omega \times {\rho }_i)$$

(5)

Then, we have

$${({\stackrel{··}{\mathcal{R}}}_i)}_{x_g{y}_g{z}_g}={\stackrel{.}{v}}_{{\rm Z}}+\omega \times v_{{\rm Z}}+{\stackrel{··}{\rho }}_i+\stackrel{.}{\omega }\times {\rho }_i+2(\omega \times {\stackrel{.}{\rho }}_i)+\omega \times (\omega \times {\rho }_i)$$

(6)

By substituting equation (6) into equation (2) under the $x_g-y_g-z_g$ frame, we have

$${\displaystyle \sum F}={\displaystyle \sum m_i}({\stackrel{.}{v}}_Z+\omega \times v_Z+{\stackrel{··}{\rho }}_i+\stackrel{.}{\omega }\times {\rho }_i+2(\omega \times {\stackrel{.}{\rho }}_i)+\omega \times (\omega \times {\rho }_i))$$

(7)

$${\displaystyle \sum F}=m({\stackrel{.}{v}}_Z+\omega \times v_Z)+{\displaystyle \sum m_i}{\stackrel{··}{\rho }}_i+2(\omega \times {\displaystyle \sum m_i}{\stackrel{.}{\rho }}_i)+\stackrel{.}{\omega }\times {\displaystyle \sum m_i}{\rho }_i+\omega \times (\omega \times {\displaystyle \sum m_i}{\rho }_i)$$

(8)

$${\displaystyle \sum F}=m({\stackrel{.}{v}}_Z+\omega \times v_Z)+m\stackrel{··}{\overline{\rho }}+2(\omega \times m\stackrel{.}{\overline{\rho }})+\stackrel{.}{\omega }\times m\overline{\rho }+\omega \times (\omega \times m\overline{\rho })$$

(9)

Under the assumption that the body of HSV is rigid, we have that $\stackrel{··}{\overline{\rho }}=\stackrel{·}{\overline{\rho }}=0$ . Equation (9) is rewritten as follows

$$\sum F=m({\stackrel{·}{v}}_Z+\omega \times v_Z)+\stackrel{·}{\omega }\times m\overline{\rho }+\omega \times (\omega \times m\overline{\rho })$$

(10)

The angular velocity, location of centroid, and the velocity at point Z are defined in $x_b-y_b-z_b$ frame to formulate the dynamic equations of HSV and the details are shown as follows

$$\omega =pi+qj+rk,\overline{\rho }=\Delta xi+\Delta yj+\Delta zk,v_Z=U_{Z}i+V_{Z}j+W_{Z}k$$

(11)

Besides, the gravity is mapped to the aircraft body frame and we have

$$W=-\mathrm{mg}\sin \theta i+\mathrm{mg}\cos \theta \sin \varphi j+\mathrm{mg}\cos \theta \cos \varphi k$$

(12)

Ultimately, the force balance equations with respect to an arbitrary point Z of HSV are obtained as follows ¹⁶

$${\stackrel{.}{U}}_Z=\frac{{\displaystyle \sum F_X}}{m}+r{V}_Z-q{W}_Z-g\sin \theta +(q^2+r^2)\Delta x-(qp-\stackrel{.}{r})\Delta y-(pr+\stackrel{.}{q})\Delta z$$

(13)

$${\stackrel{.}{V}}_Z=\frac{{\displaystyle \sum F_Y}}{m}-r{U}_Z+p{W}_Z+g\cos \theta \sin \varphi -(pq+\stackrel{.}{r})\Delta x+(p^2+r^2)\Delta y-(qr-\stackrel{.}{p})\Delta z$$

(14)

$${\stackrel{.}{W}}_Z=\frac{{\displaystyle \sum F_Z}}{m}+q{U}_Z-p{V}_Z+g\cos \theta \cos \varphi -(pr-\stackrel{.}{q})\Delta x-(qr+\stackrel{.}{p})\Delta y+(p^2+q^2)\Delta z$$

(15)

where $\varphi$ is the roll Euler angle defined in aircraft body frame with respect to the fixed inertial reference frame $(x_g-y_g-z_g)$ . $\theta$ denotes the pitch Euler angle under the same coordinate frame with $\varphi$ . The items containing $[\Delta x,\Delta y,\Delta z]$ in equations (13)–(15) represent the distance by centroid deviating from the aircraft body fixed frame, $x_b-y_b-z_b$ .

In order to obtain the rotation matrix among the four coordinate frames, the relationships of four coordinate frames and eight Euler angles need to be defined clearly, and the mapping relationships between the reference frames are described in Figure 3. For instance, the relationship of $O_b$ and $O_w$ can be described as $S_{\alpha ,\beta }^T$ , given in equation (90).

Figure 3.

The mapping relationship between different reference frames.

Force equations of HSV with centroid shift in wind coordinate frame

According to the requirements, the states describing the motion of HSV can be selected between Euler angles $(x,y,z,u_b,v_b,w_b,\varphi ,\theta ,\psi ,p,q,r)$ and aerodynamic angles $(x,y,z,V,\chi ,\gamma ,\alpha ,\beta ,\mu ,p,q,r)$ . Euler angles are built based on the aircraft body frame, while aerodynamic angles are established in wind coordinate frame. Since scramjet is rather sensitive to the changes of aerodynamic angles, angle of attack, and sideslip angle, some large changes of $\alpha$ and $\beta$ may give rise to the engine deviating from its designed state, making the engine flame out. Therefore, in order to control the aerodynamic angles accurately guaranteeing the safety of HSV, the second states are more proper and more suitable for the following study.

The velocity and its derivation in aircraft body frame need to be transformed into the wind frame through the matrix $S_{\alpha ,\beta }^T$ , as shown in Figure 3. The details are provided as follows

$$\left[\begin{array}{c}{U}_Z\hfill \\ \multicolumn{1}{c}{V_Z}\\ \multicolumn{1}{c}{W_Z}\end{array}\right]=S_{\alpha ,\beta }^T\left[\begin{array}{c}V\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right]=\left[\begin{array}{c}V\cos \alpha \cos \beta \hfill \\ \multicolumn{1}{c}{V\sin \beta }\\ \multicolumn{1}{c}{V\sin \alpha \cos \beta }\end{array}\right]$$

(16)

where $S_{\alpha ,\beta }^T$ is given in equation (90).

Then, take the derivative of equation (16) with respect to t, and one has

$$\left[\begin{array}{c}{\stackrel{·}{U}}_Z\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{V}}_Z}\\ \multicolumn{1}{c}{{\stackrel{·}{W}}_Z}\end{array}\right]=S_{\alpha ,\beta }^T\left[\begin{array}{c}\stackrel{·}{V}\hfill \\ \multicolumn{1}{c}{V\stackrel{·}{\beta }}\\ \multicolumn{1}{c}{V\cos \beta \stackrel{·}{\alpha }}\end{array}\right]$$

(17)

where $S_{\alpha ,\beta }^T$ is given in equation (90).

As is well known, the aircraft is charged by the force of aerodynamic $(f_a^g)$ , thrust $(f_T^g)$ , and gravity $(g^g)$ . The thrust vector $(f_T^g)$ is usually expressed in the aircraft body frame. The aerodynamic force vector $(f_a^g)$ is defined in wind frame. The acceleration of gravity $(g^g)$ is shown in earth-surface inertial reference frame. The specifics are as follows

$$f_T^g={[T,0,0]}^T,f_a^g=[-D,Y,-L{]}^T,g^g={[0,0,g]}^T$$

(18)

To obtain the dynamic equations about the aerodynamic angles $(\alpha ,\beta )$ , those forces need to be converted to the wind frame by the rotation matrix $S_{\alpha ,\beta }$ . The specifics are as follows

$${\displaystyle \sum F_{wind}}=S_{\alpha ,\beta }{\left[T_x,T_y,T_z\right]}_{body}^T+{\left[-D,Y,-L\right]}_{wind}^T+S_{ag}{\left[0,0,mg\right]}_{ground}^T$$

(19)

where the concrete expansion of equation (19) can be seen in equations (91)–(93) and $S_{\mathrm{ag}}$ is seen in equation equations (94).

Then, substitute equations (16)–(19) into equations (13)–(15), and the dynamic equations of HSV in wind coordinate frame are obtained as follows

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{V}=\frac{(-D-\mathrm{mg}\sin \gamma )}{m}}\\ +\frac{(T_x\cos \beta \cos \alpha +T_y\sin \beta +T_z\cos \beta \sin \alpha )}{m}\hfill \\ +\frac{(c_{11}\Delta x+c_{12}\Delta y+c_{13}\Delta z)}{m}\hfill \end{array}$$

(20)

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\beta }=p\sin \alpha -r\cos \alpha +\frac{(Y+\mathrm{mg}\cos \gamma \sin \mu -T_x\sin \beta \cos \alpha +T_y\cos \beta -T_z\sin \beta \sin \alpha )}{\mathrm{mV}}}\\ +\frac{(c_{21}\Delta x+c_{22}\Delta y+c_{23}\Delta z)}{\mathrm{mV}}\hfill \end{array}$$

(21)

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\alpha }=q-(p\cos \alpha +r\sin \alpha )\tan \beta +\frac{(-L+\mathrm{mg}\cos \gamma \cos \mu -T_x\sin \alpha +T_z\cos \alpha )}{(\mathrm{mV}\cos \beta )}}\\ +\frac{(c_{31}\Delta x+c_{32}\Delta y+c_{33}\Delta z)}{(\mathrm{mV}\cos \beta )}\hfill \end{array}$$

(22)

where $c_{i,j},i,j=1,2,3$ can be seen in equation (95).

Force equations of HSV without centroid shift in flight-path coordinate frame

In the earth-surface inertial reference frame, according to Newton’s second law again, we obtain

$$m\frac{d{V}^g}{\mathrm{dt}}=m{g}^g+f_a^g+f_T^g$$

(23)

where $(·{)}^g$ represents the states of HSV described in the earth-surface inertial reference frame.

In order to distinguish the representation of V in different coordinate systems, V is remarked as $V_p$ in this section, which means that it is located in flight-path coordinate frame

$$V^p={[V_p,0,0]}^T$$

(24)

In order to obtain the dynamic equations with respect to flight-path angles, using the formula (3), we have

$$\frac{d{V}^g}{\mathrm{dt}}=S_p^g\frac{d{V}^p}{\mathrm{dt}}+{\Omega }_p^g\times V^g$$

(25)

where $S_p^g$ is a transformation matrix to map the states of HSV in flight-path frame into the earth-surface inertial frame. The transformation matrix between flight-path frame and earth-surface inertial frame $(S_{\gamma ,\chi })$ is given in equation (97).

Then, by substituting equation (25) into equation (23), we obtain

$$m(S_p^g\frac{d{V}^p}{\mathrm{dt}}+{\Omega }_p^g\times V^g)=S_b^g{f}_T^b+S_a^g{f}_a^a+m{g}^g$$

(26)

By multiplying $S_g^p$ on both sides of equation (26), we have

$$m(S_g^p{S}_p^g\frac{d{V}^p}{\mathrm{dt}}+S_g^p{\Omega }_p^g\times V^g)=S_g^p{S}_b^g{f}_T^b+S_g^p{S}_a^g{f}_a^a+{\mathrm{mS}}_g^p{g}^g$$

(27)

where $S_g^p$ is the unit orthogonal matrix, that is, $S_g^p{S}_p^g=I$ . Then, we have $S_g^p·({\Omega }_p^g\times V^g)={\Omega }_p^g\times (S_g^p·V^g)={\Omega }_p^g\times V^p$ . Equation (27) is rewritten as follows

$$m(\frac{d{V}^p}{\mathrm{dt}}+{\Omega }_p^g\times V^p)=S_b^p{f}_T^b+S_a^p{f}_a^a+{\mathrm{mS}}_g^p{g}^g$$

(28)

where ${\Omega }_p^g$ is the angular velocity of flight-path coordinate frame with respect to earth-surface inertial reference frame, as shown in Figure 3. The specifics are as follows ¹⁵

$${\Omega }_p^g=\left[\begin{array}{c}0\\ \frac{d\gamma }{dt}\\ 0\end{array}\right]+S_y(\gamma )\left[\begin{array}{c}0\\ 0\\ d\chi /dt\end{array}\right]=\left[\begin{array}{c}0\\ \frac{d\gamma }{dt}\\ 0\end{array}\right]+\left[\begin{array}{ccc}\cos \gamma & 0& -\sin \gamma \\ 0& 1& 0\\ \sin \gamma & 0& \cos \gamma \end{array}\right]\left[\begin{array}{c}0\\ 0\\ d\chi /dt\end{array}\right]=\left[\begin{array}{c}-\frac{d\chi }{dt}\sin \gamma \\ \frac{d\gamma }{dt}\\ \frac{d\chi }{dt}\cos \gamma \end{array}\right]$$

(29)

where $S_y(\gamma )$ means that the earth-surface inertial frame rotates $\gamma$ degrees angle referring $Y-\mathrm{axis}$ . In order to facilitate subsequent deduction, we order that ${\Omega }_p^g=[{\sigma }_1,{\sigma }_2,{\sigma }_3{]}^T$ .

Then, by substituting equation (29) into equation (28), we have

$$\left[\begin{array}{c}{\stackrel{.}{V}}_p\\ 0\\ 0\end{array}\right]+\left[\begin{array}{c}-\stackrel{.}{\chi }\sin \gamma \\ \stackrel{.}{\gamma }\\ \stackrel{.}{\chi }\cos \gamma \end{array}\right]\times \left[\begin{array}{c}{V}_p\\ 0\\ 0\end{array}\right]=\frac{1}{m}\left\{S_{\mu }^p{S}_{\alpha ,\beta }{\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]}_{body}+S_{\mu }^p{\left[\begin{array}{c}-D\\ Y\\ -L\end{array}\right]}_{wind}+S_{\gamma ,\chi }\left[\begin{array}{c}0\\ 0\\ mg\end{array}\right]\right\}$$

(30)

where $S_{\mu }^p$ , $S_{\gamma ,\chi }$ , $S_{\mu }^p{S}_{\alpha ,\beta }$ are all given in equations (96)–(98), respectively. After the simplification, the force balance equations without centroid shift in flight-path frame are obtained.

Force equations of HSV with centroid shift in flight-path coordinate frame

In order to obtain the dynamic equations of HSV with centroid shift in view of the flight-path frame and the earth-surface inertial frame, take the rigid body described in Figure 3 into account again. The centroid shift in aircraft body frame $x_b-y_b-z_b$ is mapped to flight-path frame by transformation matrix $S_{\mu }^p{S}_{\alpha ,\beta }$ , given in equation (98). The flight-path frame can move freely and has relative angular velocity $({\Omega }_p^g)$ with respect to $x_g-y_g-z_g$ frame, offered in equation (29). Referring to equation (10), we have

$$m\left(\frac{d{V}^p}{dt}+{\Omega }_p^g\times V^p+{\stackrel{.}{\Omega }}_p^g\times S_{\mu }^p{S}_{\alpha ,\beta }\overline{\rho }+{\Omega }_p^g\times ({\Omega }_p^g\times S_{\mu }^p{S}_{\alpha ,\beta }\overline{\rho })\right)=S_b^p{f}_T^b+S_a^p{f}_a^a+m{S}_g^p{g}^g$$

(31)

By substituting equations (24), (29), (30) into equation (31), we have

$$\begin{array}{c}\left[\begin{array}{c}{\stackrel{·}{V}}_p\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right]+\left[\begin{array}{c}-\stackrel{·}{\chi }\sin \gamma \hfill \\ \multicolumn{1}{c}{\stackrel{·}{\gamma }}\\ \multicolumn{1}{c}{\stackrel{·}{\chi }\cos \gamma }\end{array}\right]\times \left[\begin{array}{c}{V}_p\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right]+\left[\begin{array}{c}{\stackrel{·}{\sigma }}_1\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\sigma }}_2}\\ \multicolumn{1}{c}{{\stackrel{·}{\sigma }}_3}\end{array}\right]\hfill \\ \times S_{\mu }^p{S}_{\alpha ,\beta }\overline{\rho }+\left[\begin{array}{c}{\sigma }_1\hfill \\ \multicolumn{1}{c}{{\sigma }_2}\\ \multicolumn{1}{c}{{\sigma }_3}\end{array}\right]\times (\left[\begin{array}{c}{\sigma }_1\hfill \\ \multicolumn{1}{c}{{\sigma }_2}\\ \multicolumn{1}{c}{{\sigma }_3}\end{array}\right]\times S_{\mu }^p{S}_{\alpha ,\beta }\overline{\rho })\hfill \\ =\frac{1}{m}\{S_{\mu }^p{S}_{\alpha ,\beta }{\left[\begin{array}{c}{T}_x\hfill \\ \multicolumn{1}{c}{T_y}\\ \multicolumn{1}{c}{T_z}\end{array}\right]}_{\mathrm{body}}+S_{\mu }^p{\left[\begin{array}{c}-D\hfill \\ \multicolumn{1}{c}{Y}\\ \multicolumn{1}{c}{-L}\end{array}\right]}_{\mathrm{wind}}+S_{\gamma ,\chi }\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{\mathrm{mg}}\end{array}\right]\}\hfill \end{array}$$

(32)

where $S_{\mu }^p{S}_{\alpha ,\beta }\overline{\rho }$ is remarked as $[{\overline{\rho }}_x^p,{\overline{\rho }}_y^p,{\overline{\rho }}_z^p{]}^T$ and the specifics are described in Appendix 1 in (99).

Then, by substituting equations (29), (90), (96), (97), and (99) into equation (32), the aerodynamic equations can be obtained and more details are provided in equations (100) and (101). The specifics are provided as follows

$${\stackrel{.}{V}}_p=({\sigma }_2^2+{\sigma }_3^2){\overline{\rho }}_x^p-({\sigma }_1{\sigma }_2-{\stackrel{.}{\sigma }}_3){\overline{\rho }}_y^p-({\sigma }_1{\sigma }_3+{\stackrel{.}{\sigma }}_2){\overline{\rho }}_z^p+\frac{(-D-mg\sin \gamma )}{m}+\frac{(T_x\cos \alpha \cos \beta +T_y\sin \beta +T_z\sin \alpha \cos \beta )}{m}$$

(33)

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\chi }=\frac{(-({\sigma }_1{\sigma }_2+{\stackrel{·}{\sigma }}_3){\overline{\rho }}_x^p+({\sigma }_1^2+{\sigma }_3^2){\overline{\rho }}_y^p-({\sigma }_2{\sigma }_3-{\stackrel{·}{\sigma }}_1){\overline{\rho }}_z^p)}{(V_p\cos \gamma )}}\\ +\frac{1}{m{V}_p\cos \gamma }(T_x(-\cos \alpha \sin \beta \cos \mu +\sin \alpha \sin \mu )\hfill \\ +T_y\cos \beta \cos \mu -T_z\cos \alpha \sin \mu \hfill \\ -T_z\cos \mu \sin \alpha \sin \beta )+\frac{(Y\cos \mu +L\sin \mu )}{(m{V}_p\cos \gamma )}\hfill \end{array}$$

(34)

$$\begin{array}{c}\stackrel{·}{\gamma }=\frac{-(-({\sigma }_1{\sigma }_3-{\stackrel{·}{\sigma }}_2){\overline{\rho }}_x^p-({\sigma }_2{\sigma }_3+{\stackrel{·}{\sigma }}_1){\overline{\rho }}_y^p+({\sigma }_1^2+{\sigma }_2^2){\overline{\rho }}_z^p)}{V_p}\hfill \\ -\frac{1}{m{V}_p}(T_x(-\cos \mu \sin \alpha -\cos \alpha \sin \beta \sin \mu )\hfill \\ +T_y\cos \beta \sin \mu +T_z\cos \alpha \cos \mu \hfill \\ -T_z\sin \alpha \sin \beta \sin \mu )-\frac{(Y\sin \mu -L\cos \mu +\mathrm{mg}\cos \gamma )}{(m{V}_p)}\hfill \end{array}$$

(35)

At present, the roll angle has not been involved. According to Snell, ¹⁷ we have

$$\stackrel{·}{\mu }=\stackrel{·}{\chi }\sin \gamma -\stackrel{·}{\alpha }\sin \beta +p\cos \beta \cos \alpha +r\cos \beta \sin \alpha +q\sin \beta$$

(36)

By substituting equations (22) and (34) into equation (36), one has

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\mu }=\sec \beta (p\cos \alpha +r\sin \alpha )+\frac{\tan \gamma \{-({\sigma }_1{\sigma }_2+{\stackrel{·}{\sigma }}_3){\overline{\rho }}_x^p+({\sigma }_1^2+{\sigma }_3^2){\overline{\rho }}_y^p-({\sigma }_2{\sigma }_3-{\stackrel{·}{\sigma }}_1){\overline{\rho }}_z^p\}}{V_p}}\\ -\frac{\tan \beta }{m{V}_p}(c_{31}\Delta x+c_{32}\Delta y+c_{33}\Delta z)-\frac{\tan \beta }{m{V}_p}(-T_x\sin \alpha +T_z\cos \alpha )+\frac{\tan \gamma }{m{V}_p}(Y\cos \mu +L\sin \mu )\hfill \\ -\frac{\tan \beta }{m{V}_p}(-L+\mathrm{mg}\cos \gamma \cos \mu )+\frac{\tan \gamma }{m{V}_p}\{T_x(-\cos \alpha \sin \beta \cos \mu +\sin \alpha \sin \mu )+T_y\cos \beta \cos \mu \}\hfill \\ +\frac{\tan \gamma \{T_z(-\cos \mu \sin \alpha \sin \beta -\cos \alpha \sin \mu )\}}{m{V}_p}\hfill \end{array}$$

(37)

where $c_{i,j}(i,j=1,2,3)$ and $[{\overline{\rho }}_x^p,{\overline{\rho }}_y^p,{\overline{\rho }}_z^p{]}^T$ are given in equations (95) and (99), respectively. ${\sigma }_i(i=1,2,3)$ is given in equation (29).

Navigation equations of HSV

As usual, the trace of HSV is defined in earth-surface inertial reference frame. Referring to the work by Zong et al., ¹⁵ we have

$$\left[\begin{array}{c}\stackrel{·}{x}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{y}}\\ \multicolumn{1}{c}{\stackrel{·}{z}}\end{array}\right]=\left[\begin{array}{c}{u}_g\hfill \\ \multicolumn{1}{c}{v_g}\\ \multicolumn{1}{c}{w_g}\end{array}\right]=S_{g,b}\left[\begin{array}{c}{u}_b\hfill \\ \multicolumn{1}{c}{v_b}\\ \multicolumn{1}{c}{w_b}\end{array}\right]=S_{g,b}{S}_{b,a}\left[\begin{array}{c}V\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right]=S_{g,a}\left[\begin{array}{c}V\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right]$$

(38)

where

$$S_{g,a}=\left[\begin{array}{ccc}\cos \gamma \cos \chi \hfill & \cos \gamma \sin \chi \hfill & -\sin \gamma \hfill \\ \multicolumn{1}{c}{\sin \gamma \cos \chi \sin \mu -\sin \chi \cos \mu }& \sin \gamma \sin \chi \sin \mu +\cos \chi \cos \mu \hfill & \cos \gamma \sin \mu \hfill \\ \multicolumn{1}{c}{\sin \gamma \cos \chi \cos \mu +\sin \chi \sin \mu }& \sin \gamma \sin \chi \cos \mu -\cos \chi \sin \mu \hfill & \cos \gamma \cos \mu \hfill \end{array}\right]$$

(38)

Then, the dynamic equations of navigation are achieved as follows

$$\left[\begin{array}{c}\stackrel{·}{x}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{y}}\\ \multicolumn{1}{c}{\stackrel{·}{z}}\end{array}\right]=\left[\begin{array}{c}V\cos \gamma \cos \chi \hfill \\ \multicolumn{1}{c}{V\cos \gamma \sin \chi }\\ \multicolumn{1}{c}{-V\sin \gamma }\end{array}\right]$$

(39)

Angular acceleration

According to Figure 2(b), the $x_g-y_g-z_g$ frame is still the inertial reference with respect to $O_b$ . $x_b-y_b-z_b$ frame is attached to point Z in aircraft body. The body of HSV possesses the relative angular velocity $(\omega )$ with respect to $x_g-y_g-z_g$ frame. $m_i$ is a part of mass of HSV. In this part, the momentum principle is formulated for the moment of linear momenta taken about Z

$${\Im }_Z=\sum ({\rho }_i\times m_i{v}_i)$$

(40)

where $v_i={({\stackrel{·}{\mathcal{R}}}_i)}_{x_g{y}_g{z}_g}$ is the speed of $m_i$ with respect to $x_g-y_g-z_g$ frame. Then, taking the derivation of ${\Im }_Z$ with respect to $x_g-y_g-z_g$ frame yields

$$({\stackrel{·}{\Im }}_Z{)}_{\mathrm{XYZ}}=\sum ({({\stackrel{·}{\rho }}_i)}_{\mathrm{XYZ}}\times m_i{v}_i)+\sum ({\rho }_i\times m_i{({\stackrel{·}{v}}_i)}_{\mathrm{XYZ}})$$

(41)

where ${(·)}_{x_g{y}_g{z}_g}$ represents the states of HSV with respect to $x_g-y_g-z_g$ frame. ${({\stackrel{·}{\rho }}_i)}_{x_g{y}_g{z}_g}$ and ${({\stackrel{·}{v}}_i)}_{x_g{y}_g{z}_g}$ are the derivatives with respect to time. From Figure 2(b), we can see that ${({\stackrel{·}{\rho }}_i)}_{x_g{y}_g{z}_g}={({\stackrel{·}{\mathcal{R}}}_i-{\stackrel{·}{\mathcal{R}}}_Z)}_{x_g{y}_g{z}_g}$ and ${({\stackrel{·}{v}}_i)}_{x_g{y}_g{z}_g}={({\stackrel{··}{\mathcal{R}}}_i)}_{x_g{y}_g{z}_g}$ . Due to ${({\stackrel{·}{\mathcal{R}}}_i)}_{x_g{y}_g{z}_g}\times v_i=0$ , and the whole external torques about Z are recorded as $\sum M_Z=\sum ({\rho }_i\times m_i{({\stackrel{··}{\mathcal{R}}}_i)}_{x_g{y}_g{z}_g})$ , equation (41) is formulated as follows

$$\sum M_Z={({\stackrel{·}{\Im }}_Z)}_{x_g{y}_g{z}_g}+{({\stackrel{·}{\mathcal{R}}}_Z)}_{x_g{y}_g{z}_g}\times \sum m_i{v}_i$$

(42)

$\overline{v}$ , the speed of centroid defined in $x_g-y_g-z_g$ frame, satisfies $m\overline{v}=\sum m_i{v}_i$ . From Figure 2(b), we have that $\overline{v}=v_Z+{(\stackrel{·}{\overline{\rho }})}_{x_g{y}_g{z}_g}$ , where $v_Z(={({\mathcal{R}}_Z)}_{x_g{y}_g{z}_g})$ and ${(\stackrel{·}{\overline{\rho }})}_{x_g{y}_g{z}_g}$ are defined in $x_g-y_g-z_g$ frame. Then, momentum principle equation is obtained as follows

$$\sum M_Z=({\stackrel{·}{\Im }}_Z{)}_{x_g{y}_g{z}_g}+({\stackrel{·}{\mathcal{R}}}_Z{)}_{x_g{y}_g{z}_g}\times m(\stackrel{·}{\overline{\rho }}{)}_{x_g{y}_g{z}_g}$$

(43)

If Z is the centroid of HSV, then $\overline{v}$ equals $v_Z$ and equation (42) takes up the usual form $\sum M_Z=({\stackrel{·}{\mathcal{R}}}_Z{)}_{x_g{y}_g{z}_g}$ . For transforming equation (43) with elements defined in $x_b-y_b-z_b$ frame at Z, equation (40) is formulated as follows

$${\Im }_Z=\sum ({\rho }_i\times m_i(v_Z+{\stackrel{·}{\rho }}_i+\omega \times {\rho }_i))$$

(44)

where ${\rho }_i$ is the speed of $m_i$ defined in $x_b-y_b-z_b$ frame. ${\stackrel{·}{\rho }}_i,{\rho }_i,\omega ,\mathrm{and}{v}_Z$ are all defined in $x_b-y_b-z_b$ frame. Because HSV is rigid, we have ${\stackrel{·}{\rho }}_i=0$ . Thus

$${\Im }_Z=\sum m_i{\rho }_i\times v_Z+\sum {\rho }_i\times m_i(\omega \times {\rho }_i)$$

(45)

$${\Im }_Z=m\overline{\rho }\times v_Z+I\omega$$

(46)

where I is the torque of inertia matrix of the rigid HSV, which is defined in $x_b-y_b-z_b$ at Z and is independent of time, namely, $\stackrel{·}{I}=0$ . $I\omega$ represents the integration version of $\sum {\rho }_i\times m_i(\omega \times {\rho }_i)$ in the $x_b-y_b-z_b$ frame. Because ${\Im }_Z$ and $\overline{\rho }$ are defined in $x_b-y_b-z_b$ frame, according to equation (3), equation (43) is derived from

$$({\stackrel{·}{\Im }}_Z{)}_{x_g{y}_g{z}_g}=({\stackrel{·}{\Im }}_Z{)}_{x_b{y}_b{z}_b}+\omega \times {\Im }_Z$$

(47)

$$(\stackrel{·}{\overline{\rho }}{)}_{x_g{y}_g{z}_g}=(\stackrel{·}{\overline{\rho }}{)}_{x_b{y}_b{z}_b}+\omega \times \overline{\rho }$$

(48)

By substituting equation (46) into equation (47), we have

$$({\stackrel{·}{\Im }}_Z{)}_{x_g{y}_g{z}_g}={(\frac{d}{\mathrm{dt}}(m\overline{\rho }\times v_Z+I\omega ))}_{x_b{y}_b{z}_b}+\omega \times (m\overline{\rho }\times v_Z+I\omega )$$

(49)

Then, we obtain

$${({\stackrel{.}{\Im }}_Z)}_{x_g{y}_g{z}_g}=m\stackrel{.}{\overline{\rho }}\times v_Z+m\overline{\rho }\times {\stackrel{.}{v}}_Z+I\stackrel{.}{\omega }+\omega \times (m\overline{\rho }\times v_Z)+\omega \times I\omega$$

(50)

If the referenced axes are not noticed, those states of HSV are described in the frame of $x_b-y_b-z_b$ , namely, $(·{)}_{\mathrm{xyz}}$ . As shown in equation (48), the centroid cannot move with respect to G by the assumption that HSV is rigid and we have that $\stackrel{·}{\overline{\rho }}=0$ . Then, $(\stackrel{·}{\overline{\rho }}{)}_{\mathrm{XYZ}}=\omega \times \overline{\rho }$ is obtained. By substituting this result and equation (50) into equation (43), we have

$${\displaystyle \sum M_Z}=I\stackrel{.}{\omega }+\omega \times I\omega +m\overline{\rho }\times {\stackrel{.}{v}}_Z+m\omega \times (\overline{\rho }\times v_Z)+m{v}_Z\times (\omega \times \overline{\rho })$$

(51)

By redefining $\sum M_Z$ regardless of the torque at Z caused by the gravity force working through the centroid of HSV, we have

$${\displaystyle \sum M_Z}=I\stackrel{.}{\omega }+\omega \times I\omega +m\overline{\rho }\times {\stackrel{.}{v}}_Z+m\omega \times (\overline{\rho }\times v_Z)+m{v}_Z\times (\omega \times \overline{\rho })-\overline{\rho }\times W$$

(52)

The desired torque balance equations of HSV with respect to Z and $I\omega$ are shorthand for

$$I\omega =(I_{xx}p-I_{yy}q-I_{xz}r)i+(-I_{xy}p+I_{yy}q-I_{yz}r)j+(-I_{xz}p-I_{yz}q+I_{zz}r)k$$

(53)

without the symmetry about the $x-z$ plane, namely, $I_{\mathrm{yz}}\ne 0$ , $I_{\mathrm{xy}}\ne 0$ . Equation (52) is expanded as follows ¹⁶

$$\begin{array}{c}\multicolumn{1}{c}{L=\sum M_{Z_x}=I_{\mathrm{xx}}\stackrel{·}{p}-I_{\mathrm{xy}}\stackrel{·}{q}-I_{\mathrm{xz}}\stackrel{·}{r}+I_{\mathrm{xy}}\mathrm{pr}-I_{\mathrm{xz}}\mathrm{pq}}\\ +(I_{\mathrm{zz}}-I_{\mathrm{yy}})\mathrm{qr}+(r^2-q^2)I_{\mathrm{yz}}\hfill \\ +m((p{V}_Z-q{U}_Z+{\stackrel{·}{W}}_Z-g\cos \theta \cos \varphi )\Delta y\hfill \\ +(p{W}_Z-r{U}_Z-{\stackrel{·}{V}}_Z+g\cos \theta \sin \varphi )\Delta z)\hfill \end{array}$$

(54)

$$\begin{array}{c}\multicolumn{1}{c}{M=\sum M_{Z_y}=-I_{\mathrm{xy}}\stackrel{·}{p}+I_{\mathrm{yy}}\stackrel{·}{q}-I_{\mathrm{yz}}\stackrel{·}{r}+I_{\mathrm{yz}}\mathrm{pq}-I_{\mathrm{xy}}\mathrm{qr}}\\ +(I_{\mathrm{xx}}-I_{\mathrm{zz}})\mathrm{pr}+(p^2-r^2)I_{\mathrm{xz}}\hfill \\ +m((q{U}_Z-p{V}_Z-{\stackrel{·}{W}}_Z+g\cos \theta \cos \varphi )\Delta x\hfill \\ +(q{W}_Z-r{V}_Z+{\stackrel{·}{U}}_Z+g\sin \theta )\Delta z)\hfill \end{array}$$

(55)

$$\begin{array}{c}\multicolumn{1}{c}{N=\sum M_{Z_z}=-I_{\mathrm{xz}}\stackrel{·}{p}-I_{\mathrm{yz}}\stackrel{·}{q}+I_{\mathrm{zz}}\stackrel{·}{r}+I_{\mathrm{xz}}\mathrm{qr}-I_{\mathrm{yz}}\mathrm{pr}}\\ +(I_{\mathrm{yy}}-I_{\mathrm{xx}})\mathrm{pq}+(q^2-p^2)I_{\mathrm{xy}}\hfill \\ +m((r{U}_Z-p{W}_Z+{\stackrel{·}{V}}_Z-g\cos \theta \sin \varphi )\Delta x\hfill \\ +(r{V}_Z-q{W}_Z-{\stackrel{·}{U}}_Z-g\sin \theta )\Delta y)\hfill \end{array}$$

(56)

Then, in order to facilitate the analysis, equations (54)–(56) are reshaped as follows

$$\begin{array}{c}\left[\begin{array}{ccc}{I}_{\mathrm{xx}}\hfill & -I_{\mathrm{xy}}\hfill & -I_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{-I_{\mathrm{xy}}}& I_{\mathrm{yy}}\hfill & -I_{\mathrm{yz}}\hfill \\ \multicolumn{1}{c}{-I_{\mathrm{xz}}}& -I_{\mathrm{yz}}\hfill & I_{\mathrm{zz}}\hfill \end{array}\right]\left[\begin{array}{c}\stackrel{·}{p}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{q}}\\ \multicolumn{1}{c}{\stackrel{·}{r}}\end{array}\right]=\left[\begin{array}{c}L+\mathrm{qr}(I_{\mathrm{yy}}-I_{\mathrm{zz}})+(q^2-r^2)I_{\mathrm{yz}}+\mathrm{pq}{I}_{\mathrm{xz}}-\mathrm{pr}{I}_{\mathrm{xy}}\hfill \\ \multicolumn{1}{c}{M+\mathrm{pr}(I_{\mathrm{zz}}-I_{\mathrm{xx}})+(r^2-p^2)I_{\mathrm{xz}}+\mathrm{qr}{I}_{\mathrm{xy}}-\mathrm{pq}{I}_{\mathrm{yz}}}\\ \multicolumn{1}{c}{N+\mathrm{pq}(I_{\mathrm{xx}}-I_{\mathrm{yy}})+(p^2-q^2)I_{\mathrm{xy}}+\mathrm{pr}{I}_{\mathrm{yz}}-\mathrm{qr}{I}_{\mathrm{xz}}}\end{array}\right]\hfill \\ +\left[\begin{array}{c}m((-p{V}_Z+q{U}_Z-{\stackrel{·}{W}}_Z+g\cos \theta \cos \varphi )\Delta y+(-p{W}_Z+r{U}_Z+{\stackrel{·}{V}}_Z-g\cos \theta \sin \varphi )\Delta z)\hfill \\ \multicolumn{1}{c}{m((-q{U}_Z+p{V}_Z+{\stackrel{·}{W}}_Z-g\cos \theta \cos \varphi )\Delta x+(-q{W}_Z+r{V}_Z-{\stackrel{·}{U}}_Z-g\sin \theta )\Delta z)}\\ \multicolumn{1}{c}{m((-r{U}_Z+p{W}_Z-{\stackrel{·}{V}}_Z+g\cos \theta \sin \varphi )\Delta x+(-r{V}_Z+q{W}_Z+{\stackrel{·}{U}}_Z+g\sin \theta )\Delta y)}\end{array}\right]\hfill \end{array}$$

(57)

where I and $I^{-1}$ are given in equation (102).

Torque equations of HSV with centroid shift

Once the centroid shift appears, HSV is not a symmetrical body with respect to xoz plane; therefore, the product of inertia will not maintain 0 again, namely, $I_{\mathrm{xy}}\ne 0,I_{\mathrm{yx}}\ne 0,I_{\mathrm{yz}}\ne 0,I_{\mathrm{zy}}\ne 0$ , then equation (57) is rewritten as follows

$$\begin{array}{c}\multicolumn{1}{c}{\left[\begin{array}{c}\stackrel{·}{p}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{q}}\\ \multicolumn{1}{c}{\stackrel{·}{r}}\end{array}\right]=I^{-1}\left[\begin{array}{c}L+\mathrm{qr}(I_{\mathrm{yy}}-I_{\mathrm{zz}})+(q^2-r^2)I_{\mathrm{yz}}+\mathrm{pq}{I}_{\mathrm{xz}}-\mathrm{pr}{I}_{\mathrm{xy}}\hfill \\ \multicolumn{1}{c}{M+\mathrm{pr}(I_{\mathrm{zz}}-I_{\mathrm{xx}})+(r^2-p^2)I_{\mathrm{xz}}+\mathrm{qr}{I}_{\mathrm{xy}}-\mathrm{pq}{I}_{\mathrm{yz}}}\\ \multicolumn{1}{c}{N+\mathrm{pq}(I_{\mathrm{xx}}-I_{\mathrm{yy}})+(p^2-q^2)I_{\mathrm{xy}}+\mathrm{pr}{I}_{\mathrm{yz}}-\mathrm{qr}{I}_{\mathrm{xz}}}\end{array}\right]+m{I}^{-1}\left[\begin{array}{ccc}0\hfill & -{\stackrel{·}{W}}_Z\hfill & {\stackrel{·}{V}}_Z\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{W}}_Z}& 0\hfill & -{\stackrel{·}{U}}_Z\hfill \\ \multicolumn{1}{c}{-{\stackrel{·}{V}}_Z}& {\stackrel{·}{U}}_Z\hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}\Delta x\hfill \\ \multicolumn{1}{c}{\Delta y}\\ \multicolumn{1}{c}{\Delta z}\end{array}\right]}\\ +m{I}^{-1}\left[\begin{array}{ccc}0\hfill & -p{V}_Z+q{U}_Z\hfill & -p{W}_Z+r{U}_Z\hfill \\ \multicolumn{1}{c}{p{V}_Z-q{U}_Z}& 0\hfill & -q{W}_Z+r{V}_Z\hfill \\ \multicolumn{1}{c}{p{W}_Z-r{U}_Z}& q{W}_Z-r{V}_Z\hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}\Delta x\hfill \\ \multicolumn{1}{c}{\Delta y}\\ \multicolumn{1}{c}{\Delta z}\end{array}\right]\hfill \\ +m{I}^{-1}\left[\begin{array}{ccc}0\hfill & g\cos \gamma \cos \mu \hfill & -g\cos \gamma \sin \mu \hfill \\ \multicolumn{1}{c}{-g\cos \gamma \cos \mu }& 0\hfill & -g\sin \gamma \hfill \\ \multicolumn{1}{c}{g\cos \gamma \sin \mu }& g\sin \gamma \hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}\Delta x\hfill \\ \multicolumn{1}{c}{\Delta y}\\ \multicolumn{1}{c}{\Delta z}\end{array}\right]\hfill \end{array}$$

(58)

where the details of ${[U_Z,V_Z,W_Z]}^T$ are seen in equation (16), and $[{\stackrel{·}{U}}_Z,{\stackrel{·}{V}}_Z,{\stackrel{·}{W}}_Z{]}^T$ are described in equation (17). $I^{-1}$ is given in equation (102). So far, from equation (58), we have concluded that the centroid shift creates the additional moments on the HSV, which gives rise to the strong coupling for the motion of HSV.

Analysis of the aerodynamic, propulsive forces, and torques of HSV with centroid shift

Variation of aerodynamic parameters of HSV

Referring to the previous studies,^16,18–20 the control surfaces of HSV are composed of left and right ailerons ${\delta }_e,{\delta }_a$ and the rudder ${\delta }_r$ . The aerodynamic parameters of HSV are described as follows

$$D=\widehat{q}S(\Delta C_D+C_D^{\ast }),Y=\widehat{q}S(\Delta C_Y+C_Y^{\ast }),L=\widehat{q}S(\Delta C_L+C_L^{\ast })$$

(59)

$$l_Z=\widehat{q}bS{C}_l^{\ast }+\Delta l_Z,m_Z=m_{mrc}-X_{cg}Z+\Delta m_Z,n_Z=n_{mrc}-X_{cg}Y+\Delta n_Z$$

(60)

where

$$\begin{array}{c}\hat{q}=\frac{\rho V^2}{2},m_{\mathrm{mrc}}=\hat{q}\mathrm{cS}(\Delta C_m+C_m^{*}),n_{\mathrm{mrc}}=\hat{q}\mathrm{bS}(\Delta C_n+C_n^{*}),Z=-D\sin \alpha -L\cos \alpha \hfill \\ C_L=C_{L,\alpha }+C_{L,\delta e}{\delta }_e+C_{L,\delta a}{\delta }_a,C_D=C_{D,\alpha }+C_{D,\delta e}{\delta }_e+C_{D,\delta a}{\delta }_a+C_{D,\delta r}{\delta }_r\hfill \\ C_Y=C_{Y,\beta }\beta +C_{Y,\delta e}{\delta }_e+C_{Y,\delta a}{\delta }_a+C_{Y,\delta r}{\delta }_r\hfill \\ C_l=C_{l,\beta }\beta +C_{l,\delta e}{\delta }_e+C_{l,\delta a}{\delta }_a+C_{l,\delta r}{\delta }_r+\frac{C_{l,p}\mathrm{pb}}{2V}+\frac{C_{l,r}\mathrm{rb}}{2V}\hfill \\ C_m=C_{m,a}+C_{m,\delta e}{\delta }_e+C_{m,\delta a}{\delta }_a+C_{m,\delta r}{\delta }_r+\frac{C_{m,q}\mathrm{qc}}{2V}\hfill \\ C_n=C_{n,\beta }\beta +C_{n,\delta e}{\delta }_e+C_{n,\delta a}{\delta }_a+C_{n,\delta r}{\delta }_r+C_{n,p}\frac{p·b}{2V}+C_{n,r}\frac{r·b}{2V}\hfill \\ \Delta l_A=Y_a\Delta z-Z_a\Delta y,\Delta m_A=-X_a\Delta z+Z_a\Delta x,\Delta n_A=X_a\Delta y-Y_a\Delta x\hfill \end{array}$$

(60)

and $X_Z,Y_Z,Z_Z$ are seen in equations (91)–(93), respectively. $\Delta l_Z,\Delta m_Z,\Delta n_Z$ and $\Delta C_D,\Delta C_Y,\Delta C_L$ represent the influence of centroid shift on aerodynamic force and aerodynamic moment. In this article, the superscript * represents the force and torque coefficients for HSV without damage.

Variation of the moment of inertia of HSV

As for the moment of inertia, seen in equation (57), the inertia matrix is described with respect to point Z with the fuel consumption, fuel sloshing, or some reasons. The small part mass $m_{\mathrm{tip}}$ of HSV is lost in total at $t=k\Delta T$ , which gives rise to the centroid moving ${\rho }_{\mathrm{tip}}$ with respect to A, where ${\rho }_{\mathrm{tip}}=(\Delta x_{\mathrm{tip}},\Delta y_{\mathrm{tip}},\Delta z_{\mathrm{tip}})$ . An inertial matrix formulated about the new center of mass (tip’s), $I_{\mathrm{tip},\mathrm{cg}}$ , the remaining mass of HSV (rec’s), centroid, and inertia matrix with respect to new center of mass are as follows ¹³

$$m\to m_{\mathrm{rec}}=m-m_{\mathrm{tip}},I_{\mathrm{ij}}=I_{\mathrm{ij},\mathrm{rec}}=I_{\mathrm{ij}}-I_{\mathrm{ij},\mathrm{tip}}$$

(61)

$$I_{ii}\to I_{ii,tip}=I_{ii,tip,cm}+m_{tip}(\Delta j_{tip}^2+\Delta k_{tip}^2)i,j,k=\{x,y,z\}|i\ne j1pt\mathrm{and}i\ne k$$

(62)

$$I_{\mathrm{ij}}\to I_{\mathrm{ij},\mathrm{tip}}=I_{\mathrm{ij},\mathrm{tip},\mathrm{cm}}+m_{\mathrm{tip}}\Delta i_{\mathrm{tip}}\Delta j_{\mathrm{tip}}i,j=x,y,z|i\ne j$$

(63)

So far, based on the theorem of vertical axis and parallel axis, we can analyze out the influence of the center of mass shift on the moment of inertia.

Analysis of the effect of centroid shift on the motion of HSV

In this section, we take the angular acceleration equation (64), for example, to analyze the effect of centroid shift on the rate of attitude angular rate. In order to facilitate the analysis, equation (58) can be simplified to the affine nonlinear function and the details are given as follows

$$\stackrel{·}{\omega }=f_f+\Delta f_f+(g_f^{*}+\Delta g_f)M_c,y_f=\omega$$

(64)

where $\omega =[p,q,r{]}^T$ , $M_c$ represents the control torque, $M_c=[L_c,M_c,N_c{]}^T$ . $\Delta f_f$ is the uncertain part of $f_f$ and $\Delta g_f{M}_c$ is the one of controller, driving from the abrupt centroid shift at time T.

The specific expressions of equation (64) are given as follows

$$\left[\begin{array}{c}{f}_p\hfill \\ \multicolumn{1}{c}{f_q}\\ \multicolumn{1}{c}{f_r}\end{array}\right]=I^{-1}\left[\begin{array}{c}\mathrm{qr}(I_{\mathrm{yy}}-I_{\mathrm{zz}})+\mathrm{pq}{I}_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{\mathrm{pr}(I_{\mathrm{zz}}-I_{\mathrm{xx}})+(r^2-p^2)I_{\mathrm{xz}}}\\ \multicolumn{1}{c}{\mathrm{pq}(I_{\mathrm{xx}}-I_{\mathrm{yy}})-\mathrm{qr}{I}_{\mathrm{xz}}}\end{array}\right]$$

(65)

$$\begin{array}{cc}\multicolumn{1}{c}{\Delta f_f=}& I^{-1}\left[\begin{array}{c}(q^2-r^2)I_{\mathrm{yz}}-\mathrm{pr}{I}_{\mathrm{xy}}\hfill \\ \multicolumn{1}{c}{\mathrm{qr}{I}_{\mathrm{xy}}-\mathrm{pq}{I}_{\mathrm{yz}}}\\ \multicolumn{1}{c}{(p^2-q^2)I_{\mathrm{xy}}+\mathrm{pr}{I}_{\mathrm{yz}}}\end{array}\right]+m{I}^{-1}\left[\begin{array}{ccc}0\hfill & -{\stackrel{·}{W}}_Z\hfill & {\stackrel{·}{V}}_Z\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{W}}_Z}& 0\hfill & -{\stackrel{·}{U}}_Z\hfill \\ \multicolumn{1}{c}{-{\stackrel{·}{V}}_Z}& {\stackrel{·}{U}}_Z\hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}\Delta x\hfill \\ \multicolumn{1}{c}{\Delta y}\\ \multicolumn{1}{c}{\Delta z}\end{array}\right]\hfill \\ +m{I}^{-1}\left[\begin{array}{ccc}0\hfill & -p{V}_Z+q{U}_Z\hfill & -p{W}_Z+r{U}_Z\hfill \\ \multicolumn{1}{c}{p{V}_Z-q{U}_Z}& 0\hfill & -q{W}_Z+r{V}_Z\hfill \\ \multicolumn{1}{c}{p{W}_Z-r{U}_Z}& q{W}_Z-r{V}_Z\hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}\Delta x\hfill \\ \multicolumn{1}{c}{\Delta y}\\ \multicolumn{1}{c}{\Delta z}\end{array}\right]\hfill \\ +m{I}^{-1}\left[\begin{array}{ccc}0\hfill & g\cos \gamma \cos \mu \hfill & -g\cos \gamma \sin \mu \hfill \\ \multicolumn{1}{c}{-g\cos \gamma \cos \mu }& 0\hfill & -g\sin \gamma \hfill \\ \multicolumn{1}{c}{g\cos \gamma \sin \mu }& g\sin \gamma \hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}\Delta x\hfill \\ \multicolumn{1}{c}{\Delta y}\\ \multicolumn{1}{c}{\Delta z}\end{array}\right]\hfill \end{array}$$

(66)

$$g_f^{*}={\left[\begin{array}{ccc}{I}_{\mathrm{xx}}\hfill & 0\hfill & -I_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{0}& I_{\mathrm{yy}}\hfill & 0\hfill \\ \multicolumn{1}{c}{-I_{\mathrm{xz}}}& 0\hfill & I_{\mathrm{zz}}\hfill \end{array}\right]}^{-1},\Delta g_f={\left[\begin{array}{ccc}0\hfill & -I_{\mathrm{xy}}\hfill & 0\hfill \\ \multicolumn{1}{c}{-I_{\mathrm{xy}}}& 0\hfill & -I_{\mathrm{yz}}\hfill \\ \multicolumn{1}{c}{0}& -I_{\mathrm{yz}}\hfill & 0\hfill \end{array}\right]}^{-1}$$

(67)

where the details of $[U_Z,V_Z,W_Z{]}^T$ and $[{\stackrel{·}{U}}_Z,{\stackrel{·}{V}}_Z,{\stackrel{·}{W}}_Z{]}^T$ are given in equations (16) and (17) respectively.

The affine nonlinear dynamic equations have been established. Based on equations (64)–(67), the influence of abrupt centroid shift on the motion of HSV is divided into three aspects. That is to say, the values of centroid shift are decomposed into three axes, $x-y-z$ . The details are organized as follows: case 1 describes the centroid shift just along with $x-\mathrm{axis}$ , case 2 describes the centroid shift just along with $y-\mathrm{axis}$ , and case 3 describes the centroid shift just along with $z-\mathrm{axis}$ .

Case 1: the centroid shift just along with x-axis

In this case, we assume that centroid shift moves along $x-\mathrm{axis}$ , and the values of centroid shift, $[\Delta x,\Delta y,\Delta z{]}^T$ are selected as the following five groups: $[0,0,0{]}^T,[0.08,0,0{]}^T,[0.04,0,0{]}^T,[-0.04,0,0{]}^T,[-0.08,0,0{]}^T$ , remarked as $({\alpha }_i,{\beta }_i,{\mu }_i,p_i,q_i,r_i),i=1,2,3,4,5$ , respectively. Due to the centroid shift moving along $x-\mathrm{axis}$ , the symmetry of HSV is still not destroyed, therefore, the products of inertia are still zeros, $I_{\mathrm{xy}}=I_{\mathrm{yx}}=I_{\mathrm{yz}}=I_{\mathrm{zy}}=0$ . Then, by substituting those assumptions and the values of centroid shift into equations (66) and (67), we have

$$\Delta f_f=m{I}^{-1}\left(\left[\begin{array}{c}0\\ {\stackrel{.}{W}}_Z\\ -{\stackrel{.}{V}}_Z\end{array}\right]+\left[\begin{array}{c}0\\ p{V}_Z-q{U}_Z\\ p{W}_Z-r{U}_Z\end{array}\right]+\left[\begin{array}{c}0\\ -g\cos \gamma \cos \mu \\ g\cos \gamma \sin \mu \end{array}\right]\right)\Delta x,\Delta g_f=0_{3\times 3}$$

(68)

where the details of $[U_Z,V_Z,W_Z{]}^T$ and $[{\stackrel{·}{U}}_Z,{\stackrel{·}{V}}_Z,{\stackrel{·}{W}}_Z{]}^T$ are given in equations (16) and (17), respectively.

By analyzing equation (68), we have that the effect of centroid on the motion of HSV in roll rate $(p)$ is zero. Then, referring to equations (62) and (63), by substituting the assumptions and values of centroid shift, we have

$$\begin{array}{ccc}\begin{array}{c}I{\prime }_{\mathrm{xx}}=I_{\mathrm{xx}}+m((\Delta y{)}^2+(\Delta z{)}^2)=I_{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{I{\prime }_{\mathrm{yy}}=I_{\mathrm{yy}}+m((\Delta x{)}^2+(\Delta z{)}^2)=I_{\mathrm{yy}}+m(\Delta x{)}^2}\\ \multicolumn{1}{c}{I{\prime }_{\mathrm{zz}}=I_{\mathrm{zz}}+m((\Delta x{)}^2+(\Delta y{)}^2)=I_{\mathrm{zz}}+m(\Delta x{)}^2}\end{array}\hfill & ,\hfill & \begin{array}{c}{I\prime }_{\mathrm{xz}}=I_{\mathrm{xz}}+m\Delta x\Delta z=I_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{{I\prime }_{\mathrm{xy}}=0}\\ \multicolumn{1}{c}{{I\prime }_{\mathrm{yz}}=0}\end{array}\hfill \end{array}$$

(69)

where $I{\prime }_{(·)}$ represents the moment of inertia and the product of inertia after centroid shift occurrence.

As for the variation of moment applied to HSV due to centroid shift, referring to equation (60), the details are given as follows

$$\begin{array}{c}{l}_Z=\hat{q}{\mathrm{bSC}}_l^{*}+Y_Z\Delta z-Z_Z\Delta y=\hat{q}{\mathrm{bSC}}_l^{*}\hfill \\ m_Z=m_{\mathrm{mrc}}-X_{\mathrm{cg}}Z-X_Z\Delta z+Z_Z\Delta x=m_{\mathrm{mrc}}-X_{\mathrm{cg}}Z+Z_Z\Delta x\hfill \\ n_Z=n_{\mathrm{mrc}}-X_{\mathrm{cg}}Y+X_Z\Delta y-Y_Z\Delta x=n_{\mathrm{mrc}}-X_{\mathrm{cg}}Y-Y_Z\Delta x\hfill \end{array}$$

(70)

where $X_Z,Y_Z,Z_Z$ are seen in equations (91)–(93), respectively. By analyzing equation (70), we can conclude that the centroid shift $\Delta x$ has no influence on the roll moment, $l_Z$ , of HSV.

In this article, we just consider the centroid shift on the motion of HSV and assume that the aerodynamic configuration of the aircraft has not changed; therefore, the equations of aerodynamic force in such simulation remain the same as previous, seen in equation (63). Some other initial states of HSV are given: $V=3\mathrm{km}/\mathrm{s}$ , $h=30\mathrm{km}$ , $\chi =0\circ ,\gamma =0\circ$ , $\alpha =2\circ ,\beta =0\circ ,\mu =0\circ$ , and $p=q=r=0\mathrm{rad}/\mathrm{s}$ . The simulation results are shown in Figure 4.

Figure 4.

The simulation results of HSV with centroid shift along with $x-\mathrm{axis}$ : response curves of (a) $\alpha$ , (b) $\beta$ , (c) $\mu$ , (d) p, (e) q, and (f) r.

As shown in Figure 4, we just move center of mass along $x-\mathrm{axis}$ . In Figure 4(d)–(f), combined with equation (70), we can see that the centroid shift along $x-\mathrm{axis}$ has little effect on r and p, and much larger influence on q in each group, and the same situation occurs to the attitude response curves of HSV. The variation of $\alpha$ is much larger than that of $\beta ,\mu$ , as shown in Figure 4(a)–(c). Taking the response curves of q, for example, the red curve, $q_1$ , is taken as the reference curve. The further the center of mass moves, compared with $q_2,q_3$ , the larger the pitching moment exerts on pitch channel of HSV. The same situation occurs to the opposite direction of centroid shift, seen in $q_4,q_5$ . Besides, the coupling of the system is increasing with the growth of centroid shift by comparing ${\beta }_4,{\mu }_4,p_4,r_4$ with ${\beta }_5,{\mu }_5,p_5,r_5$ in Figure 4. In conclusion, centroid shift moving along with $x-\mathrm{axes}$ has a major impact on $\alpha ,q$ , and both $\beta ,\mu$ and $p,r$ are affected. Those influences are irregular and changeable.

Case 2: the centroid shift just along with y-axis

In this case, we suppose that centroid shift moves just along $y-\mathrm{axis}$ , $[\Delta x,\Delta y,\Delta z{]}^T$ are still given in the following five groups: $[0,0,0{]}^T,[0,0.015,0{]}^T,[0,0.05,0{]}^T,[0,-0.015,0{]}^T,[0,-0.05,0{]}^T$ , remarked as $({\alpha }_i,{\beta }_i,{\mu }_i,p_i,q_i,r_i),i=1,2,3,4,5$ , respectively. The other parameters and initial conditions are the same with case 1. By substituting the data and the assumption into equations (66) and (67) again, we have

$$\Delta f_f=I^{-1}\left[\begin{array}{c}\left(q^2-r^2\right)I_{yz}-pr{I}_{xy}\\ qr{I}_{xy}-pq{I}_{yz}\\ \left(p^2-q^2\right)I_{xy}+pr{I}_{yz}\end{array}\right]+m{I}^{-1}\left(\left[\begin{array}{c}-{\stackrel{.}{W}}_{{\rm Z}}\\ 0\\ {\stackrel{.}{U}}_{{\rm Z}}\end{array}\right]+\left[\begin{array}{c}-p{V}_{{\rm Z}}+q{U}_{{\rm Z}}\\ 0\\ q{W}_{{\rm Z}}-r{V}_{{\rm Z}}\end{array}\right]+\left[\begin{array}{c}g\cos \gamma \cos \mu \\ 0\\ g\sin \gamma \end{array}\right]\right)\Delta y$$

(71)

where $\Delta g_f$ is seen in equation (67).

According to equations (62)–(63), we have

$$\begin{array}{c}{I\prime }_{\mathrm{xx}}=I_{\mathrm{xx}}+m({(\Delta y)}^2+{(\Delta z)}^2)=I_{\mathrm{xx}}+m{(\Delta y)}^2\hfill \\ \multicolumn{1}{c}{{I\prime }_{\mathrm{yy}}=I_{\mathrm{yy}}+m({(\Delta x)}^2+{(\Delta z)}^2)=I_{\mathrm{yy}}}\\ \multicolumn{1}{c}{{I\prime }_{\mathrm{zz}}=I_{\mathrm{zz}}+m({(\Delta x)}^2+{(\Delta y)}^2)=I_{\mathrm{zz}}+m{(\Delta y)}^2}\end{array},\begin{array}{c}{I\prime }_{\mathrm{xz}}=I_{\mathrm{xz}}+m\Delta x\Delta z=I_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{{I\prime }_{\mathrm{xy}}=I_{\mathrm{xy}}+m\Delta x\Delta y=I_{\mathrm{xy}}}\\ \multicolumn{1}{c}{{I\prime }_{\mathrm{yz}}=I_{\mathrm{yz}}+m\Delta y\Delta z=I_{\mathrm{yz}}}\end{array}$$

(72)

where $I{\prime }_{(·)}$ represents the moment of inertia and the product of inertia after centroid shift occurrence.

As for the variation of moment on HSV, according to equation (70), we have

$$\begin{array}{c}{l}_Z=\hat{q}{\mathrm{bSC}}_l^{*}+Y_Z\Delta z-Z_Z\Delta y=\hat{q}{\mathrm{bSC}}_l^{*}-Z_Z\Delta y\hfill \\ m_Z=m_{\mathrm{mrc}}-X_{\mathrm{cg}}Z-X_Z\Delta z+Z_a\Delta x=m_{\mathrm{mrc}}-X_{\mathrm{cg}}Z\hfill \\ n_Z=n_{\mathrm{mrc}}-X_{\mathrm{cg}}Y+X_Z\Delta y-Y_Z\Delta x=n_{\mathrm{mrc}}-X_{\mathrm{cg}}Y+X_Z\Delta y\hfill \end{array}$$

(73)

where $X_Z,Y_Z,Z_Z$ are seen in equations (91)–(93).

Referring to case 1, we also assume that the aerodynamic configuration of the aircraft has not been damaged, that is to say, the aerodynamic force applied to HSV is still the same as the previous, seen in equation (63). Some initial states of HSV are also the same with case 1 and the simulation results are shown in Figure 5.

Figure 5.

The simulation results of HSV with centroid shift along with $y-\mathrm{axis}$ : response curves of (a) $\alpha$ , (b) $\beta$ , (c) $\mu$ , (d) p, (e) q, and (f) r.

Comparing with equations (72) and (69), the variation of the moment of inertia complies with the following law. If the centroid shift just moves along with $y-\mathrm{axis}$ , it has no effect on $I_{\mathrm{yy}}$ , and the same situation occurs in $z-\mathrm{and}y-\mathrm{axes}$ . $Y\prime$ in the curves of Figure 5 means that the centroid shift just moves along with $y-\mathrm{axis}$ . From the comparison between Figure 5(d)–(f) and Figure 4(d)–(f), on the same guidance indexes, we have that the center of mass moving along with $y-\mathrm{axes}$ gives rise to the much strong coupling of system, increasing to 10 times than that of $x-\mathrm{axis}$ . As for the attitude angles of HSV, the centroid shift along with $y-\mathrm{axis}$ plays a major role in the rolling motion of HSV by comparing Figure 5(c) and Figure 5(a) and (b). Besides, the larger the centroid shift moves along with $y-\mathrm{axes}$ , the more intense the roll angle changes, as shown in ${\mu }_4,{\mu }_5$ in Figure 5(c). Although the influence of centroid shift on $\alpha ,\beta ,q,r$ is in secondary status, it still exerts a great impact on the motion of HSV and should be taken into consideration.

Case 3: the centroid shift just along with z-axis

In this case, we assume that centroid shift moves along with $z-\mathrm{axis}$ . The values of centroid along with $z-\mathrm{axis}$ are also divided into five groups, namely, $[\Delta x,\Delta y,\Delta z{]}^T:[0,0,0{]}^T,[0,0.015,0{]}^T,[0,0.05,0{]}^T,[0,-0.015,0{]}^T,[0,-0.05,0{]}^T$ , which are remarked as $({\alpha }_i,{\beta }_i,{\mu }_i,p_i,q_i,r_i),i=1,2,3,4,5$ , respectively. Analogous to cases 1 and 2, and we have

$$\Delta f_f=I^{-1}\left[\begin{array}{c}\left(q^2-r^2\right)I_{yz}-pr{I}_{xy}\\ qr{I}_{xy}-pq{I}_{yz}\\ \left(p^2-q^2\right)I_{xy}+pr{I}_{yz}\end{array}\right]+m{I}^{-1}\left({\left[\begin{array}{c}{\stackrel{.}{V}}_{{\rm Z}}\\ -{\stackrel{.}{U}}_A\\ 0\end{array}\right]}_{{\rm Z}}+\left[\begin{array}{c}-p{W}_{{\rm Z}}+r{U}_{{\rm Z}}\\ -q{W}_{{\rm Z}}+r{V}_{{\rm Z}}\\ 0\end{array}\right]+\left[\begin{array}{c}-g\cos \gamma \sin \mu \\ -g\sin \gamma \\ 0\end{array}\right]\right)\Delta z$$

(74)

where $\Delta g_f$ is seen in equation (67).

By referring to equations (62) and (63), we have

$$\begin{array}{c}{I\prime }_{\mathrm{xx}}=I_{\mathrm{xx}}+m({(\Delta y)}^2+{(\Delta z)}^2)=I_{\mathrm{xx}}+m{(\Delta z)}^2\hfill \\ \multicolumn{1}{c}{{I\prime }_{\mathrm{yy}}=I_{\mathrm{yy}}+m({(\Delta x)}^2+{(\Delta z)}^2)=I_{\mathrm{yy}}+m{(\Delta z)}^2}\\ \multicolumn{1}{c}{{I\prime }_{\mathrm{zz}}=I_{\mathrm{zz}}+m({(\Delta x)}^2+{(\Delta y)}^2)=I_{\mathrm{zz}}}\end{array},\begin{array}{c}{I\prime }_{\mathrm{xz}}=I_{\mathrm{xz}}+m\Delta x\Delta z=I_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{{I\prime }_{\mathrm{xy}}=I_{\mathrm{xy}}+m\Delta x\Delta y=I_{\mathrm{xy}}}\\ \multicolumn{1}{c}{{I\prime }_{\mathrm{yz}}=I_{\mathrm{yz}}+m\Delta y\Delta z=I_{\mathrm{yz}}}\end{array}$$

(75)

As for the variation of moment on HSV, referring to equation (90), we have

$$\begin{array}{c}{l}_Z=\hat{q}{\mathrm{bSC}}_l^{*}+Y_Z\Delta z-Z_Z\Delta y=\hat{q}{\mathrm{bSC}}_l^{*}+Y_Z\Delta z\hfill \\ m_Z=m_{\mathrm{mrc}}-X_{\mathrm{cg}}Z-X_Z\Delta z+Z_Z\Delta x=m_{\mathrm{mrc}}-X_{\mathrm{cg}}Z-X_Z\Delta z\hfill \\ n_Z=n_{\mathrm{mrc}}-X_{\mathrm{cg}}Y+X_Z\Delta y-Y_Z\Delta x=n_{\mathrm{mrc}}\hfill \end{array}$$

(76)

where $X_Z,Y_Z,Z_Z$ are seen in equations (91)–(93).

Referring to cases 1 and 2, both the assumptions and the initial states of HSV are the same, and the simulation results are shown in Figure 6.

Figure 6.

The simulation results of HSV with centroid shift along with $z-\mathrm{axis}$ : response curves of (a) $\alpha$ , (b) $\beta$ , (c) $\mu$ , (d) p, (e) q, and (f) r.

By analyzing the simulation results, the coupling of HSV due to centroid along with $z-\mathrm{axis}$ locates in the middle of cases 1 and 2, which can be clearly seen in Figure 6(b) and (c). By comparing ${\mu }_4,{\mu }_5$ of Figure 6(c) as well as $p_4,q_4,r_4$ and $p_5,q_5,r_5$ of Figure 6(d)–(f), we know that the more the centroid shift moves, the stronger is the influence on the motion of HSV, just as the same condition with cases 1 and 2. In all three cases, we know that the effect of centroid shift on the motion of HSV is coupling and possesses the different characteristics when the center of mass moves along with different axes, as shown in Figures 5 and 6. Centroid shift with $x-\mathrm{axis}$ mainly affects the motion of pitching movement. Centroid shift with $y-\mathrm{axis}$ mainly affects the motion of rolling movement; however, the characteristics of centroid shift with $z-\mathrm{axis}$ are not very clear, but it gives rise to a much strong response of $\beta ,\mu$ and greatly enhances the coupling of the system, shown in Figure 6(b) and (c).

Example: dynamic response of HSV with centroid shift

It is well known that the controller design is always following two design methods: (1) resting on the robustness of controller itself and (2) designing the auxiliary system to compensate those effects to guarantee the performance of HSV. First, the robustness of controller is reflected in itself. The TSMC is selected due to its strong robustness.^21–24 For the other, referring to the previous studies,^24–27 RBF neural network (RBFNN) observer is applied to estimate the whole uncertainties and disturbances caused by centroid shift.

The design of TSMC and RBFNN

The purpose of control is to track ${\omega }_c$ accurately and fast in finite time. In this article, referring to the previous studies,^21,22,24 the terminal sliding mode surface is designed as follows

$$s={\omega }_e+\underset{0}{\overset{t}{\int }}(a_1{\omega }_e+b_1{{\omega }_e}^{q_1/p_1})d\tau$$

(77)

where ${\omega }_e=\omega -{\omega }_c$ represents the tracking error, and $a_1,b_1,p_1,q_1$ are the design parameters, satisfying $a_1>0,b_1>0,p_1>q_1>0$ ; in addition, $p_1$ and $q_1$ are positive odd numbers.

Then, by derivating both sides of equation (77), we have

$$\stackrel{.}{s}={\stackrel{.}{\omega }}_e+a_1{\omega }_e+b_1{\omega }_e^{q_1/p_1}=f_f+\Delta f_f+(g_f+\Delta g_f)M_c-{\stackrel{.}{\omega }}_c+a_1{\omega }_e+b_1{\omega }_e^{q_1/p_1}$$

(78)

Order $\phi (\omega ,u,t)=\Delta f_f+\Delta g_f{M}_c$ , and equation (78) can be simplified as follows

$$\stackrel{·}{s}=f_f+g_f{M}_c-{\stackrel{·}{\omega }}_c+a_1{\omega }_e+b_1{{\omega }_e}^{q_1/p_1}+\phi (\omega ,u,t)$$

(79)

Then, RBFNN control is applied to estimate $\phi (\omega ,u,t)$ online and the details are obtained as follows.

Considering a class of systems with the input of $x\in R^n$ , when x belongs to a compact set ${\Pi }_x$ , the output of the neural network can be expressed as follows^24–27

$$u_{\mathrm{nn}}={\hat{W}}^T\varphi (x,\hat{\xi })$$

(80)

where $\hat{W}=[w_{\mathrm{ij}}{]}_{l\times m}$ is the weight between the hidden layer and the output layer. l and m denote the number of the hidden layer neurons and the network output dimension, respectively. $\xi =[{\xi }_1^T,{\xi }_2^T,\dots ,{\xi }_l^T{]}^T$ represents the central vector of gauss function, namely, $\varphi (x,\xi )={[{\varphi }_1,{\varphi }_2,\dots ,{\varphi }_l]}^T$ , where ${\varphi }_i$ is the Gauss basis functions, the specifics are as follows

$${\varphi }_i(x,\xi )=\exp \left[\frac{-{(x-{\xi }_i)}^T(x-{\xi }_i)}{{\eta }^2}\right]$$

(81)

where $\eta$ is the width of the Gauss function. It has been confirmed that RBFNN can approximate the bounded continuous function $\phi (\omega ,u,t)$ with arbitrary precision in the compact set ${\Pi }_x$ , namely

$$\phi (\omega ,u,t)=W^{*T}\varphi (x,{\xi }^{*})+{\epsilon }_f(x)$$

(82)

where ${\epsilon }_f(x)$ is the approximation error, meeting $||{\epsilon }_f(x)||\le \epsilon ,\epsilon >0$ . $W^{*}$ is the optimal weight and ${\xi }^{*}$ is the center of the Gauss basis function; besides, $W^{*}$ is defined as follows

$$W^{*}=\arg \underset{\hat{W}\in {\Pi }_x}{min}[\underset{x\in {\Pi }_x}{sup}\Vert f\left(x\right)-f(x|\hat{W},{\xi }^{*})\Vert ]$$

(83)

Then, the robust adaptive control law based on RBFNN compensation is designed as follows

$$M_c=g^{-1}(-f_f+{\stackrel{·}{\omega }}_c-a_1{\omega }_e-b_1{\omega }_e^{q_1/p_1}-u_{\mathrm{ns}}-\rho \mathrm{sign}(s))$$

(84)

where $u_{\mathrm{ns}}=u_{\mathrm{nn}}+u_s$ , $u_s$ is the error compensation term to improve the compensation error of RBFNN. $\rho$ is the design parameter, meeting $\rho >0$ . In order to handle the problem of chattering caused by sliding mode control, a continuous function is used to replace the symbolic function

$$\mathrm{sign}(s)=\frac{s}{\Vert s\Vert +\Delta +\Vert {\omega }_e\Vert }$$

(85)

where s is the sliding mode surface, $\Delta$ is a small enough number, satisfying $\Delta >0$ . ||.|| represents the 2-norm.

For the limitation of the length of the article, the convergence time of TSMC and the adaptive law of RBF are just provided without the proof. The more details are seen in previous studies.^24–27 The convergence time of TSMC is offered as follows

$$t_{\mathrm{ri}}=\ln ({\omega }_{\mathrm{ei}}{(0)}^{1-q_1/p_1}+\frac{b_1}{a_1})\frac{p_1}{a_1(q_1-p_1)}$$

(86)

where $t_{\mathrm{ri}}$ refers to the convergence time. The structure of the closed-loop control system is shown in Figure 7.

Figure 7.

The TSMC control structure with RBFNN of the HSV.

The weight of RBFNN, the online update law of the center, and the adaptive control law for approximation error compensation are given as follows^24–27

$$\stackrel{·}{\hat{W}}={\tau }_1\hat{\varphi }{s}^T,\stackrel{·}{\hat{\xi }}={\tau }_2(s^T{\hat{W}}^T\varphi {\prime }_{\xi }{)}^T,\stackrel{·}{\hat{\lambda }}={\tau }_3\Vert s^T\Vert$$

(87)

where ${\tau }_1,{\tau }_2,{\tau }_3$ are the regulatory factors, meeting ${\tau }_1,{\tau }_2,{\tau }_3>0$ .

In order to guarantee the bound of estimated parameters and avoid parameter drift, we constrain $\hat{W}$ , $\hat{\xi }$ , and $\hat{\lambda }$ , respectively, to reside inside compact sets ${\Omega }_W$ , ${\Omega }_{\xi }$ , and ${\Omega }_{\lambda }$ , which are defined as follows

$${\Omega }_W=\{\widehat{W}:\underset{\_}{\widehat{W}}\le \widehat{W}\le \overline{\widehat{W}}\},{\Omega }_{\xi }=\{\widehat{\xi }:\underset{\_}{\widehat{\xi }}\le \widehat{W}\le \overline{\widehat{\xi }}\},{\Omega }_{\lambda }=\{\widehat{\lambda }:\underset{\_}{\widehat{\lambda }}\le \widehat{W}\le \overline{\widehat{\lambda }}\}$$

(88)

where $\underset{¯}{\hat{W}},\overline{\hat{W}},\underset{¯}{\hat{\xi }},\overline{\hat{\xi }}$ and $\underset{¯}{\hat{\lambda }},\overline{\hat{\lambda }}$ are the design parameters, then the adaptive law is modified as follows

$$\stackrel{.}{\widehat{W}}=\mathrm{Proj}(\widehat{W},{\tau }_1\widehat{\varphi }{s}^T),\stackrel{.}{\widehat{\xi }}=\mathrm{Proj}(\widehat{\xi },{\tau }_2{(s^T{\widehat{W}}^T{{\varphi }^{\prime }}_{\xi })}^T),\stackrel{.}{\widehat{\lambda }}=\mathrm{Proj}(\widehat{\lambda },{\tau }_3|\left|s^T\right||)$$

(89)

where $\mathrm{Proj}(a,b)$ stands for $a:[\hat{W},\hat{\xi },\hat{\lambda }],b:[\stackrel{·}{\hat{W}}={\tau }_1\hat{\varphi }{s}^T,\stackrel{·}{\hat{\xi }}={\tau }_2(s^T{\hat{W}}^T\varphi {\prime }_{\xi }{)}^T,\stackrel{·}{\hat{\lambda }}={\tau }_3||s^T||$ , and $\mathrm{Proj}(a,b)=\{\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{b}\end{array}\begin{array}{c}\mathrm{if}a=\underset{¯}{\Theta },\mathrm{and}b<0\hfill \\ \multicolumn{1}{c}{\mathrm{if}a=\overline{\Theta },\mathrm{and}b>0}\\ \multicolumn{1}{c}{}\end{array}$ is a discontinuous projection operator proposed in the works by Lian et al. ²⁸ and Tang et al., ²⁹ where $\underset{¯}{\Theta }:[\underset{¯}{\hat{W}},\underset{¯}{\hat{\xi }},\underset{¯}{\hat{\lambda }}]$ , $\overline{\Theta }:[\overline{\hat{W}},\overline{\hat{\xi }},\overline{\hat{\lambda }}]$ .

Simulation

Considering that the HSV is working on a hypersonic reentry flight with the speed of $3\mathrm{km}/\mathrm{s}$ , and the flight height is $V=3\mathrm{km}/\mathrm{s}$ , $\chi =0\circ ,\gamma =0\circ$ , $\alpha =2\circ ,\beta =0\circ ,\mu =0\circ$ , $p=q=r=0\mathrm{rad}/\mathrm{s}$ are the initial attitudes and angular velocities. The reference signals are $p_c=0\circ ,r_c=0\circ$ , and $q_c$ is the square wave signal with period 10 s and amplitude 0.05 rad/s. We assume that the centroid shift occurs at the 5 s and the variation meets $\overline{\rho }:[0.03,0.015,0.025{]}^{T}m$ . Based on such assumption and initial states of HSV, the simulations are conducted as follows.

At first, the TSMC combined with RBFNN observer is applied to handle the problem. The concrete parameters are as follows: $a_1=1,b_1=1,q_1=7,p_1=9,\rho =2,\Delta =0.5$ . The parameters about RBFNN are as follows: $l=5$ , $m=3$ , ${\tau }_1=0.3,{\tau }_2=0.4,{\tau }_3=0.3$ , ${\eta }_p=0.8,{\eta }_q=1.2,{\eta }_r=0.8$ , ${\xi }_p=[-0.8,-0.4,0,0.4,0.8]$ , ${\xi }_q=[-0.2,-0.1,0,0.1,0.2]$ , ${\xi }_r=[-0.08,-0.04,0,0.04,0.08]$ . Besides, taking the robustness of observer into account, a sliding mode disturbance observer (SMDO) is adopted. With the same conditions as the former, the index order is given by $p=r=0\mathrm{rad}/\mathrm{s}$ , q is the square signal with the period 10 s and the amplitude is 0.05 rad/s, ${\varphi }_f=\mathrm{diag}\{0.35,0.20,0.15\}$ , ${\varphi }_f$ represents the bounds of some parameters of SMDO needed in this article when applying the controller of TSMC combined with SMDO, ${\delta }_p=0.015$ and ${\delta }_q=0.018,{\delta }_r=0.015$ are the small numbers and the function is described in equation (85). The simulation results are as shown in Figure 8.

Figure 8.

Simulation results of the TSMC with RBFNN: tracking effects of (a) p, (b) q, and (c) r and deflection angles of (d) ${\delta }_{\mathrm{LE}}$ , (e) ${\delta }_{\mathrm{RE}}$ , and (f) ${\delta }_{\mathrm{RUD}}$ .

From the simulations shown in Figure 8, $(·{)}_{\mathrm{nml}}$ is the index signal and $(·{)}_{\mathrm{rbf}}$ represents the simulation results under TSMC with RBFNN, while $(·{)}_{\mathrm{SMDO}}$ is the simulation result under TSMC with SMDO. ${\delta }_{\mathrm{LE}},{\delta }_{\mathrm{RE}},{\delta }_{\mathrm{RUD}}$ represent the responses of deflection angles, as shown in Figure 8(d)–(f).

The tracking effect shown in Figure 8(a)–(c) indicates that HSV with centroid shift is controllable to a certain degree under the controller of TSMC combined with RBFNN or SMDO. From the tracking effect of q in Figure 8(b), TSMC combined with RBFNN has a better tracking precision, while TSMC with SMDO has a smaller overshoot, and the same situation occurs in both Figure 8(a) and (c). The reason may be related to the sensitivity of the observer to the effects of centroid shift. As for the response rate, the TSMC with SMDO is better because SMDO with the fewer adaptive parameters and small computational can amount. From the numerical simulation results in Figures 4–6, the effects of centroid shift can be mapped some extra torque on HSV. Therefore, once the centroid shift occurs, it is bound to cause the shaking of HSV body. Then, how to stabilize HSV rapidly is crucial. From the deflection angle curves Figures 4(d)–(f) to 6(d)–(f), it can be seen that the effects of some certain centroid shift on the motion of HSV are still in the controllable range according to the simulation results in this article.

The former research works are focused on the compensation, and the other control strategy is up to the robustness of controller itself. Thus, the TSMC is adopted under the same simulation conditions, $a_1=1,b_1=1,q_1=7,p_1=9$ mentioned above, and the simulation results are shown in Figure 9.

Figure 9.

Simulation results of the TSMC: tracking effects of (a) p, (b) q, and (c) r and deflection angles of (d) ${\delta }_{\mathrm{LE}}$ , (e) ${\delta }_{\mathrm{RE}}$ , and (f) ${\delta }_{\mathrm{RUD}}$ .

As shown in Figure 9, $(·{)}_{\rho }$ represents the controller with different levels of robustness, where $\rho =0.5,2,5$ . ${\delta }_{\mathrm{LE}},{\delta }_{\mathrm{RE}},{\delta }_{\mathrm{RUD}}$ still represent the response of control deflection, as shown in Figure 9(d)–(f). From the simulation results of Figure 9(a)–(c), under the identical control conditions, it is obvious that the smaller $\rho$ (from 0.5 to 2 and then up to 5) and the smaller overshoot, shown in ( $p_{0.5},p_2,p_5$ ) Figure 9(a). The reason may be that the center shift is a smaller offset and the deflection angles of ${\delta }_e,{\delta }_a,{\delta }_r$ are all in an acceptable range, as shown in Figure 9(d)–(f). However, if the centroid shift is large, abrupt, and irregular, further research works are needed.

Conclusion

This article proposes a thorough research for the influence of centroid shift on the motion of HSV. First, based on Newton’s laws and Varignon’s theory, a mathematic model for HSV with centroid shift is built up. Then, a zero-input response test is conducted to analyze the effect of centroid shift on the motion of HSV. The results show that the unexpected centroid shift exerts an unpredictable and immeasurable torque to HSV, resulting in the strong coupling and uncertainties to HSV, which poses a great challenge to the controller. In this article, TSMC combined with RBFNN observer and TSMC itself alone is applied to handle such problems in view of robust control itself and auxiliary compensation. The unexpected centroid shift can be inhibited in a small region, but when the centroid shift is larger, irregular, and time-variant, the control effect is hard to predict. This article is just the beginning, and a lot of further research works are still needed, especially in the case of faults.