Longitudinal Trim and Stability Analysis of Hybrid Airship with Suspended Payload for Single Body and Two Body Dynamics Using Bifurcation Method

The proposed model, a winged hybrid small-sized airship flight vehicle with a suspended payload, appears in Figure 2 and will be used for modelling and analysis in this research. A winged aircraft, an airship, and aircraft elements known as the Plimp and parafoil payload delivery flight system was combined to make the model, which is in the first row and first column. Figure 3 depicts the workflow for creating airship nonlinear flight dynamics models and simulations. The associated payload is assumed to be carried by the hybrid airship.

2.png

Figure 2. Winged hybrid airship flight vehicle design concept with suspended payload [17]

3.png

Figure 3. Flow chart for mathematical modelling [17]

Table 1. Complexity chart in mathematical derivation

In comparison to an aircraft or a parafoil delivery mode, the mathematical modelling of a hybrid airship with payload is more challenging, as it indicates in the above Table 1. Due to the fact that the centres of mass or gravity and volume are situated at two different locations, as indicated in Figure 4, the mathematical derivation of a hybrid airship is challenging. At the centre of mass or centre of gravity, it is simple to obtain the moment equation and force equation. In a hybrid airship, the centre of volume is almost exactly in the middle of the hull, which is filled with helium or hydrogen gas. In a hybrid airship, the payload's attached gondola is close to the centre of mass, or centre of gravity.

4.jpg

Figure 4. Coordinate location of CG and CV [18]

2.1 3-DOF single body dynamic longitudinal model

A Hybrid Airship is modelled and simulated using assumptions as the flat Earth approximation, airship considered as rigid body, mass and volume of hybrid airship assumed as constant, centre of volume and centre of buoyancy is assumed to coincide, and atmosphere is assumed as stationary. For study or analysis purpose international atmospheric conditions are followed for flight dynamic analysis. The whole volume of airship under study is usually filled with lighter than air helium gas, that’s why centre of buoyancy is assumed to coincide with centre of volume. Equations of motion are established using Newton's and Euler's equations for general aircraft in addition of the added mass and inertia effects and buoyancy force. The geometric specifications and data of the vehicle is given in Tables 2 and 3.

Table 2. Geometry data of a HULL [18]

Table 3. Payload geometry [18]

The geometric view of the vehicle with specifications is shown in Figure 5. Figure 6 shows the aerodynamic data plot of airship taken from researches [2, 18].

The complete Non-linear 6 DOF equation of motion of a Hybrid airship is given in Eq. (1).

$\left\lceil\begin{array}{cc}m I_3+M^{\prime} & -m r_G^X \\ m r_G^X & I_0+I_0^{\prime}\end{array}\right]\left[\begin{array}{c}\dot{V}_0 \\ \dot{\omega}_0\end{array}\right]+\left[\begin{array}{c}m\left(\omega_0 \times V_0+\omega_0 \times\left(\omega_0 \times r_G\right)\right) \\ \omega_0 \times\left(I_0 \omega_0\right)+m r_G \times\left(\omega_0 \times V_0\right)\end{array}\right]=\left[\begin{array}{l}F \\ \zeta\end{array}\right]$                              (1)

where,

$r_G^X=\left[\begin{array}{ccc}0 & -a_z & a_y \\ a_z & 0 & -a_x \\ -a_y & a_x & 0\end{array}\right]$                           (2)

$\left(r_G^X\right)^T=-\left(r_G^X\right)$                       (3)

The complete Non-linear 6 DOF equation of motion for an Aircraft [19] is given in Eq. (4).

$\left[\begin{array}{cc}m I_3 & 0 \\ 0 & I_0\end{array}\right]\left[\begin{array}{l}\dot{V}_0 \\ \dot{\omega}_0\end{array}\right]+\left[\begin{array}{c}m\left(\omega_0 \times V_0\right) \\ \omega_0 \times\left(I_0 \omega_0\right)\end{array}\right]=\left[\begin{array}{l}F \\ \zeta\end{array}\right]$                         (4)

5.png

Figure 5. Top view of winged hybrid airship [18]

6a.jpg

(a) C_LV vs $\alpha$

6b.jpg

(b) C_DV vs $\alpha$

6c.jpg

(c) C_MV vs $\alpha$

Figure 6. Aerodynamic data [2, 18]

7.png

Figure 7. Flowchart for apparent mass calculation for hybrid airship

Addition of the apparent mass calculation in hybrid airship is shown in Figure 7. The hybrid airship is lighter than air vehicle and it consists of HULL filled with hydrogen/helium gas. So, apparent mass calculation, buoyancy force equation and mass distribution are important in the mathematical modelling of a hybrid airship. Contribution of mass components in hybrid airship is shown in Figure 8.

8.png

Figure 8. Mass distribution

The volume for hydrostatic lift provided in Eq. (5) and the Archimedes principle were used to calculate the lift in airship design (6).

$L_{h s t}=\operatorname{vol}\left(\rho_a-\rho_a\right)$                     (5)

where, ρ_a and ρ_g are the densities of air and lifting gas.

$w_{\text {net }}=\left(m_{\text {GTw }} \times g-L_{\text {hst }}\right)$                      (6)

The weight balanced by the hydrostatic lift is subtracted from the gross take-off weight to get the net weight, or w_net. The equation for the hybrid airship's propulsion system is given as Eq. (7a) to (7c).

$M_P=M_T=T d_z$                    (7a)

$X_P=X_T=T$                     (7b)

$Z_P=Z_T=0$                   (7c)

T is the total thrust produced by the engine, M_T is the moment owing to the propulsion system or thrust, and d_z is the distance between the propulsion system and the CV. The full nonlinear forces and moments equation for the 3DOF longitudinal dynamics is given in Eq. (8) to (11), which includes the axial force equation (U), normal force equation (W), pitching moment equation (q), and kinematic equation (θ ̇).Where X_a is axial aerodynamic force, X_S is axial buoyancy force, X_T is axial thrust force, Z_a is normal aerodynamic force, Z_S is buoyancy force in normal direction along z-axis, Z_T is normal thrust force, M_a is aerodynamic pitching moment along y-axis, M_s is buoyancy pitching moment due to buoyancy about CG and M_T is thrust moment along y-axis or pitching moment.

2.1.1 Non-linear longitudinal dynamic equations for hybrid airship

Force equations:

Axial force:

$m_x \dot{U}+\left(m a_x-X_{\dot{q}}\right) \dot{q}+m_z q W-m a_x q^2=X_a+X_S+X_T$                        (8)

Normal force:

$m_z \dot{W}-\left(m a_x+Z_{\dot{q}}\right) \dot{q}-m_x q U+m a_z q^2=Z_a+Z_S+Z_T$                     (9)

Moment equations:

Pitching moment:

$J_y \dot{q}+\left(m a_z+M_{\dot{u}}\right) \dot{U}-\left(m a_x+M_{\dot{w}}\right) \dot{W}+m a_x(q U)+m a_z(q W)=M_a+M_s+M_T$                     (10)

Kinematic equation:

$\dot{\theta}=q$                  (11)

In linearizing the nonlinear equation of motion from Eqns. (8) to (11), small perturbation theory is applied, where squares of variables and coupling derivatives are omitted and small angle assumption is used. An example of a complete expression of the aerodynamic force in x-axis is given in Eq. (12).

$X_a=X_1+X_u \cdot u+X_v \cdot v+X_w \cdot w+X_p \cdot p+X_q \cdot q+X_r \cdot r$                         (12)

where, X₁ is the trim equilibrium aerodynamic force which is defined as:

$X_1=\frac{1}{2} \rho_{a i r} V^2 s\left(C_L \sin \alpha-C_D \cos \alpha\right)$                       (13)

The rest of the components in the aerodynamic expression are the dynamic perturbations terms, which will be zero if the aircraft is in equilibrium. Some of the terms in the aerodynamic force and moment representation can be neglected if its contribution to the motion is very small [20, 21].

2.1.2 Linearized longitudinal dynamic equations for hybrid airship

Force equations:

Axial force:

$m_x \dot{U}+\left(m a_x-X_{\dot{q}}\right) \dot{q}+m_z q W_1=X_a+X_S+X_T$                           (14)

Normal force:

$m_z \dot{W}-\left(m a_x+Z_{\dot{q}}\right) \dot{q}-m_x q U_1=Z_a+Z_S+Z_T$                       (15)

Moment equations:

Pitching moment:

$J_y \dot{q}+\left(m a_z+M_{\dot{u}}\right) \dot{U}-\left(m a_x+M_{\dot{w}}\right) \dot{W}+m a_x(q U)+m a_z(q W)=M_a+M_s+M_T$                       (16)

Kinematic equation:

$\dot{\theta}=q$                     (17)

The expression of the aerodynamic force in x-axis shown in Eq. (12) is expressed below for longitudinal dynamics of the Airship.

$X_a=X_1+X_u \cdot u+X_w \cdot w+X_q \cdot q+X_{\delta_e} \cdot \delta_e$                          (18)

Table 4. Longitudinal stability derivatives

The expression for the stability derivatives is given in the Table 4. The aerodynamic coefficients (such as, $C_{m_\alpha}, C_{m_q}$ and $c_{m_{\delta_e}}$) are obtained either from wind tunnel test, numerical solutions, or semi-empirical approach [18-22]. For a complete analysis, the aerodynamic database that includes all the expression of the aerodynamic coefficients must be constructed for the vehicle. Since the hybrid airship is newly developed, such database has not been created. The main challenge here is that there is no source in the open literature that tackles the problem of estimating the aerodynamic coefficients of a winged hybrid airship, while there are several approximation methods for conventional aircrafts [23, 24]. With the assumption that the general equation for an aircraft is adequate for the winged airship application as well, the remaining aerodynamic coefficients are estimated using estimation methods and equations provided by studies [25, 26]. The estimates are only sufficient for this preliminary modelling and analysis of the winged hybrid airship. To ensure the accuracy of these aerodynamic coefficients, CFD and wind tunnel test must be performed on a model of the hybrid airship.

2.1.3 Buoyancy or static forces and moments:

The buoyancy forces and moments for hybrid airship is given in Eq. (19a) to (19c).

$X_s=-(m g-B) \sin \theta$                   (19a)

$Z_s=(m g-B) \cos \phi \cos \theta$                   (19b)

$M_s=-\left(m g a_z+B b_z\right) \sin \theta-\left(m g a_x+B b_x\right) \cos \phi \cos \theta$                   (19c)

3DOF hybrid airship Longitudinal translational and rotational dynamics Eqns. (8-11) of the model can be transformed and modelled into the Eqns. (37-44) using the wind axis representation and are given by:

$\dot{\mathrm{V}}=\left(\mathrm{T} \cos (\alpha)-\overline{\mathrm{q}} \mathrm{sC}_{\mathrm{D}}-(\mathrm{mg}-\mathrm{B}) \sin (\gamma)\right) /\left(\mathrm{m}-\mathrm{X}_{\mathrm{ud}}\right)$           (20)

$\dot{\gamma}=\left(T \sin (\alpha)+\overline{\mathrm{q}}_{\mathrm{S}} C_\mathrm{L}+\mathrm{z}_{\mathrm{q}}-(m g-B) \cos (\gamma)\right) /\left\{\left(\mathrm{m}-\mathrm{z}_{\mathrm{wd}}\right) \cdot \mathrm{V}\right\}$                      (21)

$\dot{\mathrm{q}}=\mathrm{M}_{\mathrm{z}} /\left(\mathrm{I}_{\mathrm{yy}}-\mathrm{Mq}_{\mathrm{d}}\right)$                           (22)

$\dot{\theta}=q$                   (23)

$\dot{R}=V \cos (\gamma)$                      (24)

$\dot{h}=V \sin (\gamma)$                        (25)

Considering α=θ-γ, So:

$\dot{\alpha}=q-\left[\left(T \sin (\alpha)+\overline{\mathrm{q}} s\mathrm{C}_1+\mathrm{z}_{\mathrm{q}}-(m g-B) \cos (\gamma)\right) /\left\{\left(\mathrm{m}-\mathrm{z}_{\mathrm{wd}}\right) \cdot \mathrm{V}\right\}\right.$                       (26)

$M_z=-M_s+\left(M q-m a_x u u-m a_z w w\right) \cdot q+\bar{q} S \bar{c} C_{m_{m p}}+T * 0 \cdot 1$                        (27)

2.2 4-DOF two body dynamic longitudinal model

To achieve the intended goal, a mathematical model of the hybrid airship is created by applying physical principles that describe the hybrid airship's non-linear dynamics and aerodynamics. Combining the 3DOF longitudinal dynamics of the hybrid airship with wing with the 3DOF equation of motion of the suspended payload, the 4DOF model of the hybrid airship with suspended payload has been suggested, as illustrated in Figure 9. A two-body dynamic 4DOF model of a hybrid airship with a hanging payload is presented in the research [17]. Multibody approach discussed in the research [27, 28] also. Thrust is equal to drag force in the airship dynamic model for level flight when the velocity is constant at 20 m/s. The effects of elevator at wing and moment at fin are also considered in the moment equation [17]. An apparent mass is also taken care in the modelling.

9.png

Figure 9. Longitudinal 4 DOF Two body dynamics for the hybrid airship payload system [7]

A 4-DOF dynamic model of a hybrid airship with suspended payload developed in the research [7] is shown below.

$\left[\begin{array}{cccccc}m_b \cos \theta_b & -m_b \sin \theta_b & m_b R_b & 0 & \cos \theta_b & -\sin \theta_b \\ m_b \sin \theta_b & m_b \cos \theta_b & 0 & 0 & \sin \theta_b & \cos \theta_b \\ 0 & 0 & I_b & 0 & -R_b \cos \theta_b & R_b \sin \theta_b \\ m_{p A} \cos \theta p & -m_{p A} \sin \theta p & 0 & m_{P A} R p \cos \mu & -\cos \theta_p & \sin \theta_p \\ m_{P C} \sin \theta_P & m_{P C} \cos \theta_P & 0 & -m_{p c} R_P \sin \mu & -\sin \theta_p & -\cos \theta_p \\ 0 & 0 & 0 & I_{P F} & z_{c p} & -x_{C p}\end{array}\right]\left[\begin{array}{c}\dot{u}_c \\ \dot{w}_C \\ \dot{q}_b \\ \dot{q}_p \\ F_x \\ F_z\end{array}\right]=\left\{\begin{array}{l}b_1 \\ b_2 \\ b_3 \\ b_4 \\ b_5 \\ b_6\end{array}\right\}$                       (28)

$b_1=-m_b g \sin \theta_b-\mathbb{Q}_b S_b C_{D_b} \cos \alpha_b$                    (29a)

$b_2=m_b g \cos \theta_b-Q_b S_b C_{D_b} \sin \alpha_b+m_b R_b q_b^2$                    (29b)

$b_3=0$                (29c)

$b_4=-m_p g \sin \theta_p+\mathbb{Q}_p S_p C_x+m_{p c} R_p \sin \mu q_p^2-(C-A)\left(u_C \sin \theta_P+w_C \cos \theta_p\right) q_p$                       (29d)

$b_5=m_P g \cos \theta p+Q_p S_\rho C_z+m_{p A} R_p \cos ^\mu q_P^2-(C-A)\left(U_C \cos \theta_P-w_C \sin \theta_p\right) q_p$                       (29e)

$b_6=-x_{p a} Q_p S_p C_z+Q_p S_{P C} C_M$                       (29f)

2.3 Aerodynamic model and longitudinal trim analysis

MATLAB is used to create subroutines for the model simulation. The lift force, drag force and pitching moment are expressed by L=qSC_L, D=qSC_Dand M=qcSC_m and the lift coefficient, drag coefficient and pitching moment coefficient are given by the Eqns. (30-32). These linear relations of aerodynamic coefficients given by Eqns. (30-32) are considered for the 3DOF longitudinal model simulation of hybrid airship and are also used to implement the bifurcation method with the elevator deflection $\delta_e$.

$C_L=\left\{C_{L_0}+C_{L_\alpha} \cdot \alpha+C_{L_q} \cdot \frac{q \bar{c}}{2 V_1}+C_{L_{\delta_e}} \cdot \delta_e\right\}$                   (30)

$C_D=\left\{C_{D_0}+C_{D_\alpha} \cdot \alpha+C_{D_q} \cdot \frac{q \bar{c}}{2 V_1}+C_{D_{\delta_e}} \cdot \delta_e\right\}$                          (31)

$C_m=\left\{C_{m_0}+C_{m_\alpha} \cdot \alpha+C_{m_q} \cdot \frac{q \bar{c}}{2 V_1}+C_{m_{\delta_e}} \cdot \delta_e\right\}$                       (32)

For the static longitudinal stability of the vehicle, the value of C_m must be negative, and the value of positive and negative C_m refer to nose up and nose down pitching moment respectively, C_m changes with the change of α.

The longitudinal nonlinear dynamic hybrid airship model presented with Eqns. (20-27) is used to obtain the solution for the trim or equilibrium state. These states are V, γ, q, θ and α, in which the states V, γ and θ are considered to be zero, i.e., considering the model to be flying with constant V, γ and θ values. In this situation the flight is in a straight path with the possible conditions as γ=0 for the level flight, γ=positive value for the ascend flight or γ=negative value for the descending flight, and at the same time with the constraint that the nose orientation is constant i.e., θ=constant. To achieve these trim conditions for the hybrid airship flight, we obtain the trim states by setting left-hand side of Eqns. (20-27) to zero and solving for the right-hand side, and this is given by Eqns. (33-37).

$0=\left(\mathrm{T} \cos (\alpha)-\overline{\mathrm{q}} s\mathrm{C}_{\mathrm{D}}-(\mathrm{mg}-\mathrm{B}) \sin (\gamma)\right) /\left(\mathrm{m}-\mathrm{X}_{\mathrm{ud}}\right)$                  (33)

$0=\left(T \sin (\alpha)+\overline{\mathrm{q}} \mathrm{sC}_{\mathrm{L}}+\mathrm{z}_{\mathrm{q}}-(m g-B) \cos (\gamma)\right) /\left\{\left(\mathrm{m}-\mathrm{z}_{\mathrm{wd}}\right) \cdot \mathrm{V}\right\}$                   (34)

$0=\mathrm{M}_{\mathrm{z}} /\left(\mathrm{I}_{\mathrm{yy}}-\mathrm{Mq}_{\mathrm{d}}\right)$                   (35)

$\dot{\theta}=0$                    (36)

The trim states are represented by (*) and are given by the following Eqns. (37-39).

$C_D^*=(\mathrm{T} \cos (\alpha)-(\mathrm{mg}-\mathrm{B}) \sin (\gamma)) / \bar{q} \mathrm{s}$                     (37)

$C_L^*=\left((m g-B) \cos (\gamma)-T \sin (\alpha)-\mathrm{z}_{\mathrm{q}}\right) /\bar{q} \mathrm{s}$                        (38)

$C_m^*=\left(M_s-\left(M q-m a_x u u-m a_z w w\right) \cdot q-T * 0 \cdot 1\right) / \bar{q} S c$                       (39)

From Eq. (37) we can write the thrust relation given by Eq. (40).

$\mathrm{T}=(\mathrm{mg}-\mathrm{B}) \sin (\gamma)+C_D^* \overline{\mathrm{q}}^* \mathrm{~s} / \cos (\alpha)$                        (40)

The initial values for the simulation are also used as the initial values to implement the bifurcation method.

The Bifurcation Method provides a possibility to enhance vehicle control design parameters by using the bifurcation approach to exhibit global stability. The method provides quantitative information on global stability and offers strategies to stabilize the nonlinear behavior of the aircraft. This method uses a continuation-based algorithm (CBA) to probe a nonlinear dynamic model using first-order ordinary differential equations (ODEs) and AUTO-07p [22] software to detect steady-state and various equilibrium points. ODE are represented as a function of states variables and control parameters given by, $\dot{x}=f f(x, u)$, the solution of the nonlinear equation is determined by, $\dot{x}=f f(x, u, p)=0$, by keeping the constraint p fixed with varying next parameter u, and this results in the finding of all trim states. The Standard Bifurcation Analysis (SBA) method is implemented with CBA. In this approach, the continuation algorithm computes all solutions with trim states of the system considering the varying and fixed parameters and determines the details of local stability at every trim state. The changes in the branch stability of the trim states solution are called or represented by the ‘Bifurcation Points.’ These ‘Bifurcation Points’ indicate the moment of the Eigenvalues leading to an unstable system in the stability plane i.e., the movement of the system Eigenvalues migrates from the left half to the right half of the complex plane, and prominently shows the unstable dynamical system occurrence or behaviour.

The method is implemented for the 3DOF longitudinal hybrid airship model by applying SBA for the nonlinear dynamic model [4, 7] and the state variables are given by x=[V, γ, α, θ, q, h]′. The states are further reduced for simplicity, considering constant velocity, V, constant height, h and considering very small changes in γ, thus the state variable considered for the analysis is reduced and is given by x=[α, θ, q]′. The SBA for the hybrid airship is implemented and bifurcation diagram with Eigen values for each case is obtained which shows the flight dynamics behaviours. The bifurcation diagram for velocity 20 m/s having mass of the vehicle as 10kg is shown in Figure 16. The stability information can be examined for the longitudinal dynamic model using Eigenvalues obtained from the bifurcation analysis. The system behaviour might be stable or unstable depending on these Eigenvalues. The different Eigenvalues obtained using the bifurcation method for the different cases with elevator deflection is shown in Table 5. The Eigenvalues obtained are used to plot the pole-zero map to understand the system stability behaviour and are shown in the Figure 17. The poles plot for the different simulation cases shows that Cases 1, 2, 3, 5 and 6 show stable behaviour and Cases 4 and 7 show unstable behaviour due to the poles lying on the RHS plane.

Table 5. Eigen values obtained using Bifurcation method

16a.png

(a)

16b.png

(b)

16c.png

(c)

16d.png

(d)

Figure 16. Bifurcation diagram: (a) theta, q, pitch (combined); (b) θ in degree; (c) q in degree; (d) pitch vs dele deflection in degree

17.png

Figure 17. System poles plot for the simulation cases