Analysis of the Short-Term Dynamics of Morphing Aircraft Caused by Shape Change Based on the Open-Loop Response and the Reachable Set Theory

Abstract

This work investigates the short-term dynamics caused by shape changes of morphing aircraft. We select a symmetric variable sweep morphing aircraft as the object of study and establish a six-degree-of-freedom multi-loop cascade model, and the coupling between derivative terms is eliminated by matrix transformation. Considering that the change in aerodynamic shape significantly affects the aerodynamic forces of the aircraft in a short period of time, and the variation in mass distribution generates additional aerodynamic forces and moments, we analyze the effects of these factors on the dynamic characteristics of the aircraft based on the open-loop response starting from the steady-state flight conditions. In addition, we analyze the improvement in maneuvering performance brought by morphing as an additional control input. We apply reachable set theory to multi-loop equations of motion and use the size of the reachable set to measure the maneuverability of aircraft. The results confirm that morphing can effectively improve the maneuverability of the aircraft.

1. Introduction

Morphing aircraft can change its shape to ensure optimal flight performance in different flight conditions or to provide control authority [1]. Thus, it has a broader flight envelope and better multi-mission capability than a fixed-geometry aircraft [2]. However, morphing also brings some barriers to overcome, such as additional weight, the difficulties of maintenance, and the stability challenges due to changes in aerodynamic coefficients and the variation in mass distribution during the morphing process [3].

Currently, researchers are primarily focused on the long-term benefits of morphing for aircraft, such as better long-endurance performance [4,5], broader flight envelope [6], and greater multi-mission capability [7], but they may be failing to pay sufficient attention to the short-term characteristics brought about by morphing. Note that morphing brings changes in aerodynamic forces, additional inertial forces and moments, and a change in gravity moments during the transition process. These effects may have significant impacts on aircraft dynamics. In [8,9,10], open-loop simulations of the morphing aircraft were performed to verify the effect of the wing transition process, and the results showed that the speed, height, and pitch angle of the morphing aircraft are greatly changed due to morphing. The effect of inertial forces and moments during the morphing of the longitudinal model of a morphing aircraft was investigated in [11], and the results showed that the inertial forces and moments that include the rate terms have very little effect on the dynamics of the morphing aircraft. In [12], the stability conditions for the morphing process of the longitudinal morphing aircraft model are derived by the parameterization of morphing and mathematical derivation.

Additionally, it is essential to recognize the potential for morphing mechanisms to enhance an aircraft’s maneuverability beyond the capabilities of conventional control surfaces. For example, in [13], theoretical analysis and trim calculations showed that morphing can effectively improve an aircraft’s pitching maneuverability. In [14], control allocation-based asymmetric wing telescoping was used to enhance the maneuverability of the morphing aircraft with the telescopic wing. In [15], a variable-span and cambered-span morphing UAV was given, and the analysis results showed that both stability and maneuverability were improved by shape change. However, there no study exists that can visually measure the improvement of maneuverability that morphing brings to the aircraft.

In this paper, we establish a six-degree-of-freedom multi-loop cascade morphing aircraft model and eliminate the coupling between their derivative terms by matrix transformation. Further, we analyz the effects of morphing on the dynamics based on the open-loop response of the morphing process starting from steady-state flight conditions. Furthermore, we investigate the enhancement of aircraft maneuverability by morphing based on reachable set theory. The primary contributions and novelties of this paper are summarized in the following:

• This paper provides detailed steps on converted the coupled morphing aircraft dynamics equations into a form that can be calculated using a numerical solver;
• Unlike the longitudinal model in the existing works, this paper investigates the dynamics of the morphing process of the morphing aircraft based on a six-degree-of-freedom model and compares the forces and moments of inertial, aerodynamic, and gravity during this process;
• We use the size of the backward reachable set to measure the maneuverability and response speed of the morphing aircraft and present a clear visualization of the improvement of the aircraft’s maneuverability due to morphing.

The rest of this paper is organized as follows. The model of the morphing aircraft is given in Section 2. The analysis of morphing dynamics based on open-loop response characteristics is described in Section 3. The maneuvering performance analysis of morphing aircraft based on reachable set theory is presented in Section 4, and the conclusion of this paper is given in Section 5.

2. Model of Morphing Aircraft

The Firebee UAV is selected as the baseline aircraft in this study [16], and we assume that it can morph according to the morphing concept proposed by NextGen Aeronautics [17]. As shown in Figure 1, the shape and area of its wings vary significantly with the sweep angle $\zeta$ ranging between 16 and 55 degrees.

The choice of this particular morphing methodology is based on its ability to significantly alter the wing area and mass distribution compared to variable wingspan and variable camber morphing methodologies. As such, it is a more representative morphing methodology for analyzing the short-term dynamics of a morphing aircraft.

2.1. Multi-Loop Cascaded Dynamic Equations

The equations of motion for the morphing aircraft are established based on the morphing dynamics in [16] and the cascade form of dynamic equations in [18], as shown below. The equations involve transformations between the coordinate systems of the flight path($F_H$), wind-axes($F_V$), body-fixed($F_b$), and earth-fixed($F_I$), which can be found in [4].

Based on Newton’s second law, the equations of motion in the flight path loop are given as

$$\left[\begin{array}{c}m\dot{V}\\ mVcos\gamma \dot{\chi }\\ -V\dot{\gamma }\end{array}\right]={\mathit{C}}_{\mathit{HB}}\left(\left[\begin{array}{c}{f}_{ix}\\ f_{iy}\\ f_{iz}\end{array}\right]+\left[\begin{array}{c}T\\ 0\\ 0\end{array}\right]\right)+{{\mathit{C}}_{\mathit{HV}}}^{}\left[\begin{array}{c}-D\\ 0\\ -L\end{array}\right]+{\mathit{C}}_{\mathit{HI}}\left[\begin{array}{c}0\\ 0\\ mg\end{array}\right].$$

(1)

where m is the mass of the aircraft, V is the velocity of the aircraft, $\chi$ is the kinematic azimuth angle, $\gamma$ is the flight path angle, T is the engine thrust, ${\mathit{C}}_{\mathit{HB}}$, ${\mathit{C}}_{\mathit{HV}}$, and ${\mathit{C}}_{\mathit{HI}}$ are the rotation matrices [19] between F_H, F_b, F_V and F_I, and their expressions are given as

$${\mathit{C}}_{\mathit{HB}}=\left[\begin{array}{ccc}cos\alpha cos\beta & sin\beta & sin\alpha cos\beta \\ -cos\alpha sin\beta cos\mu +sin\alpha sin\mu & cos\beta cos\mu & -sin\alpha sin\beta cos\mu -cos\alpha sin\mu \\ cos\alpha sin\beta sin\mu -sin\alpha cos\mu & cos\beta sin\mu & -sin\alpha sin\beta sin\mu +cos\alpha cos\mu \end{array}\right],$$

$${\mathit{C}}_{\mathit{HV}}=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\mu & -sin\mu \\ 0& sin\mu & cos\mu \end{array}\right],$$

and

$${\mathit{C}}_{\mathit{HI}}=\left[\begin{array}{ccc}cos\chi cos\gamma & -sin\chi & cos\chi sin\gamma \\ sin\chi cos\gamma & cos\chi & sin\chi sin\gamma \\ -sin\gamma & 0& cos\gamma \end{array}\right].$$

where $\mu$ is the bank angle, $\alpha$ is the angle of attack, and $\beta$ is the angle of the sideslip. Further, $D$ and $L$ denote drag and lift, respectively; they can be formulated as

$$\begin{array}{c}D=\overline{q}{S}_w\left(\zeta \right)\left(C_D^{{\alpha }^2}\left(\zeta ,V,h\right){\alpha }^2+C_D^{\alpha }\left(\zeta ,V,h\right)\alpha +C_D^0\left(\zeta ,V,h\right)\right),\hfill \\ L=\overline{q}{S}_w\left(\zeta \right)\left(C_L^{\alpha }\left(\zeta ,V,h\right)\alpha +C_L^0\left(\zeta ,V,h\right)\right).\hfill \end{array}$$

(2)

where $\zeta$ is the sweep angle, $\overline{q}$ is the dynamic pressure, and $S_w\left(\zeta \right)$ is the wing area. $C_D^{{\alpha }^2}$, $C_D^{\alpha }$, $C_D^0$, $C_L^{\alpha }$, and $C_L^0$ are aerodynamic coefficients, and they are all functions of the altitude $h$, velocity $V$, and sweep angle $\zeta$. Furthermore, $f_{\mathit{Ix}}$, $f_{\mathit{Iy}}$, and $f_{\mathit{Iz}}$ are the components of the inertial force in $F_b$, and their expressions are obtained using the Kane method [16], which are presented below.

$$\left[\begin{array}{c}{f}_{\mathit{Ix}}\\ f_{\mathit{Iy}}\\ f_{\mathit{Iz}}\end{array}\right]=\left[\begin{array}{c}\dot{q}{S}_z-\dot{r}{S}_y+2(q{\dot{S}}_z-r{\dot{S}}_y)+q\left(p{S}_y-q{S}_x\right)-r\left(r{S}_x-p{S}_z\right)+{\ddot{S}}_x\\ \dot{r}{S}_x-\dot{p}{S}_z+2(r{\dot{S}}_x-p{\dot{S}}_z)+r\left(q{S}_z-r{S}_y\right)-p\left(p{S}_y-q{S}_x\right)+{\ddot{S}}_x\\ \dot{q}{S}_y-\dot{r}{S}_x+2(q{\dot{S}}_y-r{\dot{S}}_x)+p\left(r{S}_x-p{S}_z\right)-q\left(q{S}_z-p{S}_y\right)+{\ddot{S}}_z\end{array}\right]m.$$

(3)

where ${\left[\begin{array}{ccc}p& q& r\end{array}\right]}^{\mathrm{T}}$ is the angular rate, ${\left[\begin{array}{ccc}{S}_x& S_y& S_z\end{array}\right]}^{\mathrm{T}}$ represents the center of gravity of the plane, and its position is also a function of $\zeta$. For the sake of clarity and conciseness, we will simply use $S_x$, $S_y$, and $S_z$ instead of $S_x\left(\zeta \right)$, $S_y\left(\zeta \right)$, and $S_z\left(\zeta \right)$ in our writing.

Based on the rotation relationship between the reference coordinate systems, we can establish the equations of motion as

$$\left[\begin{array}{c}\dot{\mu }\\ \dot{\alpha }\\ \dot{\beta }\end{array}\right]={\left[\begin{array}{ccc}cos\alpha cos\beta & 0& sin\alpha \\ sin\beta & 1& 0\\ sin\alpha cos\beta & 0& -cos\alpha \end{array}\right]}^{-1}\left({\mathit{C}}_{\mathit{HB}}^{\mathrm{T}}\left[\begin{array}{c}-\dot{\chi }sin\gamma \\ \dot{\gamma }\\ \dot{\chi }cos\gamma \end{array}\right]+\left[\begin{array}{c}p\\ q\\ r\end{array}\right]\right).$$

(4)

According to the momentum theorem, the equations of motion in the angular rate loop are obtained as given below.

$$\mathit{J}\left[\begin{array}{c}\dot{p}\\ \dot{q}\\ \dot{r}\end{array}\right]=\left(\left[\begin{array}{c}{m}_{\mathit{Ax}}+m_{\mathrm{Gx}}+m_{\mathit{Ix}}\\ m_{\mathit{Ay}}+m_{\mathit{Gy}}+m_{\mathit{Iy}}\\ m_{\mathit{Az}}+m_{\mathit{Gz}}+m_{\mathit{Iz}}\end{array}\right]-\left[\begin{array}{c}p\\ q\\ r\end{array}\right]\times \mathit{J}\left[\begin{array}{c}p\\ q\\ r\end{array}\right]\right).$$

(5)

where

$$\mathit{J}=\left[\begin{array}{ccc}{J}_x& 0& -J_{\mathit{xz}}\\ 0& J_y& 0\\ -J_{\mathit{xz}}& 0& J_z\end{array}\right]$$

(6)

is the inertia matrix of the aircraft. The elements of the inertia matrix $\mathit{J}$ are also functions of $\zeta$. To simplify the writing, we use letter symbols to represent variables directly, omitting their functional forms. Further,

$$\begin{array}{c}{m}_{\mathit{Ax}}=C_{l\beta }\left(\zeta ,V,h\right)\beta +C_{\mathit{lp}}\left(\zeta ,V,h\right)\frac{\mathit{pb}\left(\zeta \right)}{2V}+C_{\mathit{lr}}\left(\zeta ,V,h\right)\frac{\mathit{rb}\left(\zeta \right)}{2V}\hfill \\ C_{l{\delta }_a}\left(\zeta ,V,h\right){\delta }_a+C_{l{\delta }_r}\left(\zeta ,V,h\right){\delta }_r\hfill \\ m_{\mathit{Ay}}=C_{m0}\left(\zeta ,V,h\right)+C_{m\alpha }\left(\zeta ,V,h\right)\alpha +C_{\mathit{mq}}\left(\zeta ,V,h\right)\frac{q\overline{C}\left(\zeta \right)}{2V}+C_{m{\delta }_e}\left(\zeta ,V,h\right){\delta }_e\hfill \\ m_{\mathit{Az}}=C_{n\beta }\left(\zeta ,V,h\right)\beta +C_{\mathit{np}}\left(\zeta ,V,h\right)\frac{\mathit{pb}\left(\zeta \right)}{2V}+C_{\mathit{nr}}\left(\zeta ,V,h\right)\frac{\mathit{rb}\left(\zeta \right)}{2V}\hfill \\ C_{n{\delta }_a}\left(\zeta ,V,h\right){\delta }_a+C_{n{\delta }_r}\left(\zeta ,V,h\right){\delta }_r\hfill \end{array}$$

(7)

is the array of the components of aerodynamic moments in $F_b$, where $C_{l\beta }$, $C_{\mathit{lp}}$, $C_{\mathit{lr}}$, $C_{l{\delta }_a}$, $C_{l{\delta }_r}$, $C_{m0}$, $C_{m\alpha }$, $C_{\mathit{mq}}$, $C_{m{\delta }_e}$, $C_{n\beta }$, $C_{\mathit{np}}$, $C_{\mathit{nr}}$, $C_{n{\delta }_a}$, and $C_{n{\delta }_r}$ are all the aerodynamic moment coefficients, and they are all functions of $V$, $h$, and $\zeta$. $b\left(\zeta \right)$ and $\overline{c}\left(\zeta \right)$ are the mean aerodynamic chord and wingspan, respectively, which are functions of $\zeta$. Furthermore,

$$\left[\begin{array}{c}{M}_{\mathit{Gx}}\\ M_{\mathit{Gy}}\\ M_{\mathit{Gz}}\end{array}\right]=\left[\begin{array}{c}{S}_x\\ S_y\\ S_z\end{array}\right]\times \left({\mathit{C}}_{B/I}\left[\begin{array}{c}0\\ 0\\ \mathit{mg}\end{array}\right]\right)$$

(8)

is the array of the components of moment caused by gravity in $F_b$. Furthermore, ${\left[\begin{array}{ccc}{m}_{\mathit{Ix}}& m_{\mathit{Iy}}& m_{\mathit{Iz}}\end{array}\right]}^{\mathrm{T}}$ is the array of components of moment caused by the inertia force in $F_b$. According to the results in [14], its expression is given below.

$$\begin{array}{c}{m}_{\mathit{Ix}}={\dot{J}}_{x}p-{\dot{J}}_{\mathit{xy}}q-{\dot{J}}_{\mathit{xz}}r+m{S}_y\left(\mathit{pv}-\mathit{qu}\right)-m{S}_z\left(\mathit{ru}-\mathit{pw}\right)+m\left(S_y{\ddot{S}}_z-S_z{\ddot{S}}_y\right)\hfill \\ +\mathit{mq}\left(S_x{\dot{S}}_y-S_y{\dot{S}}_x\right)-\mathit{mr}(S_z{\dot{S}}_x-S_x{\dot{S}}_z)+m{S}_y\dot{w}-m{S}_z\dot{v},\hfill \\ m_{\mathit{Iy}}=-{\dot{J}}_{\mathit{xy}}p+{\dot{J}}_{y}q-{\dot{J}}_{\mathit{yz}}r+m{S}_z\left(\mathit{qw}-\mathit{rv}\right)-m{S}_x\left(\mathit{pv}-\mathit{qu}\right)+m\left(S_z{\ddot{S}}_x-S_x{\ddot{S}}_z\right)\hfill \\ +\mathit{mr}\left(S_y{\dot{S}}_z-S_z{\dot{S}}_y\right)-\mathit{mp}(S_x{\dot{S}}_y-S_y{\dot{S}}_x)+m{S}_z\dot{u}-m{S}_x\dot{w},\hfill \\ m_{\mathit{Iz}}=-{\dot{J}}_{\mathit{xz}}p+{\dot{J}}_{\mathit{yz}}q-{\dot{J}}_{z}r+m{S}_x\left(\mathit{ru}-\mathit{pw}\right)-m{S}_y\left(\mathit{qw}-\mathit{rv}\right)+m\left(S_x{\ddot{S}}_y-S_y{\ddot{S}}_x\right)\hfill \\ +\mathit{mp}\left(S_z{\dot{S}}_x-S_x{\dot{S}}_z\right)-\mathit{mq}(S_y{\dot{S}}_z-S_z{\dot{S}}_y)+m{S}_x\dot{v}s.-m{S}_y\dot{u}.\hfill \end{array}$$

(9)

where $u$, $v$, and $w$ is the components of the array of the velocity vector in the body-fixed coordinate system, and their expressions are

$$\begin{array}{c}u=Vcos\alpha cos\beta ,\hfill \\ v=Vsin\beta ,\hfill \\ w=Vsin\alpha cos\beta .\hfill \end{array}$$

(10)

Among these three loops, the angular rate loop is the innermost loop with the fastest response, the attitude loop is the second, and the flight path loop has the slowest response. The three loops together form the dynamic equations of the morphing aircraft.

Morphing alters the shape of the aircraft, which subsequently changes the aerodynamic coefficients. This causes variations in the aerodynamic forces and moments that the aircraft experiences. Furthermore, morphing alters the mass distribution of the aircraft, leading to the emergence of additional inertial and gravitational forces and moments. All the differences between the dynamics of the morphing aircraft and those of the fixed-geometry aircraft can be traced to differences in the flight path and angular rate loops. It is important to note that the attitude loops of the morphing aircraft and the fixed-geometry aircraft have the same dynamics.

2.2. Model Decoupling

Unlike the horizontal and vertical decoupling that separate the equations of motion, the model decoupling in this context refers to the algebraic separation of the variables, which aims to make the equation solvable by numerical algorithms. Note that commonly used numerical tools for solving ordinary differential equations require the form $\dot{\mathit{y}}=\mathit{f}\left(\mathit{y}\right)$ to calculate $\mathit{f}\left(\mathit{y}\right)$ at discrete points. However, in the flight path loop, $\dot{p}$, $\dot{q}$, and $\dot{r}$ appear on the right-hand side of the dynamic equations, and in the angular rate loop, $\dot{u}$, $\dot{v}$, and $\dot{w}$ appear on the right-hand side of the equations of motion. Because of coupling, numerical calculations are not possible, and we need to decouple the equations before we can simulate or analyze the aircraft’s dynamics.

Remark 1.

Numerical methods used for solving systems of ordinary differential equations involve the discretization of equations and iterative techniques to generate solutions. These solutions are approximated at discrete time intervals using numerical methods. The selection of suitable numerical methods with convergence properties and appropriate time-step sizes holds the key to achieving accurate and stable solutions, thereby preventing error accumulation. For our study, we have employed the ode4 solver with a fixed timestep of 0.02 s.

According to Equation (10)

$$\begin{array}{c}\dot{u}=\dot{V}cos\alpha cos\beta -Vsin\alpha cos\beta \dot{\alpha }-Vcos\alpha sin\beta \dot{\beta },\hfill \\ \dot{v}s.=\dot{V}sin\beta +Vcos\beta \dot{\beta },\hfill \\ \dot{w}=Vsin\alpha cos\beta +Vcos\alpha cos\beta \dot{\alpha }-Vsin\alpha sin\beta \dot{\beta }.\hfill \end{array}$$

(11)

Moving the terms in Equation (1) and (5) that contain $\dot{p}$, $\dot{q}$, $\dot{r}$, $\dot{V}$, $\dot{\alpha }$, and $\dot{\beta }$ to the left-hand side, we obtain

$$\left[\begin{array}{cc}{\mathit{I}}_{\mathbf{6}\mathbf{\times }\mathbf{6}}& {\mathit{M}}_{\mathbf{12}}\\ {\mathit{M}}_{\mathbf{21}}& \mathit{J}\end{array}\right]\left[\begin{array}{c}{\dot{\mathit{s}}}_{\mathbf{1}}\\ {\dot{\mathit{s}}}_{\mathbf{2}}\end{array}\right]=\left[\begin{array}{c}{\mathit{f}}_{\mathbf{1}}\\ {\mathit{f}}_{\mathbf{2}}\end{array}\right].$$

(12)

where ${\mathit{s}}_{\mathbf{1}}={\left[\begin{array}{cccccc}V& \chi & \gamma & \mu & \alpha & \beta \end{array}\right]}^{\mathrm{T}}$, $s_2={\left[\begin{array}{ccc}p& q& r\end{array}\right]}^{\mathrm{T}}$,

$${\mathit{M}}_{\mathbf{12}}=\left[\begin{array}{ccc}\frac{-1}{m}\left(-C_{12}{S}_z+C_{13}{S}_y\right)& \frac{-1}{m}\left(-C_{11}{S}_z+C_{13}{S}_x\right)& \frac{-1}{m}\left(-C_{11}{S}_y+C_{12}{S}_x\right)\\ \frac{-1}{\mathit{mV}c\gamma }\left(-C_{22}{S}_z+C_{23}{S}_y\right)& \frac{-1}{\mathit{mV}\mathrm{c}\gamma }\left(-C_{21}{S}_z+C_{23}{S}_x\right)& \frac{-1}{\mathit{mV}\mathrm{c}\gamma }\left(-C_{11}{S}_y+C_{22}{S}_x\right)\\ \frac{1}{\mathit{mV}}\left(-C_{32}{S}_z+C_{13}{S}_y\right)& \frac{1}{\mathit{mV}}\left(-C_{31}{S}_z+C_{33}{S}_x\right)& \frac{1}{\mathit{mV}}\left(-C_{31}{S}_y+C_{32}{S}_x\right),\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right],$$

where

$$\left[\begin{array}{ccc}{C}_{11}& C_{12}& C_{13}\\ C_{21}& C_{22}& C_{23}\\ C_{31}& C_{32}& C_{33}\end{array}\right]=\left[\begin{array}{ccc}cos\alpha cos\beta & sin\beta & sin\alpha cos\beta \\ \begin{array}{c}\hfill -cos\alpha sin\beta cos\mu \hfill \\ \hfill +sin\alpha sin\mu \hfill \end{array}& cos\beta cos\mu & \begin{array}{c}\hfill -sin\alpha sin\beta cos\mu \hfill \\ \hfill -cos\alpha sin\mu \hfill \end{array}\\ \begin{array}{c}\hfill cos\alpha sin\beta sin\mu \hfill \\ \hfill -sin\alpha cos\mu \hfill \end{array}& cos\beta sin\mu & \begin{array}{c}\hfill -sin\alpha sin\beta sin\mu \hfill \\ \hfill +cos\alpha cos\mu \hfill \end{array}\end{array}\right]$$

and

$${\mathit{M}}_{\mathbf{21}}=\left[\begin{array}{ccc}{S}_z\mathrm{s}\beta -S_y\mathrm{s}\alpha \mathrm{c}\beta & -S_{y}V\mathrm{c}\alpha \mathrm{c}\beta & S_{z}V\mathrm{c}\beta +S_{y}V\mathrm{s}\alpha \mathrm{s}\beta \\ -S_z\mathrm{c}\alpha \mathrm{c}\beta +S_x\mathrm{s}\alpha \mathrm{c}\beta & S_{z}V\mathrm{s}\alpha \mathrm{c}\beta +S_{x}V\mathrm{c}\alpha \mathrm{c}\beta & S_{z}V\mathrm{c}\alpha \mathrm{s}\beta -S_{x}V\mathrm{s}\alpha \mathrm{s}\beta \\ -S_x\mathrm{s}\beta +S_y\mathrm{c}\alpha \mathrm{c}\beta & -S_{y}V\mathrm{s}\alpha \mathrm{c}\beta & -S_{y}V\mathrm{c}\alpha \mathrm{s}\beta -S_{x}V\mathrm{c}\beta \end{array}\right].$$

where $s\alpha$, $c\alpha$, $s\beta$, and $c\beta$ denote $sin\alpha$, $cos\alpha$, $sin\beta$, and $cos\beta$, respectively. Further

$${\mathit{f}}_{\mathbf{1}}=\left[\begin{array}{c}{\mathit{f}}_{\mathbf{11}}\\ {\mathit{f}}_{\mathbf{12}}\end{array}\right]=\left[\begin{array}{c}{C}_{\mathit{HB}}^{}\left(\left[\begin{array}{c}{f}_{\mathit{ixr}}\\ f_{\mathit{iyr}}\\ f_{\mathit{izr}}\end{array}\right]+\left[\begin{array}{c}T\\ 0\\ 0\end{array}\right]\right)+C_{\mathit{HV}}^{}\left[\begin{array}{c}-D\\ 0\\ -L\end{array}\right]+C_{\mathit{HI}}\left[\begin{array}{c}0\\ 0\\ \mathit{mg}\end{array}\right]\\ {\left[\begin{array}{ccc}cos\alpha cos\beta & 0& sin\alpha \\ sin\beta & 1& 0\\ sin\alpha cos\beta & 0& -cos\alpha \end{array}\right]}^{-1}\left(C_{\mathit{HB}}^T\left[\begin{array}{c}-\dot{\chi }sin\gamma \\ \dot{\gamma }\\ \dot{\chi }cos\gamma \end{array}\right]+\left[\begin{array}{c}p\\ q\\ r\end{array}\right]\right)\end{array}\right],$$

(13)

$${\mathit{f}}_{\mathbf{2}}=\left[\begin{array}{c}{m}_{\mathit{Ax}}+m_{\mathrm{Gx}}+m_{\mathit{Ixr}}\\ m_{\mathit{Ay}}+m_{\mathit{Gy}}+m_{\mathit{Iyr}}\\ m_{\mathit{Az}}+m_{\mathit{Gz}}+m_{\mathit{Izr}}\end{array}\right]-\left[\begin{array}{c}p\\ q\\ r\end{array}\right]\times J\left[\begin{array}{c}p\\ q\\ r\end{array}\right].$$

(14)

where

$$\begin{array}{c}\hfill f_{\mathit{Ixr}}=f_{\mathit{Ix}}-\left(\dot{q}{S}_z-\dot{r}{S}_y\right),\hfill \\ \hfill f_{\mathit{Iyr}}=f_{\mathit{Iy}}-\left(\dot{r}{S}_x-\dot{p}{S}_z\right),\hfill \\ \hfill f_{\mathit{Izr}}=f_{\mathit{Iz}}-\left(\dot{p}{S}_y-\dot{q}{S}_x\right),\hfill \end{array}$$

(15)

and

$$\begin{array}{c}{m}_{\mathit{Ixr}}=m_{\mathit{Ix}}-m{S}_y\left(\dot{V}sin\alpha cos\beta +Vcos\alpha cos\beta \dot{\alpha }-Vsin\alpha sin\beta \dot{\beta }\right)+\hfill \\ +m{S}_z\left(\dot{V}sin\beta +Vcos\beta \dot{\beta }\right)\hfill \\ m_{\mathit{Iyr}}=m_{\mathit{Iy}}-m{S}_z\left(\dot{V}cos\alpha cos\beta -Vsin\alpha cos\beta \dot{\alpha }-Vcos\alpha sin\beta \dot{\beta }\right)\hfill \\ +m{S}_x\left(\dot{V}sin\alpha cos\beta +Vcos\alpha cos\beta \dot{\alpha }-Vsin\alpha sin\beta \dot{\beta }\right),\hfill \\ m_{\mathit{Izr}}=m_{\mathit{Iz}}+m{S}_y\left(\dot{V}cos\alpha cos\beta -Vsin\alpha cos\beta \dot{\alpha }-Vcos\alpha sin\beta \dot{\beta }\right)\hfill \\ -m{S}_x\left(\dot{V}sin\beta +Vcos\beta \dot{\beta }\right).\hfill \end{array}$$

(16)

In this way, the equations of motion of the morphing aircraft become a form that can be easily solved using numerical methods. It should be noted that there is only an algebraic transformation in this subsection, thus the resulting Equation (13) is mathematically equivalent to the known results in [4,20].

Remark 2.

Note that Equation (13) includes $\dot{\chi }$ and $\dot{\gamma }$ on the right-hand side. These terms are not moved to the left side due to three reasons. First, they are not related to morphing and are also present in the equations of the fixed-geometry aircraft’s attitude loop dynamics. Second, obtaining $\dot{\chi }$ and $\dot{\gamma }$ using Equation (1) before solving Equation (12) can help avoid any interference during the solving process. Third, moving them to the left-hand side requires matrix inversion and multiplication, resulting in a burdensome algebraic workload.

2.3. Parameter Setting

The total weight of the Firebee UAV is 907 kg, the wing weight is 120 kg, the length of the fuselage is 6.98 m, and the center of gravity of the fuselage is 3.05 m from the nose. We set the incidence angle of the wing to 0 degrees and utilized the NASA0010 airfoil. The detailed geometric parameters including the wing root chord, wing tip chord, wing center of gravity position, and trailing edge sweep for various configurations can be obtained from previously published literature [9,16].

Further, the values of the inertia parameters including $S_x$, $S_y$, $S_z$, $J_x$, $J_y$, $J_z$, $J_{\mathit{xz}}$, $S_{\mathit{ref}}$, $b$, and $\overline{c}$ under different configurations can be found in [9,20]. We utilized the method of polynomial fitting to derive functions of these parameters with respect to $\zeta$. In this process, we assumed $S_y$ and $S_z$ to be zero and fitted $S_x$, $b$, and $\overline{c}$ as primary polynomial functions while fitting $J_x$, $J_y$, $J_z$, $J_{\mathit{xz}}$, and $S_{\mathit{ref}}$ as quadratic polynomial functions. Note that the advantage of adopting the method of polynomial fitting is that the derivatives of these inertial parameters can be derived analytically for a given morphing process.

Furthermore, the aerodynamic coefficients for different altitudes, Mach numbers, and sweep angles can be found in [21], and they form a three-dimensional lookup table, and the interpolation is used to construct the aerodynamic derivatives during the transition.

3. Analysis of Morphing Dynamics Based on Open-Loop Response Characteristics

This section addresses the case of the morphing aircraft undergoing a change in shape from steady-state flight. Specifically, we will examine two types of steady-state flight conditions: steady wing-level flight and steady turning flight.

3.1. Steady Wing-Level Flight

When considering steady wing-level flight conditions, all the derivative terms at the left-hand side of Equation (13) are equal to 0. It is crucial to note that there are nine equations and fifteen variables, which include $V$, $\chi$, $\gamma$, $\mu$, $\alpha$, $\beta$, $p$, $q$, $r$, $\zeta$, $T$, ${\delta }_a$, ${\delta }_e$, ${\delta }_r$, and $h$. As long as six of these variables are specified, the trim point can be computed. To illustrate the application of our model, we consider a typical cruising condition for an aircraft with a speed of $V=200$ m/s and altitude of $h=5000$ m. Combining the conditions of $p=q=r=0$ under the steady wing-level flight condition, the trim point is derived as

$$\begin{array}{c}\left[V^t,{\chi }^t,{\gamma }^t,{\mu }^t,{\alpha }^t,{\beta }^t,p^t,q^t,r^t\right]=\left[250\mathrm{m}/\mathrm{s},0^{\circ },0^{\circ },0^{\circ },0.504^{\circ },0^{\circ },0^{\circ },0^{\circ },0^{\circ }\right],\hfill \\ \left[T^t,{{\delta }_a}^t,{{\delta }_e}^t,{{\delta }_r}^t,{\zeta }^t,h^t\right]=\left[2349.2N,0^{\circ },-{1.42}^{\circ },0^{\circ },{35}^{\circ },5000\mathrm{m}\right].\hfill \end{array}$$

The control input is fixed as $T\equiv T^t,{\delta }_a\equiv {{\delta }_a}^t,{\delta }_e\equiv {{\delta }_e}^t,{\delta }_r\equiv {{\delta }_r}^t$ throughout the simulation, while the aircraft undergoes a morphing process as shown in Figure 2. The sweep angle changes as follows: it increases uniformly from an initial angle of ${35}^{\circ }$ to ${55}^{\circ }$ at a rate of ${10}^{\circ }$/s between 10 s and 12 s, remains constant until 17 s, and then decreases uniformly at a rate of ${10}^{\circ }$/s to the final angle of ${35}^{\circ }$ between 17 s and 19 s. The aerodynamic forces, inertial forces, aerodynamic moments, inertia moments, and gravity moments in the body-fixed coordinated system during the open-loop response are shown in Figure 3 and Figure 4. The results indicate that the aircraft’s dynamic characteristics are minimally affected by inertial forces, but significant effects result from both the inertial and gravitational moments caused by the displacement of the centroid.

Figure 5 shows the state response of the aircraft during the morphing process. It can be seen that the velocity of the aircraft initially increases and then slowly returns to the trim point due to drag reduction resulting from the configuration morph to $\zeta ={55}^{\circ }$. The flight path angle $\gamma$ initially decreases and then gradually returns to 0. These slow changes in $V$ and $\gamma$ reflect the long-period (phugoid) mode of the aircraft, and the ability of the aircraft to return to the trim point demonstrates its longitudinal dynamic stability. Further, the value of $\alpha$ changes slowly with the variation in $\gamma$ because $\alpha =\theta -\gamma$, and note that the response of $\alpha$ exhibits two peaks during the transitions of the sweep angle from ${35}^{\circ }$ to ${55}^{\circ }$ and from ${55}^{\circ }$ to ${35}^{\circ }$. This response is attributed to the rapid changes in $q$, causing $\alpha$ to follow suit based on Equation (4). Furthermore, as the aircraft is designed to be longitudinally statically stable, $\alpha$ returns to the trim point and generates fast oscillations. These oscillations reflect the excitation of the short-period mode of the aircraft. Additionally, $\chi$, $\mu$, $\beta$, $q$, and $r$ remain constant at 0 because there are no lateral forces acting on the aircraft.

To quantify the effects of inertia forces, moments of inertia and gravitational moments due to morphing, we define $r{f}_{\mathit{Ii}}$, $r{m}_{\mathit{Ii}}$, and $r{m}_{\mathit{Gi}}$ as follows

$$\begin{array}{c}r{f}_{\mathit{Ii}}=\underset{t}{max}\left|f_{\mathit{Ii}}\right|/\underset{t}{\mathit{means}}\left(\left|f_{\mathit{Ai}}\right|\right),i=x,y,z,\hfill \\ r{m}_{\mathit{Ii}}=\underset{t}{max}\left|m_{\mathit{Ii}}\right|/\underset{t}{\mathit{means}}\left(\left|m_{\mathit{Ai}}\right|\right),i=x,y,z,\hfill \\ r{m}_{\mathit{Gi}}=\underset{t}{max}\left|m_{\mathrm{Gi}}\right|/\underset{t}{\mathit{means}}\left(\left|m_{\mathit{Gi}}\right|\right),i=x,y,z.\hfill \end{array}$$

In open-loop simulations, these values indicate how the maximum inertial force, inertia moment, and gravitational moment of the aircraft compare to the average aerodynamic force and moment over the course of the simulation. Specifically, they represent the ratio between these two quantities. We perform open-loop simulations in various flight states and record these parameters, which are presented in

Table 1. Effect of morphing in various flight conditions under steady wing-level flight (part 1).

Flight Condition	Lose Stability	$\mathit{r}{\mathit{f}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{f}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{f}}_{\mathit{Iz}}$
$h=2000$ m, $V=150$ m/s	No	$5.23\times {10}^{-5}$	0	$3.74\times {10}^{-4}$
$h=2000$ m, $V=200$ m/s	No	$1.49\times {10}^{-5}$	0	$2.57\times {10}^{-4}$
$h=2000$ m, $V=250$ m/s	No	$5.35\times {10}^{-6}$	0	$4.64\times {10}^{-4}$
$h=5000$ m, $V=150$ m/s	No	$7.89\times {10}^{-5}$	0	$3.32\times {10}^{-4}$
$h=5000$ m, $V=200$ m/s	No	$1.85\times {10}^{-5}$	0	$2.13\times {10}^{-4}$
$h=5000$ m, $V=250$ m/s	No	$5.99\times {10}^{-6}$	0	$1.47\times {10}^{-4}$
$h=8000$ m, $V=150$ m/s	No	$8.82\times {10}^{-5}$	0	$3.17\times {10}^{-4}$
$h=8000$ m, $V=200$ m/s	No	$1.83\times {10}^{-5}$	0	$1.82\times {10}^{-4}$
$h=8000$ m, $V=250$ m/s	No	$4.64\times {10}^{-6}$	0	$1.16\times {10}^{-4}$

and

Table 2. Effect of morphing in various flight conditions under steady wing-level flight (part 2).

Flight Condition	$\mathit{r}{\mathit{m}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iz}}$	$\mathit{r}{\mathit{m}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iz}}$
$h=2000$ m, $V=150$ m/s	0	$0.2469$	0	0	$1.19$	0
$h=2000$ m, $V=200$ m/s	0	$0.2338$	0	0	$1.10$	0
$h=2000$ m, $V=250$ m/s	0	$0.2180$	0	0	$1.12$	0
$h=5000$ m, $V=150$ m/s	0	$0.2382$	0	0	$1.18$	0
$h=5000$ m, $V=200$ m/s	0	$0.2201$	0	0	$1.09$	0
$h=5000$ m, $V=250$ m/s	0	$0.2006$	0	0	$1.030$	0
$h=8000$ m, $V=150$ m/s	0	$0.2281$	0	0	$1.17$	0
$h=8000$ m, $V=200$ m/s	0	$0.2004$	0	0	$1.08$	0
$h=8000$ m, $V=250$ m/s	0	$0.1686$	0	0	$1.01$	0

. Additionally, we investigate different morphing processes by designing two other configurations with initial angles of $\zeta ={16}^{\circ }$ and $\zeta ={55}^{\circ }$, as shown in Figure 6 and Figure 7, respectively. We perform open-loop simulations of these morphing processes under level flight conditions at an altitude of $h=5000$ m and a velocity of $V=250$ m/s. The obtained results for $r{f}_{\mathit{Ii}}$, $r{M}_{\mathit{Ii}}$, and $r{M}_{\mathit{Gi}}$ are presented in

Table 3. Effect of morphing in different morphing processes under steady wing-level flight (part 1).

Initial Configuration	Lose Stability	$\mathit{r}{\mathit{f}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{f}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{f}}_{\mathit{Iz}}$
$\zeta ={16}^{\circ }$	No	$4.68\times {10}^{-6}$	0	$1.61\times {10}^{-4}$
$\zeta ={35}^{\circ }$	No	$1.49\times {10}^{-5}$	0	$2.57\times {10}^{-4}$
$\zeta ={55}^{\circ }$	No	$1.35\times {10}^{-5}$	0	$1.40\times {10}^{-4}$

and

Table 4. Effect of morphing in different morphing processes under steady wing-level flight (part 2).

Initial Configuration	$\mathit{r}{\mathit{m}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iz}}$	$\mathit{r}{\mathit{m}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iz}}$
$\zeta ={16}^{\circ }$	0	$0.225$	0	0	$1.16$	0
$\zeta ={35}^{\circ }$	0	$0.234$	0	0	$1.10$	0
$\zeta ={55}^{\circ }$	0	$0.425$	0	0	$2.08$	0

.

The simulation results indicate that the impact of inertial forces is relatively small during open-loop simulations under different flight conditions. However, the effects of inertial and gravitational moments are more substantial, and they have a more significant influence on the dynamic behavior of the aircraft.

3.2. Steady Turning Flight

During stable turning flight, only the derivative term $\dot{\chi }$ on the left-hand side of Equation (13) is constant, while the other terms are set to 0. Similarly, specifying six variables allows for the calculation of the trim point. Let us consider an aircraft with an airspeed of $V=150$ m/s, altitude $h=5000$ m, and a turning rate of $\dot{\chi }=5^{\circ }$/s. Given a bank angle of $\zeta ={35}^{\circ }$ and a yaw angle of $\chi =0^{\circ }$, along with a pitch angle of $\gamma =0^{\circ }$ and a sideslip angle of $\beta =0^{\circ }$ under the steady turning flight condition, we can derive the trim point as

$$\begin{array}{c}\left[V^t,{\chi }^t,{\gamma }^t,{\mu }^t,{\alpha }^t,{\beta }^t,p^t,q^t,r^t\right]=\left[150\mathrm{m}/\mathrm{s},0^{\circ },0^{\circ },{53.19}^{\circ },4.0068^{\circ },0^{\circ },\right.\hfill \\ \left.-{0.209}^{\circ }/\mathrm{s},4.0047^{\circ }/s,{2.99}^{\circ }/\mathrm{s}\right],\hfill \\ \left[T^t,{{\delta }_a}^t,{{\delta }_e}^t,{{\delta }_r}^t,{\zeta }^t,h^t\right]=\left[1683.2N,-0.9074^{\circ },-6.486^{\circ },{5.444}^{\circ },{35}^{\circ },5000\mathrm{m}\right].\hfill \end{array}$$

The morphing process still starts at 10 s, and the sweep angle of the aircraft increases uniformly from ${35}^{\circ }$ to ${55}^{\circ }$ at 12 s, and then maintains $\zeta ={55}^{\circ }$. The open loop responses of the morphing aircraft dynamic are shown in Figure 8, Figure 9 and Figure 10. As can be seen, after the onset of the morphing, the aircraft loses stability within 3 s if the control input is kept fixed. This is mainly due to the additional moments generated by morphing while the inertial force remains small.

To obtain the values of $r{f}_{\mathit{Ii}}$, $r{M}_{\mathit{Ii}}$, and $r{M}_{\mathit{Gi}}$ under various steady turning conditions, we can follow the same procedures outlined previously. The results are shown in

Table 5. Effect of morphing in various flight conditions under steady turning flight (part 1).

Flight Condition	Lose Stability	$\mathit{r}{\mathit{f}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{f}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{f}}_{\mathit{Iz}}$
$h=2000$ m, $V=150$ m/s	Yes	$0.0264$	$0.0050$	$4.71\times {10}^{-4}$
$h=2000$ m, $V=200$ m/s	Yes	$0.0100$	$0.0050$	$4.83\times {10}^{-4}$
$h=2000$ m, $V=250$ m/s	Yes	$0.0091$	$0.0045$	$5.17\times {10}^{-4}$
$h=5000$ m, $V=150$ m/s	Yes	$0.0275$	$0.0050$	$5.50\times {10}^{-4}$
$h=5000$ m, $V=200$ m/s	Yes	$0.0206$	$0.0041$	$6.54\times {10}^{-4}$
$h=5000$ m, $V=250$ m/s	Yes	$0.0155$	$0.0033$	$7.48\times {10}^{-4}$
$h=8000$ m, $V=150$ m/s	Yes	$0.0280$	$0.0050$	$5.77\times {10}^{-4}$
$h=8000$ m, $V=200$ m/s	Yes	$0.0214$	$0.0040$	$7.18\times {10}^{-4}$
$h=8000$ m, $V=250$ m/s	Yes	$0.0164$	$0.0033$	$8.66\times {10}^{-4}$

and

Table 6. Effect of morphing in various flight conditions under steady turning flight (part 2).

Flight Condition	$\mathit{r}{\mathit{m}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iz}}$	$\mathit{r}{\mathit{m}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iz}}$
$h=2000$ m, $V=150$ m/s	$0.3371$	$5.017$	$0.9271$	0	$1.0853$	$1.0074$
$h=2000$ m, $V=200$ m/s	$0.6695$	$6.015$	$0.8834$	0	$0.6532$	$0.8781$
$h=2000$ m, $V=250$ m/s	$0.5682$	$5.364$	$0.8833$	0	$0.3309$	$0.8741$
$h=5000$ m, $V=150$ m/s	$0.4419$	$4.877$	$0.9295$	0	$1.2685$	$0.9980$
$h=5000$ m, $V=200$ m/s	$0.3245$	$5.705$	$0.9182$	0	$0.5963$	$0.9918$
$h=5000$ m, $V=250$ m/s	$0.2699$	$5.743$	$0.9078$	0	$0.2812$	$0.9979$
$h=8000$ m, $V=150$ m/s	$0.5476$	$4.808$	$0.9294$	0	$1.4168$	$0.9924$
$h=8000$ m, $V=200$ m/s	$0.4043$	$5.575$	$0.9185$	0	$0.6727$	$0.9890$
$h=8000$ m, $V=250$ m/s	$0.3051$	$5.498$	$0.9082$	0	$0.3165$	$0.9776$

. Moreover, we also conducted simulations for initial configurations with $\zeta ={16}^{\circ }$ and $\zeta ={55}^{\circ }$, as depicted in Figure 6 and Figure 7, respectively.

Table 7. Effect of morphing in different morphing processes under steady turning flight (part 1).

Initial Configuration	Lose Stability	$\mathit{r}{\mathit{f}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{f}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{f}}_{\mathit{Iz}}$
$\zeta ={16}^{\circ }$	No	$5.23\times {10}^{-5}$	0	$3.74\times {10}^{-4}$
$\zeta ={35}^{\circ }$	No	$1.49\times {10}^{-5}$	0	$2.57\times {10}^{-4}$
$\zeta ={55}^{\circ }$	No	$5.35\times {10}^{-6}$	0	$4.64\times {10}^{-4}$

and

Table 8. Effect of morphing in different morphing processes under steady turning flight (part 2).

Initial Configuration	$\mathit{r}{\mathit{m}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iz}}$	$\mathit{r}{\mathit{m}}_{\mathit{Ix}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iy}}$	$\mathit{r}{\mathit{m}}_{\mathit{Iz}}$
$\zeta ={16}^{\circ }$	0	$0.2469$	0	0	$1.19$	0
$\zeta ={35}^{\circ }$	0	$0.2338$	0	0	$1.10$	0
$\zeta ={55}^{\circ }$	0	$0.2180$	0	0	$1.12$	0

show the resulting $r{f}_{\mathit{Ii}}$, $r{M}_{\mathit{Ii}}$, and $r{M}_{\mathit{Gi}}$ values for each of the three morphing cases when the flight altitude is $h=5000$ m and the velocity is $V=250$ m/s.

The results show that similar to the steady wing-level flight condition, the effects of inertia moment and gravity moment are still dominant, whereas the effect of the inertial force is very small.

The open-loop response simulation results show that the morphing process has a considerable impact on the dynamical characteristics of the aircraft. In particular, the changes in aerodynamic forces and moments, inertia moments, and gravity moments due to morphing have significant effects. In contrast, the effects of inertia forces are generally negligible. Besides, the aircraft is more likely to lose stability due to morphing in the turning state than in the level flight state.

4. Maneuvering Performance Analysis of Morphing Aircraft Based on Reachable Set Theory

4.1. Preliminaries of Reachable Set Theory

Consider a continuous time system

$$\dot{x}=f\left(x,u\right),x\in R^n,u\left(·\right)\in U_{\left[0,T\right]}$$

where $x$ and $u$ are the state variable and control input of the system, respectively. $U_{\left[0,T\right]}$ denote the input constraint. If we assume that $f$ is Lipschitz continuous [22], then there exists a unique solution of the system, i.e., for any given initial state $x_0$, initial moment $t_0$, and control signal $u\left(·\right)\in U_{\left[0,T\right]}$, the solutions of the system are located on the unique solution trajectory, denoted by ${\xi }_f\left(t;x_0,t_0,u\left(·\right)\right)$.

Based on the above definitions and assumptions, the backward reachable set for the target set $\kappa$ is defined as given below.

$$R_b=\left\{x\in R^n\right.|\exists u\left(·\right)\in U_{\left[0,T\right]},\exists \tau \in \left[0,T\right],\left.{\xi }_f\left(\tau ;x_0,t_0,u\left(·\right)\right)\in \kappa \right\}.$$

It represents the set of initial states where each state point finds a control input that allows the system’s trajectory to enter the target set $\kappa$ in a specified time.

Remark 3.

A larger range of reachable set indicates the system’s ability to transfer a wider variety of initial states to the desired target set within a set timeframe. Therefore, the size of the backward reachable set is a measure of the system’s ability to move into a predetermined area. With respect to aircraft control, the ability to move through state space represents the aircraft’s maneuvering capabilities. Hence, in the following subsections, we analyze how morphing can enhance the aircraft’s maneuvering performance by studying the size of its backward reachable set.

The level-set method is a practical algorithm for solving the backward reachable set, and its workflow is given below.

Describe the target set as the form of a level set as

$$\kappa =\left\{x\in R^n,l\left(x\right)>0\right\},$$

(17)

and defining the level set function

$$V_1=\underset{u\in U_{\left[t,T\right]}}{inf}\underset{\tau \in \left[t,T\right]}{min}l\left({\xi }_f\left(\tau ;x,t,u\left(·\right)\right)\right).$$

(18)

Then, there exists the following relationship between $R_b$ and $V_1$ as

$$R_b\left(t,{\kappa }^c\right)=\left\{x\in R^n|V_1\left(x,t\right)<0\right\}.$$

(19)

It can be seen that the computation of the reachable set can be converted to the computation of $V_1$. According to the result in [22], the numeric solution of $V_1$ is the viscosity solution of the following Hamilton–Jacobi partial differential equation as

$$\begin{array}{c}\frac{\partial V_1}{\partial t}\left(x,t\right)+min\left(0,H_1^{*}\right)=0,\hfill \\ H_1^{*}=\underset{u\in U}{inf}\frac{\partial V_1}{\partial x}f\left(x\right),\hfill \\ V_1\left(x,T\right)=l\left(x\right).\hfill \end{array}$$

(20)

Therefore, the solution of the reachable set is transformed into the problem of solving the Hamilton–Jacobi equation, whose analytical solution is difficult to find, but the numerical solution can be approximated at discrete grid points with the grid-based finite difference method. By utilizing the Lax–Friedrich method [23] available in the MATLAB Level Set Toolbox [24], we can attain the numerical solution for the Hamilton–Jacobi partial differential equation presented in Equation (20).

4.2. Reachable Set of the Flight Path Loop

In this subsection, we will compare the size of the reachable sets in the flight path loop between the morphing aircraft and the fixed-geometry aircraft with $\zeta ={35}^{\circ }$. To accurately represent the morphing process, it is essential to account for the duration of the shape change. Therefore, we model the dynamic behavior of the morphing mechanism with the following equation:

$$\dot{\zeta }=10\left({\zeta }_{\mathit{des}}-\zeta \right).$$

(21)

Therefore, Equations (1) and (21) together form the dynamics equations of the flight path loop. Set $\beta =0^{\circ }$ and consider $\mu$, $\alpha$ and ${\zeta }_{\mathit{des}}$ as control inputs satisfying $\alpha \in [-{10}^{\circ },{15}^{\circ }]$, $\mu \in \left[-{60}^{\circ },{60}^{\circ }\right]$, ${\zeta }_{\mathit{des}}\in \left[{16}^{\circ },{55}^{\circ }\right]$. Further, based on the results in Section 3, we neglect the effect of the additional inertial forces. Then, the optimum value of its Hamiltonian function is given as

$$\begin{array}{c}{H}_{\mathit{fp}}^{*}=\left(\frac{-\rho S_{\mathit{ref}}{V}^2{C}_{D0}}{2m}-gsin\gamma +\frac{Tcos{\alpha }^{*}}{m}-\frac{\rho S_{\mathit{ref}}{V}^2}{2m}\left(C_{D\alpha }{\alpha }^{*}+C_{D{\alpha }^2}{\alpha }^{*2}\right)\right)p_1\hfill \\ +\left(\frac{\rho S_{\mathit{ref}}V}{2mcos\gamma }\left(C_{L0}+C_{L\alpha }{\alpha }^{*}\right)sin{\mu }^{*}+\frac{Tsin{\alpha }^{*}sin{\mu }^{*}}{\mathit{mV}cos\gamma }\right)p_2\hfill \\ -\left(\frac{g}{Vcos\gamma }+\frac{\rho S_{\mathit{ref}}V}{2m}\left(C_{L0}+C_{L\alpha }{\alpha }^{*}\right)cos{\mu }^{*}+\frac{Tsin{\alpha }^{*}cos{\mu }^{*}}{\mathit{mV}}\right)p_3\hfill \\ +10\left({\zeta }_{\mathit{des}}^{*}-\zeta \right)p_4.\hfill \end{array}$$

(22)

where $p_1=\frac{\partial V_{\mathit{fp}1}}{\partial V}$, $p_2=\frac{\partial V_{\mathit{fp}1}}{\partial \chi }$, $p_3=\frac{\partial V_{\mathit{fp}1}}{\partial \gamma }$, $V_{\mathit{fp}1}$ is the level-set function in the flight path loop.

$${\alpha }^{*}=\left\{\begin{array}{c}{\alpha }_{min},\mathrm{if}{p}_1>0,\widehat{\alpha }>\overline{\alpha }\mathrm{or}{p}_1=0,p_3>0,\mathrm{or}{p}_1<0,\widehat{\alpha }<{\alpha }_{min}\hfill \\ {\alpha }_{max},\mathrm{if}{p}_1>0,\widehat{\alpha }\le \overline{\alpha }\mathrm{or}{p}_1=0,p_3\le 0,\mathrm{or}{p}_1<0,\widehat{\alpha }\ge {\alpha }_{min}\hfill \\ \widehat{\alpha },\mathrm{if}{p}_1<0,{\alpha }_{min}<\widehat{\alpha }<{\alpha }_{max}\hfill \end{array}\right..$$

where $\overline{\alpha }=\frac{{\alpha }_{min}+{\alpha }_{max}}{2}$, $\widehat{\alpha }=\frac{1}{2p_{1}V{C}_{D{\alpha }^2}}\left(p_{3}cos\mu C_{L\alpha }-p_{1}V{C}_{D\alpha }+\frac{p_2{C}_{L\alpha }sin\mu }{cos\gamma }\right)$,

$${\mu }^{*}=\left\{\begin{array}{c}{\mu }_{min},\mathrm{if}{p}_2>0\hfill \\ {\mu }_{max},\mathrm{if}{p}_2<0\hfill \\ 0,\mathrm{if}{p}_2=0\hfill \end{array}\right.,$$

$${\zeta }_{\mathit{des}}^{*}=\left\{\begin{array}{c}{\zeta }_{min},\mathrm{if}{p}_4>0\hfill \\ {\zeta }_{max},\mathrm{if}{p}_4<0\hfill \\ \zeta ,\mathrm{if}{p}_4=0\hfill \end{array}\right..$$

Remark 4.

The optimal Hamiltonian function for a fixed-geometry aircraft is treated in a similar way, simply by removing Equation (21) from the dynamics equations. We omit this part for the sake of brevity.

We set the target set as

$$\begin{array}{c}{\kappa }_1=\left[V_{{\kappa }_{1}min}=140m/s,V_{{\kappa }_{1}max}=160m/s\right]\times \left[{\chi }_{{\kappa }_{1}min}=-5^{\circ },{\chi }_{{\kappa }_{1}max}=5^{\circ }\right]\hfill \\ \times \left[{\gamma }_{{\kappa }_{1}min}={30}^{\circ },{\gamma }_{{\kappa }_{1}max}={40}^{\circ }\right]\hfill \end{array}$$

and

$$\begin{array}{c}{\kappa }_2=\left[V_{{\kappa }_{2}min}=140m/s,V_{{\kappa }_{2}max}=160m/s\right]\times \left[{\chi }_{{\kappa }_{2}min}=-5^{\circ },{\chi }_{{\kappa }_{2}max}=5^{\circ }\right]\hfill \\ \times \left[{\gamma }_{{\kappa }_{2}min}=-{40}^{\circ },{\gamma }_{{\kappa }_{2}max}=-{30}^{\circ }\right]\hfill \end{array}$$

respectively for the reachable set computation, and their sizes reflect the ability of the aircraft to maneuver to the pull-up and dive states. We take the slices $\zeta ={35}^{\circ }$ of the morphing aircraft’s reachable sets to compare with the fixed-geometry aircraft. The results of the reachable set for $T=3$ s are shown in Figure 11 and Figure 12. It can be seen that the climbing reachable sets of the morphing aircraft roughly exhibit an elliptical cylinder shape, with $V$ range from 120 m/s to 180 m/s, $\chi$ range from −50${}^{\circ }$ to 50${}^{\circ }$, and $\gamma$ range from 20 ${}^{\circ }$ to 60 ${}^{\circ }$. Meanwhile, the diving reachable sets mainly focus on the region with velocities below the target set, with a smaller size compared to the climbing reachable sets. This is due to the inherent difficulty of reducing velocity during a dive. By counting the number of grid points in the reachable set, we estimate the size of the reachable set for the morphing aircraft in Figure 11 and Figure 12 to be 6.43% and 5.93% larger than that of the conventional configuration aircraft, respectively. This indicates that the morphing aircraft can slightly increase the response speed of the flight path loop by morphing, even if the difference in the response speed of the inner loop is not considered.

4.3. Reachable Set of the Angular Rate Loop

In the previous section, we discussed the response time and maneuverability of the flight path loop without considering the inner loop dynamics. In this section, we will analyze the dynamics of the inner angular rate loop. It is worth noting that the dynamics of the angular velocity loop are still dependent on the outer loop variables, such as $V$, $\mu$, $\alpha$, and $\beta$. According to the principle of time-scale separation, the slow outer loop variables can be treated as constants when analyzing the fast inner loop variables. Thus, we set the following parameters for the outer loop: $V=200\mathrm{m}/\mathrm{s}$, $\gamma =0^{\circ }$, $\alpha =0.5^{\circ }$, $\mu =0^{\circ }$, and $\beta =0^{\circ }$. The control inputs for the angular rate loop are ${\delta }_a$, ${\delta }_e$, ${\delta }_r$, and ${\zeta }_{\mathit{des}}$, with constraints set at ${\delta }_a\in \left[-{30}^{\circ },{30}^{\circ }\right]$, ${\delta }_e\in \left[-{30}^{\circ },{30}^{\circ }\right]$, ${\delta }_r\in \left[-{30}^{\circ },{30}^{\circ }\right]$, and ${\zeta }_{\mathit{des}}\in \left[-{16}^{\circ },{55}^{\circ }\right]$.

According to Equations (5) and (20), the optimal Hamiltonian function in the angular rate loop is given as

$$\begin{array}{c}{H}_a^{*}=\frac{1}{\Gamma }\left(J_{\mathit{xz}}\left(J_x-J_y+J_z\right)\mathit{pq}-\left(J_z\left(J_z-J_y\right)+J_{\mathit{xz}}^2\right)\mathit{qr}+J_z\left(l_a+l_g+l_i\right)+J_{\mathit{xz}}(n_a+n_g+n_i)\right)p_5\hfill \\ +\frac{1}{J_y}\left(\left(J_z-J_x\right)\mathrm{pr}-J_{\mathit{xz}}\left(p^2-r^2\right)+m_a+m_g+m_i\right)p_6\hfill \\ +\frac{1}{\Gamma }\left(\left(\left(J_x-J_y\right)J_x+J_{\mathit{xz}}^2\right)\mathit{pq}-J_{\mathit{xz}}\left(J_x-J_y+J_z\right)\mathit{qr}+J_{\mathit{xz}}\left(l_a+l_g+l_i\right)+J_x(n_a+n_g+n_i)\right)p_7\hfill \\ +10\left({\zeta }_{\mathit{des}}^{*}-\zeta \right)p_8.\hfill \end{array}$$

where $\Gamma =J_x{J}_z-{J_{x}z}^2$, $p_5=\frac{\partial V_{a1}}{\partial p}$, $p_6=\frac{\partial V_{a1}}{\partial q}$, $p_7=\frac{\partial V_{a1}}{\partial r}$, $p_8=\frac{\partial V_{a1}}{\partial \zeta }$, and $V_{a1}$ is the level set function in the angular rate loop, Further, the optimal control inputs is set as given below.

$${\delta }_a^{*}=\left\{\begin{array}{c}{\delta }_{amin},\mathrm{if}{p}_5>0\hfill \\ {\delta }_{amax},\mathrm{if}{p}_5<0\hfill \\ 0,\mathrm{if}{p}_5=0\hfill \end{array}\right.,$$

$${\delta }_e^{*}=\left\{\begin{array}{c}{\delta }_{emin},\mathrm{if}{p}_6>0\hfill \\ {\delta }_{emax},\mathrm{if}{p}_6<0\hfill \\ 0,\mathrm{if}{p}_6=0\hfill \end{array}\right.,$$

$${\delta }_r^{*}=\left\{\begin{array}{c}{\delta }_{rmin},\mathrm{if}{p}_6>0\hfill \\ {\delta }_{rmax},\mathrm{if}{p}_6<0\hfill \\ 0,\mathrm{if}{p}_6=0\hfill \end{array}\right.,$$

$${\zeta }_{\mathit{des}}^{*}=\left\{\begin{array}{c}{\zeta }_{min},\mathrm{if}{p}_8>0\hfill \\ {\zeta }_{max},\mathrm{if}{p}_8<0\hfill \\ 0,\mathrm{if}{p}_8=0\hfill \end{array}\right..$$

We set the target set as

$$\begin{array}{c}{\kappa }_3=\left[p_{{\kappa }_{3}min}={50}^{\circ }/s,p_{{\kappa }_{3}max}={60}^{\circ }/s\right]\times \left[q_{{\kappa }_{3}min}=-5^{\circ }/s,q_{{\kappa }_{3}max}=5^{\circ }/s\right]\hfill \\ \times \left[r_{{\kappa }_{3}min}=-5^{\circ }/s,r_{{\kappa }_{3}max}=5^{\circ }/s\right],\hfill \end{array}$$

$$\begin{array}{c}{\kappa }_4=\left[p_{{\kappa }_{4}min}=5^{\circ }/s,p_{{\kappa }_{4}max}=-5^{\circ }/s\right]\times \left[q_{{\kappa }_{3}min}={50}^{\circ }/s,q_{{\kappa }_{3}max}={60}^{\circ }/s\right]\hfill \\ \times \left[r_{{\kappa }_{3}min}=-5^{\circ }/s,r_{{\kappa }_{3}max}=5^{\circ }/s\right],\hfill \end{array}$$

and

$$\begin{array}{c}{\kappa }_5=\left[p_{{\kappa }_{5}min}=5^{\circ }/s,p_{{\kappa }_{5}max}=-5^{\circ }/s\right]\times \left[q_{{\kappa }_{3}min}=5^{\circ }/s,q_{{\kappa }_{3}max}=-5^{\circ }/s\right]\hfill \\ \times \left[r_{{\kappa }_{3}min}={50}^{\circ }/s,r_{{\kappa }_{3}max}={60}^{\circ }/s\right]\hfill \end{array}$$

for the computation of the reachable set. The size of these can reflect the maneuvering ability around three axes in the angular rate loop. Take the slices $\zeta ={35}^{\circ }$ of the morphing aircraft’s reachable sets to compare with that of the fixed-geometry aircraft. The results of the computations for $T=0.5$ s are given in Figure 13, Figure 14 and Figure 15.

It can be seen that in the angular rate loop, the reachable set region of the morphing aircraft is significantly larger than that of the fixed geometry aircraft; this is caused by the significant effect of morphing on the aerodynamic moments and the moments of inertia. Therefore, the response speed of the inner loop can be significantly improved by morphing.

According to the computation of the reachable sets of the morphing aircraft and the fixed-geometry aircraft, we found that the response speed of the angular rate loop of the morphing aircraft is significantly faster than that of the fixed-geometry aircraft, while the response speed of the flight path loop of the morphing aircraft is slightly faster than that of the morphing aircraft under the condition of the inner loop dynamic is not considered. However, it is worth pointing out that the faster response speed of the inner loop also has a significant impact on the response speed of the outer loop. Considering the above factors, we can conclude that the aircraft’s maneuverability and response speed can be improved through reasonable morphing.

Remark 5.

Note that the reachable set analysis of the attitude loop is not given. This is because the equations of motion of the morphing aircraft in the attitude loop are identical to those of the fixed-geometry aircraft such that the reachable sets of both are also the same.

5. Conclusions

In this study, we examined the impact of morphing on the short-term dynamic characteristics of the aircraft. We first obtained a decoupled six-degree-of-freedom morphing aircraft model and we then analyzed the effects of morphing on the aircraft’s dynamics through open-loop simulations. Finally, we investigated the effects of morphing on the aircraft’s maneuverability using the reachable set theory. Based on our findings, the following conclusions can be drawn:

• Morphing can significantly impact the dynamics of the aircraft, especially during turns where a change in shape can result in a loss of stability.
• Among the various factors affected by morphing, the change in aerodynamic forces, additional inertia moment, and variation in the additional gravitational moment appear to have the most significant effects on dynamics, while inertial forces are generally negligible compared to the aerodynamic forces.
• Morphing can substantially improve the response speed of the angular velocity loop, while having a smaller impact on the response speed of the flight path loop. As a result, the maneuverability of the aircraft can be significantly improved through morphing.

In our future research, we plan to explore the effects of asymmetric morphing on the short-term dynamics of aircraft using more accurate aerodynamic models. We will focus on the crucial role of asymmetric deformation in generating roll moment and investigate its impact on the aircraft’s small-time local controllability and global controllability. In addition, we intend to examine the effects of morphing in various scenarios by incorporating control laws or typical control sequences. By conducting these analyses, we will gain a more comprehensive understanding of the impact of morphing on modern aircraft performance.