Modeling and Control of an Articulated Multibody Aircraft

Abstract

Insects use dynamic articulation and actuation of their abdomen and other appendages to augment aerodynamic flight control. These dynamic phenomena in flight serve many purposes, including maintaining balance, enhancing stability, and extending maneuverability. The behaviors have been observed and measured by biologists but have not been well modeled in a flight dynamics framework. Biological appendages are generally comparatively large, actuated in rotation, and serve multiple biological functions. Technological moving masses for flight control have tended to be compact, translational, internally mounted and dedicated to the task. Many flight characteristics of biological flyers far exceed any technological flyers on the same scale. Mathematical tools that support modern control techniques to explore and manage these actuator functions may unlock new opportunities to achieve agility. The compact tensor model of multibody aircraft flight dynamics developed here allows unified dynamic and aerodynamic simulation and control of bioinspired aircraft with wings and any number of idealized appendage masses. The demonstrated aircraft model was a dragonfly-like fixed-wing aircraft. The control effect of the moving abdomen was comparable to the control surfaces, with lateral abdominal motion substituting for an aerodynamic rudder to achieve coordinated turns. Vertical fuselage motion achieved the same effect as an elevator, and included potentially useful transient torque reactions both up and down. The best performance was achieved when both moving masses and control surfaces were employed in the control solution. An aircraft with fuselage actuation combined with conventional control surfaces could be managed with a modern optimal controller designed using the multibody flight dynamics model presented here.

1. Introduction

The potential applications of small unmanned aerial vehicles (SUAVs) for both military and civil applications continue to drive interest in their design. Insect-inspired maneuverability and agility are desirable characteristics for small SUAVs that might be expected to operate in cluttered environments [1].

Reproductive efficiency and the survival of many animals depends on locomotion for migration, territorial defense, predation and escape [2]. During aerial locomotion, animals actively change their body posture, and hence mass distribution, to produce forces and moments for propulsion and control [3,4]. Active maneuvering, therefore, relies on the ability of the animal to optimize the benefits of three main physical mechanisms to control posture: aerodynamic lift and drag on wings and body [5,6], control of the center of body mass and thus distance to aerodynamic force vectors [3,7], and body moments of inertia quantities and rates through active changes in body shape [8].

For mass-distribution modification, air and space craft have used internally moving masses as an addition to, or replacement of, traditional aerodynamic and moment generating actuators. Moving masses generate moments due to gravity or change the inertia properties for the purpose of controlling the motion of the vehicle. Most efforts have focused on employing movement of internal masses for trim, attitude, stability margin, detumbling of spinning spacecraft, trajectory, orbit, glide path and depth control.

The biological systems exploiting moving masses have similarities to technological examples and substantial differences. In addition, existing formulations of flight dynamics do not adequately capture the effect of comparatively large swinging appendages (masses) observed in insect flight. A convenient mathematical representation appropriate to aerodynamic systems connected to comparatively large, possibly numerous, mass elements connected by joints will be presented. This formulation allows modern multiple-input-multiple-output (MIMO) control techniques to be implemented to manage winged craft with actuated joints and weighted appendages.

1.1. Related Work

This study draws together two specific parts of the literature. One part replicates mechanisms and control paradigms that have application to small aircraft in the vicinity of insect scale. Another part involves developing flight dynamics and control models for aircraft, using conventional nomenclature and extensions of established technique, to simulate operation and use modern control methods.

1.1.1. Moving Masses in Insect Flight

In recent years, there has been an increased interest in the development of insect-inspired micro-aerial vehicles (MAVs), but biological research into insect flight has a long history. The development of bioinspired systems requires the understanding of the physics of biological systems and behaviors, not just copying the external shape [9].

Prior research on the importance of body appendages for force and moment generation in actively flying animals was largely done on insects, including fruit flies [10,11,12], house flies [13], orchid bees [14], moths [15,16], honeybees [17] and butterflies [18]. Several videos by Rüppell show dragonfly abdominal movements during flight [19]. The study conducted by Bode-Oke [20] also found some abdominal movement from kinematic analysis of dragonfly flight. This theory was tested using visual stimulation that mimicked yaw rotations, in which flying tethered flies bent their rear legs and abdomen horizontally to the inside of the intended turn [10,11,13,21]. In contrast, visual stimulation that mimicked body-pitching caused the abdomen to bend upward during upward motion of the visual pattern [22,23,24]. Abdominal steering is adequate for maintenance of body posture in the hawkmoth, according to mathematical models [15,16,22,25]. Some insects, such as desert locusts, flex their abdomens in response to changing air flow conditions in addition to the visual stimulus [26]. This behavior is said to resemble an aerodynamic rudder that aids the animal in orienting itself into the direction of the wind during flight [27,28]. In the flight of monarch butterflies, abdomen undulation has been shown to improve stability, and reduce energy and power consumption [29,30].

Many flight dynamics models used in studies rely upon a conventional single-body approach [31,32,33] to explore limited aspects of insect flight. Yet, to adequately capture the effects of movement of anatomical parts, the multibody modeling approach provides the necessary insight to insect flight [34]. Various methods have been used to derive multibody equations of motion for air vehicles including Newton–Euler [35,36,37,38], Lagrange’s energy-based methods [29,30,39,40], D’Alembert’s principle [41,42], and Kane’s method [43,44]. Regarding insect-inspired aircraft, most of the multibody flight dynamics models in the literature focus on the degrees of freedom associated with wing motion. Although some strictly only consider the relative movement of the wings [35,43,45,46,47,48,49], a few consider the effects of wing mass and inertia [38,50]. Fewer studies have considered the effect or role of abdominal articulation. A nonlinear model of a dragonfly-like flapping-wing MAV was developed by Du [51]. The study considered the effect of the head and abdomen in developing the equations of motion of a dragonfly-like MAV; however, their movements were not investigated with regard to their effects on flight, stabilization or control of the model. Dhyr et al. [22] developed a model to examine the active role of the abdomen in the hawkmoth. However, the model was obtained using a system identification approach. Another study by Tejaswi et al. [29,30] developed a model of the dynamics of a flapping-wing flyer, inspired by monarch butterflies. The model was, however, developed using Euler–Lagrange equations to examine the use of abdomen undulation to improve stability, and reduce energy and power consumption.

Dynamic moving masses in insects are typically a relatively large part of the aircraft mass, or have quite high moment due to length, or both. Actuation is generally in rotation about a joint with at least two degrees of rotational freedom, or involve controlled flexibility, or might have multiple joints along its length. Actuation can often be rapid, such that reaction torques are also significant. The mass is generally multifunctional, serving both anatomical and flight control functions.

1.1.2. Moving-Mass Control in Spacecraft and Aircraft

Mass distribution plays an important role in vehicle dynamics, and varying the vehicle’s mass distribution affects the system’s dynamics, stability and balance [52,53,54,55,56,57,58,59,60,61,62]. In the early days of aviation, pilots would modify their body posture, thus moving the aircraft’s center of mass to accomplish control or stability [63,64]. The concept is of primary importance in the control of hang-gliders [7].

In [65], Qin and Yang designed an attitude control system for a transatmospheric missile using three internal moving masses. Each mass was able to move along its respective axis to change the attitude of the missile. Simulation results showed the method was able to effectively control the missile, achieving static stability and attitude control. Other applications for missile control can be found in [66,67].

Rogers and Mark [68,69,70] developed a seven-degree-of-freedom (7DoF) flight dynamics model for trajectory control of a projectile with an internal translating vibrating mass to study the potential of this control mechanism. Significant control authority was achieved with a mass of about 1% of the total projectile mass by vibrating the mass normal to the axis of symmetry and at the roll rate frequency. Nonlinear equations of a spinning vehicle with two internal moving elements were derived by Wang et al. [71] to investigate the use of linearly moving-mass actuators for trajectory control. Results showed that the moving-mass controller was able to effectively control the vehicle’s trajectory.

Chen et al. [72] investigated a stratospheric airship with various control systems, including traditional aerodynamic control surfaces, vector thrust, ballonet, and moving mass. One mass moved laterally, controlling roll angle, while the other moved longitudinally to regulate pitch angle. The results showed that moving-mass control (MMC) outperformed control surfaces, since the moving masses had no effect on the drag of the airship.

The use of internal moving masses for flight control in SUAVs has also been demonstrated. Ertuk [36] was able to demonstrate the use of translational internal moving masses as a feasible alternative moment generation mechanism to conventional aerodynamic control surfaces for SUAVs. Another study by Vengate [73] showed through flight test that without ailerons, a small unmanned airplane was able to control roll using two laterally moving internal masses.

Moving-mass control has been implemented in rotorcraft as well. For the control of a large quadrotor, Haus et al. [74,75] implemented four translational movable masses in the quadrotor’s four arms to change the center of gravity and control the pitch and roll of the craft. Bouabdallah [76] used MMC for passive roll and pitch stabilization of an indoor coaxial helicopter. The mechanism involved the use of two semi-circular guides for mass movement. A study presented by Bermes [77] successfully implemented a spherical moving-mass steering mechanism to improve the passive stability of a mini-coaxial helicopter in hover.

The most effective techniques for managing spacecraft center of mass (CoM) and attitude trajectory are reaction forces and torques, using internal mass ejection and rotating masses. Command torques are also generated by reaction forces that do not pass through the CoM. As such, only propulsion, orbit and attitude control can be achieved [78]. Rotating masses, on the other hand, produce a variation in angular momentum in any direction, consequently commanding a torque for attitude control [78]. The masses are rotated by electric motors and the corresponding actuators are either fixed-axis variable-speed motors (reaction wheels) or control moment gyros with gimbaled fixed-speed motors that may be oriented in any desired direction. At low orbital altitudes of less than about 450 km, spacecraft suffer from aerodynamic disturbances. Momentum generated by moving-mass actuation has been observed to counteract the effects of these disturbances and stabilize the spacecraft [79].

Moving masses are attractive for exo-atmospheric situations, where aerodynamic forces are not available. Spacecraft can also suffer from wobbling. A study carried out by Childs [80] identified the movement of astronauts (moving mass) inside the space station as the main source of wobbling motion. Accordingly, Childs designed and implemented a momentum exchange mechanism for attitude stabilization with a single moving mass [81,82]. The mechanism was comprised of two cylindrical rigid bodies axially connected by cables that was effective in generating torque to eliminate wobbling.

When angular velocity is high, docking and reentry of spaceships is unsafe for astronauts. Accidental collisions can cause uncontrolled tumbling, which is also dangerous. Edwards and Kaplan using a movable mass control concept, designed a mechanism to convert tumbling motion into simple spin and demonstrated it through simulation. Results from the study showed that a large space station tumbling after collision was able to stabilize using a movable mass of about 1% of its total mass.

Due to the low density of the upper atmosphere, aerodynamic controllers such as a rudder may be inadequate for reentry vehicles [83]. Because aerodynamic control surfaces are prone to ablation and degradation by hypersonic flow, internal moving-mass actuation has been used to control reentry vehicles for decades [84]. Gao et al. [85] investigated the combined use of moving mass and aileron control. Pitch and yaw were controlled by moving masses along two orthogonal rails, while roll was controlled by ailerons. The study demonstrated the benefits of combining these control methods. Su et al. [86] also demonstrated roll control for a non-axisymmetric reentry vehicle using a moving mass. The moving mass was able to provide robust stability and performance throughout the flight without the need to modify controller gains. For other applications of moving mass in air and spacecraft, see [36]. Fluid transfer is another mass-shifting mechanism that has been used for control. For instance, because the Concorde could fly at supersonic speeds, changing the internal mass distribution by moving fuel helps with trimming the aircraft in transonic conditions [87,88,89].

For aircraft, adjustment of balance has by far been the preferred moving-mass application. The choice of mechanism for moving the masses varies with application. Internal solid mass actuation is done with linear motion in most applications [65,68,70,72,86,90,91,92,93,94] and along a circular path in some other applications [76,77,95].

From the literature on aircraft control, it is clear that any moving-mass control is an exception rather than the norm. For automated systems, mass movement has usually been slow (apart from gyroscope reaction torque-based control), generally dedicated to the function, usually linear, and typically internal.

1.1.3. Control Design for Multibody Aircraft

One of the many attributes that emerged from the application of internal moving masses in vehicles is the concept of multiple effectors. The availability of multiple effectors with overlapping functions introduces the problem of control allocation in an overactuated system. Feedback control system theory forms the backbone to address the challenges associated with MIMO systems. By comparing the reference signal to the actual response of the system, the controller calculates the control effort required [96,97]. The response of an aircraft system can be represented by nonlinear equations of motion with 6DoF, which can further be divided into longitudinal and lateral–directional parts, respectively describing pitching and rolling–yawing motions. Using the dynamics of each of these parts, individual flight controllers can be developed. Analysis of linear systems contain concepts and tools that facilitate control of such systems. Compared to nonlinear system-control theory, linear system-control theory presents a tractable solution and has been used extensively in aircraft flight control [96]. The Proportional Integral Derivative (PID) controller is the most widely used controller in real-world applications, and has been used for flight control as in [98,99,100,101,102,103]. However, the gain of the controller needs to be iteratively tuned, which consumes time and effort. Of the more robust control methods, $H_{\infty }$ control and sliding mode control (SMC) are some popular choices [104,105,106,107,108]. However, $H_{\infty }$ control is conceptually difficult and challenging to implement because it requires in-depth mathematical knowledge and model information. On the other hand, SMC is susceptible to chattering due to switching of the high-frequency control. Dynamic inversion is another common control method used for flight control [98,109]. Although this method is intended for nonlinear systems, estimation of the analytical inversion of nonlinear systems online is very difficult.

Optimal control is another control design approach, which has the advantage of being able to provide the best possible solution for a given cost function. Popular optimal control schemes include linear quadratic regulator (LQR), linear quadratic integrator (LQI), and model predictive control (MPC). Such schemes provide the optimal solution, and can also include multiple loop feedback at the same time [97,102,110,111]. Control Allocation (CA) is another strategy that has been developed to handle overactuated systems, which answers the question of how to distribute control action among a redundant set of effectors to achieve a control design. In control allocation, actuator selection is separated from the controller task so that the controller only handles virtual control [112]. Common methods for solving the control allocation problem include optimization-based CA [113,114,115], direct CA [116] and daisy chaining CA [117,118]. A detailed description of these methods can be found in [112,116,119].

1.2. A Dragonfly Biomimetic Aircraft Model

The conceptual SUAV developed in this study is inspired by the impressive flight characteristics of the dragonfly (of order Odonata, infraorder Anisoptera) [120]. Dragonflies can both hover and achieve high performance and efficiency in forward flight with cruise speeds of up to 54 km/h, with some species documented to migrate up to 800 km [121]. Dragonflies have excellent glide properties, with lift-to-drag ratios ranging from 3.5 to over 10 [122,123,124]. They are characterized by a large head with large eyes, a dense thorax, two pairs of independently controlled wings and a long abdomen. The abdomen of dragonflies constitutes a significant fraction (30–35%) of the total mass compared to the wings, which weigh less than 2% of the total mass. The dragonfly can be considered a tailless aircraft due to the lack of conventionally arranged aileron/elevator/rudder aerodynamic control surfaces for roll/pitch/yaw control. Most tailless aircraft are swept, with some dihedral. Dragonflies in fixed-wing flying mode, however, operate with a nearly straight-wing (zero quarter-chord sweep) configuration and are convenient insects for this study because they arguably best represent conventional aircraft. In many ways, they have made the fewest concessions to terrestrial operation, with limbs so adapted to grappling that they can hardly be used for walking. Long, unencapsulated, permanently horizontal wings sprout from a heavily muscled thorax, with most weight concentrated close to the line of bilateral symmetry. This leads to an anatomical simplicity relative to other flying insects. Comparing a mosquito to a dragonfly in Figure 1, it is apparent that the number of appendages available for the control of a mosquito is high.

1.3. Scope and Contributions

In this paper, an increased understanding of the various abdominal motions for control in insect flight is pursued through mathematical modeling and attitude control of a dragonfly-inspired straight-wing aircraft (DISWA), augmented to include an articulated abdomen (see Figure 2b). The flight of some insects consists of both flapping- and fixed-wing episodes. Since this is an isolated study of the effects of abdominal motions, the approach used will provide insight into the use of appendage motion for control. Technological systems have included dedicated active mass-distribution mechanisms for control, and in most cases, the slowly moving elements have been enclosed within a craft, creating space constraints. In the case of animals, all cases mentioned had comparatively fast, multifunctional, rotational or telescoping, external structures with physiological as well as control functions. This is indicative of the high level of optimization of biological systems, presenting inspiration for the design of technology.

The multibody nonlinear equations of motion were developed using the Newton–Euler method. Flight dynamics involves the use of many reference frames and respective coordinate systems. Changing the reference frame, for example, from inertial to body frame, is sometimes advantageous. The change of frame is governed by Euler’s generalized transformation. Classical Euler transformation factors the change in reference frame using ordinary time derivatives, where additional terms appear in the equations of motion (EoMs) to account for the time-dependent coordinate transformations. Using the classical formulation, the resulting expression is not a tensor due to the extra terms from time-dependent changes in coordinates, resulting in equations that are often described as complex and cumbersome [42,50]. The formulation used in this study preserves the tensor formulation using the rotational time derivative operator, offering a few major advantages. It is tensor-based and compact. It is in an invariant tensor form, which is independent of any coordinate system and therefore allows us to derive the mathematical model first without consideration of coordinate systems. Points and frames are used instead of classical radius vectors and coordinate systems.

Since the fore and hind wings are combined to form single wings with elevons, the equations of motion are simplified to a two-body problem; however, the method can be extended to include more than two bodies. In addition, the model of the articulated abdomen for control was derived for a specific regime, which allowed the linearization of the multibody model for stable equilibrium (level cruising flight). A linear quadratic integral controller (LQI) was designed and implemented to optimally maintain attitude (pitch and yaw) of the DISWA. The LQI controller offers the advantage of simplicity in design, and application of multiple-input-multiple-output (MIMO) systems in terms of control effort allocation. The developed model will enable the study of the flight dynamics of insect-like flapping aerial vehicles and a determination of the relative importance of abdominal actuation to insect flight.

The model we will derive treats appendages as mass elements that do not create lift or drag. Depending on the appendage, this might be a reasonable approximation. Future work could consider moments from drag on appendages. In this first formulation, it would replace one approximation (dragless slender bodies) with another (drag on slender bodies).

2. Development of the Multibody Equations of Motion

In this section, the equations of motion (EoMs) for the DISWA are derived using the simplest case of a two-body model. The dragonfly head/thorax/wings combination is abstracted as a rigid body and will be referred to as the central body. The central body typically has six degrees of freedom—three translational and three rotational. The abdomen is the second rigid body which has three additional degrees of freedom, constrained to move with respect to the central body about a joint. Additionally, the abdomen will be interchangeably referred to as the tail in this paper. It is common practice to derive equations of motion referenced to the combined center of gravity ($cg$) of an aircraft; however, to properly reflect the effects associated with dynamic changes in combined $cg$ position, the dynamic equations presented in this study are referenced to point b, which is the center of mass of the central body.

2.1. Preliminaries

This section briefly introduces the concept of the rotational time derivative operator, which is used throughout the text. A brief description is presented below; more details about the rotational time derivative and its application in modeling can be found in [125]. Consider two arbitrary frames A and O, with relative angular velocity of ${\omega }^{OA}$, and a vector a, associated with frame A. For clarity, a vector associated with a frame implies it is constant, as it is fixed relative to the associated frame. In addition, vector a transforms similar to a first-order tensor given the coordinate systems of frames A and O as in

$${\left[a\right]}^O={\left[R\right]}^{OA}{\left[a\right]}^A$$

(1)

where ${\left[R\right]}^{OA}$ is the transformation matrix from the coordinate system of frame A to frame O. Transformation matrices are applicable to second-order tensors (equivalent to matrices) as well.

Given that the time rate of change of vector a with respect to frame A is given by ${\left[\frac{da}{dt}\right]}^A$, the time rate of change of vector a, expressed in frame O, can be obtained by applying the classical Euler transformation using ordinary time derivatives as

$${\left[\frac{da}{dt}\right]}^A={\left[\frac{da}{dt}\right]}^O+{\Omega }^{OA}a$$

(2)

where ${\Omega }^{OA}$ is the skew symmetric matrix of the angular velocity vector ${\omega }^{OA}$. The resulting expression on the right-hand side (RHS) of Equation (2) does not preserve the tensorial formulation due to the time-dependent nature of the derivative of vector a.

To preserve the tensor formulation during Euler transformation, a time operator, called the rotational time derivative, is used instead of the ordinary time derivative, and was first introduced in Wrede’s textbook [126] on vector and tensor analysis. The operator depends only on a frame of reference such that the rotational time derivative of a vector a regarding frame A is expressed as $D^{A}a$. A change in reference frame of vector a from frame A to O is obtained through Euler transformation [125] with the rotational time derivative, using the following relationship

$$D^{A}a=D^{O}a+{\Omega }^{OA}a$$

(3)

Thus, the vector a can be expressed in the coordinate system of any frame due to the tensor properties the rotational time derivative possesses as in

$${\left[D^{A}a\right]}^O={\left[R\right]}^{OA}{\left[D^{A}a\right]}^A$$

(4)

The time rate of change of vector a with respect to its associated frame A can be written as ${\left[D^{A}a\right]}^A$. It is important to note that rotational time derivatives expressed in their own associated frame become ordinary time derivatives ready for numerical integration as in

$${\left[D^{A}a\right]}^A=\frac{da}{dt}.$$

(5)

In the case of a multibody aircraft, the use of fixed body reference frames and displacement vectors between points in the associated reference frames enables the equations to be derived in tensor form, and to be both compact and coordinate-independent (invariant).

2.2. Reference Frames and Coordinate Systems

The reference frames and associated coordinate systems used for the development of equations of motion are as follows and shown in Figure 2b.

(a)

The inertial reference frame $I\{X_I,Y_I,Z_I\}$: The Earth frame is assumed to be the inertial frame with its origin, I fixed at an arbitrary point relative to the Earth’s surface. The orientation of the inertial frame is such that the $X_I$ axis is positive facing North, $Y_I$ axis is positive facing East and $Z_I$ axis is positive downwards towards Earth’s center of gravity.

(b)

The body-fixed reference frame $B\{x_B,y_B,z_B\}$: The origin is located at the center of mass of the central body, point b. The body frame is oriented such that $x_B$ axis lies on the plane of symmetry of the aircraft and points in the forward direction towards the head of the aircraft. The $y_B$ axis is perpendicular to the $x_B$ axis, pointing towards the right side of the aircraft and the $z_B$ axis is positive downwards and lies on the plane of symmetry.

(c)

The abdominal/tail reference frame $T\{x_T,y_T,z_T\}$: This reference frame originates from the center of mass of the tail with its orientation the same as that of the body frame coordinate system when the tail is not deflected.

2.3. Reference Points

Reference points are introduced for each individual part of the aircraft model and are shown in Figure 2b. The reference point for the central body is its center of mass, point b, associated with frame B. The reference point for the tail is its center of mass, point t, associated with frame T. The tail joint, j, serves as an intermediate reference point. It will be mostly associated with frame B except in cases where it is stated to be associated with frame T. The center of mass of the whole aircraft (rigid body C) is located at point c and is associated with frame B.

2.4. Orientation and Transformation Matrices

The numerical simulation framework will rely on the appropriate transformation matrices to ultimately yield the EoMs in the body-fixed coordinate system. The standard transformation matrices of a rotation angle $\mu$, about x, y and z axes of the coordinate systems are expressed as

$$R_{{\mu |}_x}=\left[\begin{array}{ccc}1& 0& 0\\ 0& c\mu & s\mu \\ 0& -s\mu & c\mu \end{array}\right]R_{{\mu |}_y}=\left[\begin{array}{ccc}c\mu & 0& -s\mu \\ 0& 1& 0\\ s\mu & 0& c\mu \end{array}\right]R_{{\mu |}_z}=\left[\begin{array}{ccc}c\mu & s\mu & 0\\ -s\mu & c\mu & 0\\ 0& 0& 1\end{array}\right]$$

(6)

where $s\mu$ and $c\mu$ represent the sin and cosine functions of $\mu$, respectively. Given any two arbitrary coordinate systems, the transformation matrices can be obtained using relationships between the rotation angles and axes of rotation in Equation (6).

In this paper, the attitude of the central body with respect to the inertial frame is described by the Euler angles $\psi$ (yaw), $\theta$ (pitch) and $\varphi$ (roll). The transformation matrix from inertial to body coordinate system, ${\left[R\right]}^{BI}$, is obtained using the z-y-x rotational sequence of the Euler angles as in

$${\left[R\right]}^{BI}=R_{{\psi |}_z}{R}_{{\theta |}_y}{R}_{{\varphi |}_x}$$

(7)

The transformation matrix from body to inertial coordinate system is obtained as the transpose of the ${\left[R\right]}^{BI}$ matrix, ${\left[R\right]}^{IB}={\left[R{}^{\prime }\right]}^{BI}$. Euler angles ${\psi }_T$ (tail yaw), ${\theta }_T$ (tail pitch) and ${\varphi }_T$ (tail roll) are used to specify the orientation of the abdomen relative to the central body. The transformation matrix from the body to tail coordinate system also follows the z-y-x rotational sequence as in

$${\left[R\right]}^{TB}=R_{{\psi }_T{|}_z}{R}_{{\theta }_T{|}_y}{R}_{{\varphi }_T{|}_x}$$

(8)

2.5. Kinematics

The velocities of each of the bodies will be defined in this section using the rotational time derivative operator and displacement vector between points. The velocities are ultimately expressed in the preferred B frame. Given that the displacement vector of the central body relative to the inertial frame is $s_{{}_{bI}}$, the translational velocity of the central body regarding the inertial frame is expressed as $D^I{s}_{{}_{bI}}$. The translational velocity of the central body in the body frame can then be obtained using the inertial to body transformation matrix

$${\left[D^B{s}_{{}_{bI}}\right]}^B={\left[R\right]}^{BI}{\left[D^I{s}_{{}_{bI}}\right]}^I$$

(9)

where the term ${\left[D^B{s}_{{}_{bI}}\right]}^B=[u,v,w]$, and $u,v,w$ are the body velocity components.

Using the abdominal joint j as the reference point, the displacement of the joint with respect to the inertial frame can be expressed as a summation of displacement vectors as in

$$s_{{}_{jI}}=s_{{}_{jb}}+s_{{}_{bI}}$$

(10)

Hence, the relative translational velocity of the joint regarding the inertial frame is expressed as

$$D^I{s}_{{}_{jI}}=D^I{s}_{{}_{jb}}+D^I{s}_{{}_{bI}}$$

(11)

Since $D^I{s}_{{}_{jb}}$ is not readily available for measurement, a simpler expression is established by applying a Euler transformation of $D^I{s}_{{}_{jb}}$. Equation (11) becomes

$$D^I{s}_{{}_{jI}}=D^B{s}_{{}_{jb}}+{\Omega }^{BI}{s}_{{}_{jb}}+D^I{s}_{{}_{bI}}$$

(12)

Noting that $D^B{s}_{{}_{jb}}=0$, since $s_{{}_{jb}}$ is associated with frame B, simplifies Equation (12) to

$${\left[D^I{s}_{{}_{jI}}\right]}^I={\Omega }^{BI}{s}_{{}_{jb}}+D^I{s}_{{}_{bI}}$$

(13)

Similar to Equation (9), the abdominal joint relative velocity can be obtained using the inertial to body transformation matrix.

Rotational kinematic equations are mathematical representations of the relationships between the rate of change of kinematic parameters and the angular velocities. Let the Euler angle derivatives for the central body and the abdomen be ($\dot{\varphi },\dot{\theta },\dot{\psi }$) and (${\dot{\varphi }}_T,{\dot{\theta }}_T,{\dot{\psi }}_T$), respectively. The angular velocity of the central body relative to the inertial frame, expressed in the B frame, is introduced as ${\left[{\omega }^{BI}\right]}^B=[p,q,r]$, while that of the abdomen relative to the body frame and expressed in the T frame is introduced as ${\left[{\omega }^{TB}\right]}^T=[p_{{}_T},q_{{}_T},r_{{}_T}]$. The relationships between the rate of change of kinematic parameters and the angular velocities follow the sequence of rotations and are calculated according to

$${\left[{\omega }^{BI}\right]}^B=\left[\begin{array}{c}\dot{\varphi }\\ 0\\ 0\end{array}\right]+R_{{\varphi |}_x}\left[\begin{array}{c}0\\ \dot{\theta }\\ 0\end{array}\right]+R_{{\varphi |}_x}{R}_{{\theta |}_y}\left[\begin{array}{c}0\\ 0\\ \dot{\psi }\end{array}\right]$$

(14)

$${\left[{\omega }^{TB}\right]}^T=\left[\begin{array}{c}{\dot{\varphi }}_T\\ 0\\ 0\end{array}\right]+R_{{\varphi }_T{|}_x}\left[\begin{array}{c}0\\ {\dot{\theta }}_T\\ 0\end{array}\right]+R_{{\varphi }_T{|}_x}{R}_{{\theta }_T{|}_y}\left[\begin{array}{c}0\\ 0\\ {\dot{\psi }}_T\end{array}\right]$$

(15)

2.6. Multibody Equations of Motion

Next, the equations of motion (EoMs) are derived using Newtonian and Eulerian mechanics. The equations will be derived using points defined in their associated reference frames and related using invariant displacement vectors. The approach presented in this section allows the tensor nature of the equations to be preserved. EoMs are required for translational, rotational, and abdominal motion.

2.6.1. Translational Dynamics

For a system of k clustered rigid bodies, the translational dynamic equations are developed as the sum of their individual contributions, considering their respective centers of mass at $B_k$, with respect to the inertial frame I as follows [38,125]:

$$D^I\left(\sum _k{m}^{B_k}{D}^I{s}_{{}_{B_{k}I}}\right)=\sum _k{F}_k$$

(16)

where ${\displaystyle \sum _k}{F}_k$ is the contributing sum of all the forces acting on the aircraft from each rigid body, which will be discussed in Section 2.7.

Recall three previously introduced reference points b, j and t. Separating the individual contributions of the central body and the tail with $s_{{}_{B_{k}I}}=s_{{}_{bI}}+s_{{}_{tI}}$ and $s_{{}_{tI}}=s_{{}_{tj}}+s_{{}_{jb}}+s_{{}_{bI}}$ gives

$$D^I\left(m^B{D}^I{s}_{{}_{bI}}+m^T{D}^I{s}_{{}_{tj}}+m^T{D}^I{s}_{{}_{jb}}+m^T{D}^I{s}_{{}_{bI}}\right)=\sum _k{F}_k$$

(17)

The first term on the left-hand side of Equation (17) represents the contribution of the central body, while the second, third and fourth terms represent the contributions of the tail. Merging the first and fourth terms of the left-hand side of Equation (17) and introducing $D^I{s}_{{}_{bI}}=V_b^I$ results in

$$\underset{\mathrm{Term}\mathrm{I}}{\underbrace{m{D}^I{V}_b^I}}+\underset{\mathrm{Term}\mathrm{II}}{\underbrace{m^T{D}^I{s}_{{}_{tj}}}}+\underset{\mathrm{Term}\mathrm{III}}{\underbrace{m^T{D}^I{D}^I{s}_{{}_{jb}}}}=\underset{\mathrm{Term}\mathrm{IV}}{\underbrace{\sum _k{F}_k}}$$

(18)

Since $TermI$ of Equation (18) relates to the body frame, we use the Euler transformation to represent this term in the body-fixed frame as

$$m{D}^I{V}_b^I=m(D^B{V}_b^I+{\Omega }^{BI}{V}_b^I)$$

(19)

$TermII$ of Equation (18) relates to the tail frame, therefore, the rotational time derivative $D^I{s}_{{}_{tj}}$ is transformed to the tail frame, using Equation (3) as in

$$m^T{D}^I{D}^I{s}_{{}_{tj}}=m^T{D}^I(D^T{s}_{{}_{tj}}+{\Omega }^{TI}{s}_{{}_{tj}})$$

(20)

The vector $s_{{}_{tj}}$ is constant in the tail frame, hence, $D^T{s}_{{}_{tj}}=0$ in Equation (20). Additionally, rewriting ${\Omega }^{TI}={\Omega }^{TB}+{\Omega }^{BI}$ results in

$$m^T{D}^I{D}^I{s}_{{}_{tj}}=m^T{D}^I({\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{BI}{s}_{{}_{tj}})$$

(21)

In addition, the rotational time derivatives of ${\Omega }^{TB}$ and ${\Omega }^{BI}$ should be transformed into their associated frames, which are the tail and body frames, respectively. Therefore, Equation (21) becomes

$$\begin{array}{cc}\hfill m^T{D}^I{D}^I{s}_{{}_{tj}}=& m^T[D^T\left({\Omega }^{TB}{s}_{{}_{tj}}\right)+{\Omega }^{TI}{\Omega }^{TB}{s}_{{}_{tj}}+D^B\left({\Omega }^{BI}{s}_{{}_{tj}}\right)+{\Omega }^{BI}{\Omega }^{BI}{s}_{{}_{tj}}]\hfill \\ \hfill & \begin{array}{c}=m^T\left[\begin{array}{c}{D}^T{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{TB}{D}^T{s}_{{}_{tj}}+{\Omega }^{TI}{\Omega }^{TB}{s}_{{}_{tj}}\\ +D^B\left({\Omega }^{BI}{s}_{{}_{tj}}\right)+{\Omega }^{BI}{D}^B{s}_{{}_{tj}}+{\Omega }^{BI}{\Omega }^{BI}{s}_{{}_{tj}}\end{array}\right]\end{array}\hfill \end{array}$$

(22)

Again, since $D^T{s}_{{}_{tj}}=0$, the second term on the RHS of Equation (22) vanishes. Additionally, rewriting ${\Omega }^{TI}={\Omega }^{TB}+{\Omega }^{BI}$ and using the Euler transformation $D^B{s}_{{}_{tj}}=D^W{s}_{{}_{tj}}+{\Omega }^{TB}{s}_{{}_{tj}}$ results in

$$\begin{array}{cc}\hfill m^T{D}^I{D}^I{s}_{{}_{tj}}=& m^T\left[\begin{array}{c}{D}^T{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{TB}{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{BI}{D}^T{s}_{{}_{tj}}+D^B{\Omega }^{BI}{s}_{{}_{tj}}\\ +2{\Omega }^{BI}{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{BI}{\Omega }^{BI}{s}_{{}_{tj}}\end{array}\right]\hfill \end{array}$$

(23)

However, $D^T{s}_{{}_{tj}}=0$; therefore, the final expression for $TermII$ in Equation (18) is

$$\begin{array}{cc}\hfill m^T{D}^I{D}^I{s}_{{}_{tj}}=& m^T\left[\begin{array}{c}{D}^T{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{TB}{\Omega }^{TB}{s}_{{}_{tj}}\\ +D^B{\Omega }^{BI}{s}_{{}_{tj}}+2{\Omega }^{BI}{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{BI}{\Omega }^{BI}{s}_{{}_{tj}}\end{array}\right]\hfill \end{array}$$

(24)

$TermIII$ of Equation (18) should be transformed to the body frame as

$$m^T{D}^I{D}^I{s}_{{}_{jb}}=m^T{D}^I(D^B{s}_{{}_{jb}}+{\Omega }^{BI}{s}_{{}_{jb}})$$

(25)

Considering $D^B{s}_{{}_{jb}}=0$ since $s_{{}_{jb}}$ is fixed in the body frame and transforming the rotational derivative to the body frame gives

$$m^T{D}^I{D}^I{s}_{{}_{jb}}=m^T[D^B{\Omega }^{BI}{s}_{{}_{jb}}+{\Omega }^{BI}{\Omega }^{BI}{s}_{{}_{jb}}]$$

(26)

Finally, $TermIV$ of Equation (18) represents the sum of all external forces acting on the aircraft. Substituting Equation (19), (24) and (26) into Equation (18) gives the translational dynamic equations of motion in an invariant tensor form as

$$\begin{array}{cc}\hfill m(D^B{V}_b^I+{\Omega }^{BI}{V}_b^I)+m^T\left[\begin{array}{c}{D}^T{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{TB}{\Omega }^{TB}{s}_{{}_{tj}}+D^B{\Omega }^{BI}(s_{{}_{tj}}+s_{{}_{jb}})\\ +{\Omega }^{BI}{\Omega }^{BI}(s_{{}_{tj}}+s_{{}_{jb}})+2{\Omega }^{BI}{\Omega }^{TB}{s}_{{}_{tj}}\end{array}\right]=& \sum _k{F}_k\hfill \end{array}$$

(27)

To use this equation in a simulation framework, all terms in Equation (27) must be written as body coordinate frames using transformation matrices.

$$\begin{array}{cc}\hfill & m\left({\left[{\dot{V}}_B^I\right]}^B+{\left[{\Omega }^{BI}\right]}^B{\left[V_b^I\right]}^B\right)+m^T({\left[R\right]}^{BT}\left({\left[{\dot{\Omega }}^{TB}\right]}^T{\left[s_{tj}\right]}^T+{\left[{\Omega }^{TB}\right]}^T{\left[{\Omega }^{TB}\right]}^T{\left[s_{tj}\right]}^T)\right)\hfill \\ \hfill & +{\left[{\dot{\Omega }}^{BI}\right]}^B({\left[R\right]}^{BT}{\left[s_{tj}\right]}^T+{\left[s_{jb}\right]}^B)+\left({\left[{\Omega }^{BI}\right]}^B{\left[{\Omega }^{BI}\right]}^B({\left[R\right]}^{BT}{\left[s_{tj}\right]}^T+{\left[s_{jb}\right]}^B)\right)\hfill \\ \hfill & +\left(2{\Omega }^{BI}{\left[R\right]}^{BT}\left({\left[{\Omega }^{TB}\right]}^T{\left[s_{tj}\right]}^T\right)\right))=\sum _k{\left[F_k\right]}^B\hfill \end{array}$$

(28)

2.6.2. Rotational Dynamics

For a system of k clustered rigid bodies, the rotational dynamic equations are developed as the sum of their individual contributions based on Euler’s law, as follows [38,125]:

$$\underset{\mathrm{Term}\mathrm{I}}{\underbrace{\sum _k{D}^I\left(I_{B_k}^{B_k}{\omega }^{{}_{B_{k}I}}\right)}}+\underset{\mathrm{Term}\mathrm{II}}{\underbrace{\sum _k{D}^I\left(m^{B_k}{S}_{{}_{B_{k}I}}{D}^I{s}_{{}_{B_{k}I}}\right)}}=\sum _k{M}_k+\sum _k{S}_{{}_{B_{k}I}}{F}_k$$

(29)

where the expression ${\displaystyle {\sum }_k}{M}_k+{\displaystyle {\sum }_k}{S}_{{}_{B_{k}I}}{F}_k$ on the right-hand side (RHS) of Equation (29) is the contributing sum of all the moments acting on the aircraft from each rigid body, which will be discussed in Section 2.7. For the system being considered in this study, the tail and body contributions are separated in $TermI$ in Equation (29) as

$$\begin{array}{c}\hfill \sum _k{D}^I\left(I_{B_k}^{B_k}{\omega }^{{}_{B_{k}I}}\right)=D^I\left(I_B^B{\omega }^{{}_{BI}}\right)+D^I\left(I_T^T{\omega }^{{}_T}\right)\end{array}$$

(30)

Applying the chain rule and Euler transformation to each term on the RHS of Equation (30) leads to

$$\begin{array}{c}\hfill \sum _k{D}^I\left(I_{B_k}^{B_k}{\omega }^{{}_{B_{k}I}}\right)={\omega }^{BI}{D}^B{I}_b^B+I_b^B{D}^B{\omega }^{BI}+{\Omega }^{BI}{I}_b^B{\omega }^{BI}\\ \hfill +{\omega }^{TI}{D}^T{I}_t^T+I_t^T{D}^T{\omega }^{TI}+{\Omega }^{TI}{I}_t^T{\omega }^{TI}\end{array}$$

(31)

Since the moment of inertia of the body and wing are constant in their corresponding frame, the first and fourth terms on the RHS of Equation (31) vanish. Additionally, considering $D^T{\omega }^{TI}=D^T{\omega }^{TB}+D^T{\omega }^{BI}$ and ${\Omega }^{BT}=-{\Omega }^{TB}$, and performing the Euler transformation to establish $D^T{\omega }^{BI}$ in the body frame results in

$$\begin{array}{cc}\hfill \sum _k{D}^I\left(I_{B_k}^{B_k}{\omega }^{{}_{B_{k}I}}\right) =& I_b^B{D}^B{\omega }^{BI}+{\Omega }^{BI}{I}_b^B{\omega }^{BI}\hfill \\ \hfill & +I_t^T\left(D^T{\omega }^{TB}+D^B{\omega }^{BI}-{\Omega }^{TB}{\omega }^{BI}\right)\hfill \\ \hfill & +\left(({\Omega }^{TB}+{\Omega }^{BI})I_t^T({\omega }^{TB}+{\omega }^{BI})\right)\hfill \end{array}$$

(32)

Applying the chain rule to $TermII$ in Equation (29) yields

$$\sum _k{D}^I\left(m^{B_k}{S}_{{}_{B_{k}I}}{D}^I{s}_{{}_{B_{k}I}}\right) =\sum _k{m}^{B_k}\left(D^I{S}_{{}_{B_{k}I}}{D}^I{s}_{{}_{B_{k}I}}+S_{{}_{B_{k}I}}{D}^I{D}^I{s}_{{}_{B_{k}I}}\right)$$

(33)

The cross-product of the vector derivatives is zero; therefore, the first term on the RHS of Equation (33) is zero. Additionally, the body and abdominal contributions are separated, leading to

$$\sum _k{D}^I\left(m^{B_k}{S}_{{}_{B_{k}I}}{D}^I{s}_{{}_{B_{k}I}}\right)=m^B{S}_{{}_{BI}}{D}^I{D}^I{s}_{{}_{bI}}+m^T{S}_{{}_{tI}}{D}^I{D}^I{s}_{{}_{tI}}$$

(34)

Using Euler transformation, the first term on the RHS of Equation (34) becomes

$$m^B{S}_{{}_{bI}}{D}^I{D}^I{s}_{{}_{bI}}=m^B{S}_{{}_{bI}}(D^B{V}_b^I+{\Omega }^{BI}{V}_b^I)$$

(35)

The second term on the RHS of Equation (34) could be developed using the vector expansion $s_{{}_{tI}}=s_{{}_{tj}}+s_{{}_{jb}}+s_{{}_{bI}}$ and performing Euler transformation, results in

$$m^T{S}_{{}_{tI}}{D}^I{D}^I{s}_{{}_{tI}}=m^T{S}_{{}_{tI}}{D}^I\left(D^T{s}_{{}_{tj}}+{\Omega }^{TI}{s}_{{}_{tj}}+D^B{s}_{{}_{jb}}+{\Omega }^{BI}{s}_{{}_{jb}}+V_b^I\right)$$

(36)

As $s_{{}_{tj}}$ and $s_{{}_{jb}}$ are constant in their corresponding frames, the first and third terms of the equation above vanish. Additionally, the rotational time derivatives of the last two terms can be transformed to the body frame as follows

$$D^I\left({\Omega }^{BI}{s}_{{}_{jb}}\right)=D^B{\Omega }^{BI}{s}_{{}_{jb}}+{\Omega }^{BI}{\Omega }^{BI}{s}_{{}_{jb}}\mathrm{and}{D}^I{V}_b^I=D^B{V}_b^I+{\Omega }^{BI}{V}_b^I$$

(37)

Only the rotational derivative of the second term on the RHS of Equation (36) needs to be calculated. Using the expansion of the angular velocity tensor ${\Omega }^{TI}={\Omega }^{TB}+{\Omega }^{BI}$ and performing the Euler transformation results in

$$D^I\left({\Omega }^{TI}{s}_{{}_{tj}}\right)=D^T{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{TI}{\Omega }^{TB}{s}_{{}_{tj}}+D^B\left({\Omega }^{BI}{s}_{{}_{tj}}\right)+{\Omega }^{BI}{\Omega }^{BI}{s}_{{}_{tj}}$$

(38)

${\Omega }^{TI}$ in the second term of RHS the above equation should be substituted using ${\Omega }^{TI}={\Omega }^{TB}+{\Omega }^{BI}$. For the third term on the RHS of the above equation, we take the rotational time derivative of the term, apply the Euler transformation, $D^B{s}_{{}_{tj}}=D^T{s}_{{}_{tj}}+{\Omega }^{TB}{D}^B{s}_{{}_{tj}}$ and consider that $D^T{s}_{{}_{tj}}=0$. These simplify the third term to

$$D^B\left({\Omega }^{BI}{s}_{{}_{tj}}\right)=D^B\left({\Omega }^{BI}{s}_{{}_{tj}}\right)+{\Omega }^{TI}{\Omega }^{TB}{s}_{{}_{tj}}$$

(39)

Substituting Equation (39) into Equation (38) and then substituting Equations (37) and (38) into Equation (36) gives

$$\begin{array}{cc}\hfill m^T{S}_{{}_{tI}}{D}^I{D}^I{s}_{{}_{tI}}=& m^T{S}_{{}_{tI}}\left[\begin{array}{c}{D}^T{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{TB}{\Omega }^{TB}{s}_{{}_{tj}}\\ +D^B{\Omega }^{BI}{s}_{{}_{tj}}+2{\Omega }^{BI}{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{BI}{\Omega }^{BI}{s}_{{}_{tj}}\\ +D^B{\Omega }^{BI}{s}_{{}_{jb}}+{\Omega }^{BI}{\Omega }^{BI}{s}_{{}_{jb}}+D^B{V}_b^I+{\Omega }^{BI}{V}_b^I\end{array}\right]\hfill \end{array}$$

(40)

Finally, substituting Equations (35) and (40) into Equation (34) and further substituting Equations (32) and (34) into Equation (29), we obtain

$$\begin{array}{cc}\hfill & m^T{S}_{{}_{tI}}\left[\begin{array}{c}{D}^T{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{TB}{\Omega }^{TB}{s}_{{}_{tj}}\\ +D^B{\Omega }^{BI}{s}_{{}_{tj}}+2{\Omega }^{BI}{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{BI}{\Omega }^{BI}{s}_{{}_{tj}}\\ +D^B{\Omega }^{BI}{s}_{{}_{jb}}+{\Omega }^{BI}{\Omega }^{BI}{s}_{{}_{jb}}+D^B{V}_b^I+{\Omega }^{BI}{V}_b^I\end{array}\right]\hfill \\ \hfill & +I_t^T\left(D^T{\omega }^{TB}+D^B{\omega }^{BI}-{\Omega }^{TB}{\omega }^{BI}\right)+\left(({\Omega }^{TB}+{\Omega }^{BI})I_t^T({\omega }^{TB}+{\omega }^{BI})\right)\hfill \\ \hfill & +I_b^B{D}^B{\omega }^{BI}+{\Omega }^{BI}{I}_b^B{\omega }^{BI}+m^B{S}_{{}_{bI}}(D^B{V}_b^I+{\Omega }^{BI}{V}_b^I)=\sum _k{M}_k+\sum _k{S}_{{}_{B_{k}I}}{F}_k\hfill \end{array}$$

(41)

The second term on the RHS of this equation can be written as the sum of contributions from individual parts

$$\begin{array}{ccc}\hfill \sum _k{S}_{{}_{b_{k}I}}{F}_k=& S_{{}_{bI}}{F}_B+\left(S_{{}_{TB}}+S_{{}_{bI}}\right)F_T\hfill & \hfill =S_{{}_{bI}}\sum _k{F}_k+S_{{}_{tb}}{F}_T\end{array}$$

(42)

The expression for ${\displaystyle {\sum }_k}{F}_k$ can be obtained from Equation (27). Multiplying both sides of Equation (27) by $S_{{}_{bI}}$ and inserting into Equation (34) results in corresponding terms being canceled from both sides. Therefore, the rotational dynamic equations of motion in an invariant tensor form are

$$\begin{array}{cc}\hfill & m^T{S}_{{}_{tb}}\left[\begin{array}{c}{D}^T{\Omega }^{TB}{s}_{{}_{tj}}+{\Omega }^{TB}{\Omega }^{TB}{s}_{{}_{tj}}\\ +D^B{\Omega }^{BI}(s_{{}_{tj}}+s_{{}_{jb}})+{\Omega }^{BI}{\Omega }^{BI}(s_{{}_{tj}}+s_{{}_{jb}})\\ 2{\Omega }^{BI}{\Omega }^{TB}{s}_{{}_{tj}}+D^B{V}_b^I+{\Omega }^{BI}{V}_b^I\end{array}\right]+I_t^T\left(D^T{\omega }^{TB}+D^B{\omega }^{BI}-{\Omega }^{TB}{\omega }^{BI}\right)\hfill \\ \hfill & +\left(({\Omega }^{TB}+{\Omega }^{BI})I_t^T({\omega }^{TB}+{\omega }^{BI})\right)+I_b^B{D}^B{\omega }^{BI}+{\Omega }^{BI}{I}_b^B{\omega }^{BI}=\sum _k{M}_k+S_{{}_{tb}}{F}_T\hfill \end{array}$$

(43)

For numerical simulations, these equations must be expressed in the body frame coordinate system as follows

$$\begin{array}{cc}\hfill & m^T({\left[R\right]}^{BT}{\left[S_{tj}\right]}^T+{\left[S_{jb}\right]}^B)({\left[R\right]}^{BT}\left({\left[{\dot{\Omega }}^{TB}\right]}^T{\left[s_{tj}\right]}^T+{\left[{\Omega }^{TB}\right]}^T{\left[{\Omega }^{TB}\right]}^T{\left[s_{tj}\right]}^T\right)\hfill \\ \hfill & +{\left[{\dot{\Omega }}^{BI}\right]}^B({\left[R\right]}^{BT}{\left[s_{tj}\right]}^T+{\left[s_{jb}\right]}^B)+\left({\left[{\Omega }^{BI}\right]}^B{\left[{\Omega }^{BI}\right]}^B({\left[R\right]}^{BT}{\left[s_{tj}\right]}^T+{\left[s_{jb}\right]}^B)\right)\hfill \\ \hfill & +\left(2{\Omega }^{BI}{\left[R\right]}^{BT}\left({\left[{\Omega }^{TB}\right]}^T{\left[s_{tj}\right]}^T\right)\right)+\left({\left[{\dot{V}}_B^I\right]}^B+{\left[{\Omega }^{BI}\right]}^B{\left[V_b^I\right]}^B\right))\hfill \\ \hfill & +{\left[R\right]}^{BT}{\left[I_t^T\right]}^T\left({\left[{\dot{\omega }}^{TB}\right]}^T+{\left[R{}^{\prime }\right]}^{BT}{\dot{\omega }}^{BI}-{\left[{\Omega }^{TB}\right]}^T{\left[R{}^{\prime }\right]}^{BT}{\left[{\omega }^{BI}\right]}^B\right)\hfill \\ \hfill & +\left(({\left[R\right]}^{BT}{\left[{\Omega }^{TB}\right]}^T+{\left[{\Omega }^{BI}\right]}^B){\left[R\right]}^{BT}{\left[I_t^T\right]}^T({\left[{\omega }^{TB}\right]}^T+{\left[R{}^{\prime }\right]}^{BT}{\left[{\omega }^{BI}\right]}^B)\right)\hfill \\ \hfill & +{\left[I_b^B\right]}^B{\left[{\dot{\omega }}^{BI}\right]}^B+{\left[{\Omega }^{BI}\right]}^B{\left[I_b^B\right]}^B{\left[{\omega }^{BI}\right]}^B=\sum _k{\left[M_k\right]}^B+\left[{\left[S_{{}_{tb}}\right]}^B{\left[F_T\right]}^B\right]\hfill \end{array}$$

(44)

2.6.3. Abdominal Motion

The 3 DoF actuation of the abdomen can be considered to be three revolute joints for each DoF. Consider the pitch angular motion ${\theta }_T$ of the abdomen about point j, relative to the body. The corresponding pitch angular motion is given by

$$I_{jy}^T{\ddot{\theta }}_T+\kappa {\dot{\theta }}_T={\tau }_{ly}+{\tau }_{ay}$$

(45)

where $I_{jy}^T$ is the mass moment of inertia of the abdomen about y-axis at the joint j, and $\kappa$ is the damping constant. ${\tau }_{ly}$ is the load torque due to gravity and ${\tau }_{ay}$ is the torque applied at the joint to actuate the abdomen.

2.7. Forces and Moments

The resultant forces and moments acting on the aircraft, expressed as contributions from the central body and the abdomen, are represented by the RHS of Equations (28) and (43) respectively. The forces and moments acting on the aircraft are due to aerodynamics ($F_A,M_A$), gravity ($F_G,M_G$), and propulsion ($F_P,M_P$), which ultimately must be represented in the body frame coordinate system as follows

$$\sum _k{\left[F_k\right]}^B={\left[F_A\right]}^B+{\left[F_G\right]}^B+{\left[F_P\right]}^B$$

(46)

$$\sum _k{\left[M_k\right]}^B+\left[{\left[s_{{}_{tb}}\right]}^B\times {\left[F_T\right]}^B\right]={\left[M_A\right]}^B+{\left[M_G\right]}^B$$

(47)

The propulsive force is assumed to act parallel to the aircraft $x_B$ axis and produces no moments, which yields

$${\left[F_P\right]}^B=\left[\begin{array}{c}{T}_n\\ 0\\ 0\end{array}\right]$$

(48)

The resultant gravitational forces are estimated by summing the contributions from the central body and the abdomen

$${\left[F_G\right]}^B={\left[R\right]}^{BI}\left[\begin{array}{c}0\\ 0\\ mg\end{array}\right]$$

(49)

where g is the acceleration due to gravity.

As the gravitational force acting on the central body produces no moment, the resultant gravitational moment is produced by the abdomen

$${\left[M_G\right]}^B=({\left[R\right]}^{BT}{\left[s_{tj}\right]}^T+s_{jb})\times {\left[R\right]}^{BI}\left[\begin{array}{c}0\\ 0\\ m^{T}g\end{array}\right]$$

(50)

The resultant aerodynamic forces and moments produced for fixed-wing aircraft are well established [96,127] and are represented as

$${\left[F_A\right]}^B=\frac{1}{2}\rho V^2{S}_{ref}\left[\begin{array}{c}{C}_x\\ C_y\\ C_z\end{array}\right]\mathrm{and}{\left[M_A\right]}^B=\frac{1}{2}\rho V^2{S}_{ref}\left[\begin{array}{c}{C}_l{b}_{ref}\\ C_m{c}_{ref}\\ C_m{b}_{ref}\end{array}\right]$$

(51)

where $\rho$ is the density of air, $c_{ref}$, $b_{ref}$ and $S_{ref}$ are the respective reference chord, span and area of the wing. $C_{(x,y,z,l,m,n)}$ are the dimensionless aerodynamic coefficients.

2.8. Control Effectors

For propulsion, the control input is the thrust $\left(T_n\right)$. The control inputs for abdominal motion are the abdominal mass angular deflections $({\varphi }_T,{\theta }_T,{\psi }_T)$, relative to the body frame. The aerodynamic control surfaces, which are the left and right elevons, denoted by ${\eta }_l$ and ${\eta }_r$ respectively, function as elevators or ailerons, depending on the desired flight regime. The combined deflections of the elevons as elevators ${\delta }_e$ or ailerons ${\delta }_a$ are according to [127]

$$\left(\begin{array}{c}{\delta }_e\\ {\delta }_a\end{array}\right)=\left(\begin{array}{cc}1& 1\\ -1& 1\end{array}\right)\left(\begin{array}{c}{\eta }_r\\ {\eta }_l\end{array}\right)$$

(52)

Following the sign convention for control effectors in flight dynamics [127], where applicable, a downward deflection is positive and a deflection to the right is positive, with negative to the left.

3. Control System Design

In this section, the control system is defined for a linearized system. Let the control objective be for the output $\mathbf{y}$ to track a constant reference signal $\mathbf{r}$, such that $\mathbf{y}=\mathbf{r}$ is reached asymptotically. The high-level block diagram of the control system is shown in Figure 3.

3.1. System Description

The nonlinear equations of motion have been derived in Section 2.6 and Section 2.7, which can be represented as a function of states $\mathbf{x}$ and inputs $\mathbf{u}$ in the form,

$$\begin{array}{c}\hfill \dot{\mathbf{x}}=f(\mathbf{x},\mathbf{u})\\ \hfill \mathbf{y}=g(\mathbf{x},\mathbf{u})\end{array}$$

(53)

where $\mathbf{x}\in {\mathbb{R}}^n$ is state of the system, $\mathbf{u}\in {\mathbb{R}}^m$ and $\mathbf{y}\in {\mathbb{R}}^p$ is the output of the system.

Linearizing the nonlinear model developed around a steady-state condition (trimmed flight) produces a linear system of the form

$$\begin{array}{c}\hfill \dot{\mathbf{x}}=\mathbf{Ax}+\mathbf{B}\mathbf{u}\\ \hfill \mathbf{y}=\mathbf{Cx}+\mathbf{Du}\end{array}$$

(54)

where matrix $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ is the system matrix, $\mathbf{B}\in {\mathbb{R}}^{n\times m}$ is the input matrix, $\mathbf{C}\in {\mathbb{R}}^{l\times n}$ is the output matrix, $\mathbf{D}\in {\mathbb{R}}^{l\times m}$ is the feedthrough matrix, and $\left(\mathbf{A},\mathbf{B}\right)$ can be stabilized. It is assumed that all the states $\mathbf{x}$ are measurable.

Controllability and Observability

In designing linear control systems, controllability and observability of the system should be investigated. Controllability determines the possibility of achieving the desired response of the system states $\mathbf{x}$ using the input $\mathbf{u}$. To check the controllability ($\mathcal{C}$) of a system, the matrix given in Equation (55) is used. If the matrix is full rank, the system is controllable [96].

$$\mathcal{C}=[\mathbf{B},\mathbf{AB},{\mathbf{A}}^2\mathbf{B},\cdots ,{\mathbf{A}}^{n-1}\mathbf{B}]$$

(55)

Observability ($\mathcal{O}$), on the other hand, is to know the possibility of determining all the states from the output and input signals. The system is observable if the matrix in Equation (56) is full rank [96].

$$\mathcal{O}=\left[\begin{array}{c}\mathbf{C}\\ \mathbf{CA}\\ ⋮\\ {\mathbf{CA}}^{n-1}\end{array}\right]$$

(56)

3.2. Optimal Linear Quadratic Regulator Control Theory

The control law for an optimal linear quadratic regulator (LQR) is presented in this section. For the system described in Equation (54), consider a cost function $\mathbf{J}$ defined as

$$\mathbf{J}={\int }_0^{\infty }({\mathbf{x}}^{\prime }\mathbf{Q}\mathbf{x}+{\mathbf{u}}^{\prime }\mathbf{R}\mathbf{u})dt$$

where $\mathbf{Q}={\mathbf{Q}}^{{}^{\prime }},\mathbf{Q}\in {\mathbb{R}}^{n\times n}$ is a positive semidefinite weighting matrix for the states and $\mathbf{R}={\mathbf{R}}^{{}^{\prime }},\mathbf{R}\in {\mathbf{R}}^{m\times m}$ is a positive definite weighting matrix for the control inputs. The optimal control input is given by

$$\begin{array}{cc}\hfill \mathbf{u}\left(t\right)& =-\mathbf{K}\mathbf{x}\left(t\right)\hfill \\ \hfill \mathrm{where}\mathrm{the}\mathrm{gain}\mathbf{K}& ={\mathbf{R}}^{-1}{\mathbf{B}}^{{}^{\prime }}\mathbf{P}\hfill \end{array}$$

(57)

$\mathbf{P}$ represents the unique positive semidefinite and symmetric solution to the Algebraic Riccati Equation (ARE)

$$\begin{array}{c}\hfill \mathbf{A}{}^{\prime }\mathbf{P}+\mathbf{P}\mathbf{A}+\mathbf{Q}-\mathbf{P}\mathbf{B}{\mathbf{R}}^{-1}\mathbf{B}{}^{\prime }\mathbf{P}=0\end{array}$$

(58)

If the pair $(\mathbf{A},\mathbf{B})$ is stabilizable and the pair $(\mathbf{A},\mathbf{Q})$ is detectable, then Equation (58) has a unique positive definite solution $\mathbf{P}$ such that the optimal control

$$\begin{array}{c}\hfill \mathbf{u}=-{\mathbf{R}}^{-1}{\mathbf{B}}^{{}^{\prime }}\mathbf{P}\mathbf{x}\end{array}$$

(59)

is asymptotically stabilizing [128,129].

3.3. Linear Quadratic Integral (LQI) Control

Since the aim of the controller is to track a non-zero reference input, integral control is also added to improve the tracking performance of the LQR controller, hence the term linear quadratic integral (LQI) control. LQI is, therefore, an augmented version of the LQR controller that designs an optimal controller, using full state feedback to optimize the quadratic cost function. Additional states are introduced into the system as the integrals of the error e, between the reference command and the actual output of the system. For the output equation in Equation (54), the LQI state feedback control law is given by

$$\mathbf{u}=-\mathbf{K}[\mathbf{x};{\mathbf{x}}_i]=-\mathbf{K}\mathbf{z}$$

(60)

where $\mathbf{z}$ are the states of the augmented system and ${\mathbf{x}}_i$ are the additional integral states. Upon introducing the integral states, the state–space representation of the augmented system becomes

$$\left[\begin{array}{c}\dot{\mathbf{x}}\\ {\dot{\mathbf{x}}}_i\end{array}\right]=\left[\begin{array}{cc}\mathbf{A}& 0\\ -\mathbf{C}& 0\end{array}\right]\left[\begin{array}{c}\dot{\mathbf{x}}\\ {\dot{\mathbf{x}}}_i\end{array}\right]+\left[\begin{array}{c}\mathbf{B}\\ -\mathbf{D}\end{array}\right]\mathbf{u}$$

(61)

Equation (61) above implies that

$$\left\{\begin{array}{cc}\dot{\mathbf{z}}\hfill & ={\widehat{\mathbf{A}}}_{aug}\mathbf{z}+{\widehat{\mathbf{B}}}_{aug}\mathbf{u}\hfill \\ \mathbf{y}\hfill & \mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u}\hfill \end{array}\right.$$

(62)

Hence, the general cost function for the augmented system is the given by

$$\mathbf{J}={\int }_0^{\infty }({\mathbf{z}}^{\prime }\mathbf{Q}\mathbf{z}+{\mathbf{u}}^{\prime }\mathbf{R}\mathbf{u})dt$$

where $\mathbf{z},\mathbf{u},\mathbf{Q}$, and $\mathbf{R}$ represent the augmented state vector, control vector, state weight matrix and input weight matrix, respectively. For the augmented system, $\mathbf{Q}$ and $\mathbf{R}$ are design parameters that can be adjusted to meet the design requirements of the controller.

4. Simulation and Results

The flight dynamics model of the conceptual DISWA was simulated in the MATLAB/Simulink environment [130].

4.1. Aircraft Specification

The DISWA model used in this study was obtained from a previous study [131]. The properties of the aircraft are summarized in

Table 1. Aircraft model physical properties [131].

Parameter	Value	Parameter	Value
$m^B$ (kg)	0.325	$m^T$ (kg)	0.06
Body length, $l_B$ (m)	0.3	Tail length, $l_{T_{max}}$ (m)	0.4
Max. body diameter, $d_{B_{max}}$ (m)	0.14	Tail diameter, $d_T$ (m)	0.05
${I_{xx}}_b^B$ (kg$·$${\mathrm{m}}^2$)	0.00187	$b_{ref}$ (m)	1.4
${I_{yy}}_b^B$ (kg$·$${\mathrm{m}}^2$)	0.01117	$c_{ref}$ (m)	0.19434
${I_{zz}}_b^B$ (kg$·$${\mathrm{m}}^2$)	0.00934	$S_{ref}$ (${\mathrm{m}}^2$)	0.26865
$c{g}_b^B$ (m)	[−0.064; 0; 0.003]	$ARP$ (m)	[0.025; 0; 0]

. The aerodynamic properties of the aircraft are assumed to be solely produced by the wings. The Athena Vortex Lattice (AVL) tool developed by Drela [132] was used to obtain the aerodynamic data for the DISWA model. AVL produces acceptable aerodynamic data that has been used for real aircraft applications [133]. See Appendix A for the AVL geometry file used in this study.

4.2. Model Validation: Single-Body vs. Multibody

The verification/validation of the proposed multibody model presented in Section 2.6 and Section 2.7 are presented in this section. The states and control inputs for the aircraft considered are:

$$\mathbf{x}={(u,v,w,p,q,r,\varphi ,\theta ,\psi ,X,Y,Z)}^T,$$

(63)

$$\mathbf{u}={({\delta }_e,{\delta }_a,T_n,{\varphi }_T,{\theta }_T,{\psi }_T)}^T,$$

(64)

The single-body and multibody model of the same aircraft should present the same longitudinal trim results. The translational and rotational flight dynamics equations of motion for a single rigid body not referenced to the center of mass have been derived by Bacon and Gregory in [134] as

$$\begin{array}{cc}\hfill & m\left({\left[{\dot{V}}_b^I\right]}^B+{\left[{\Omega }^{BI}\right]}^B{\left[V_b^I\right]}^B\right)+\left({\left[{\dot{\Omega }}^{BI}\right]}^{B}m{\left[s_{cb}\right]}^B\right)\hfill \\ \hfill & +\left({\left[{\Omega }^{BI}\right]}^B\left({\left[{\Omega }^{BI}\right]}^{B}m{\left[s_{cb}\right]}^B\right)\right)={\left[F\right]}^B\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill & \left(m{\left[S_{cb}\right]}^B{\left[{\dot{V}}_b^I\right]}^B\right)+\left(m{\left[{\Omega }^{BI}\right]}^B\left({\left[S_{cb}\right]}^B{\left[V_b^I\right]}^B\right)\right)\hfill \\ \hfill & \left(m{\left[V_b^I\right]}^B\left({\left[{\Omega }^{BI}\right]}^B{\left[S_{cb}\right]}^B\right)\right)+{\left[I_b^C\right]}^B{\left[{\dot{\omega }}^{BI}\right]}^B\hfill \\ \hfill & +\left({\left[{\Omega }^{BI}\right]}^B{\left[I_b^C\right]}^B{\left[{\omega }^{BI}\right]}^B\right)={\left[M\right]}^B+\left[{\left[S_{cb}\right]}^B{\left[F\right]}^B\right]\hfill \end{array}$$

(66)

The longitudinal trim results of the single-body model in [134] and the multibody model presented in this paper are calculated and compared. For the single-body model, the aircraft rotation of the abdomen was not considered during the estimation for longitudinal trim. The longitudinal trim using the single-body model can be described as follows: for a given airspeed V and a fixed abdominal pitch angle ${\theta }_T$, the objective is to estimate the appropriate longitudinal pitch angle $\theta$ and controls elevator ${\delta }_e$ and thrust $T_n$, to make the longitudinal rates, $\dot{u},\dot{w},\dot{q}$ and $\dot{h}$ approach zero.

In the case of the multibody model, the longitudinal trim can be described as follows: for a given airspeed V and a fixed abdominal pitch angle ${\theta }_T$, the objective is to estimate the appropriate longitudinal pitch angle $\theta$ and controls, elevator ${\delta }_e$ and thrust $T_n$, required to keep the abdomen pitch angle ${\theta }_T$ at a given value to make the longitudinal rates $\dot{u},\dot{w},\dot{q},\dot{h}$ and ${\left[{\dot{\omega }}^{TB}\right]}^T$ approach zero. Because the abdominal motion dynamics are included in the multibody model, the torque ${\tau }_{ay}$ required to keep the abdominal pitch inclination angle at the specified value can also be estimated.

The longitudinal trim problem is solved in MATLAB for an airspeed $V_0=10$ m/s, a height $H_0=100$ m and abdominal pitch inclination angles ${\theta }_{T_0}=0,-10$ and −30${}^{\circ }$. The results are presented in

Table 2. Single-body vs multibody trim results at steady cruise trim condition.

${\mathit{V}}_0$ (m/s)	${\mathit{H}}_0$ (m)	${\mathit{\theta }}_{{\mathit{T}}_0}$ (${}^{\circ }$)	Model Type	${\mathit{\theta }}_0$ (${}^{\circ }$)	${\mathit{\delta }}_{\mathit{e}0}$ (${}^{\circ }$)	${\mathit{T}}_{{\mathit{n}}_0}$ (N)	${\mathit{\tau }}_{\mathit{ay}0}$ (Nm)
10	100	0	Single-body	−11.35	0.228	0.679	-
			Multibody	−11.35	0.228	0.679	−0.235
		−10	Single-body	−10.66	0.0707	0.683	-
			Multibody	−10.66	0.0707	0.683	−0.232
		−30	Single-body	−6.02	−0.96	0.708	-
			Multibody	−6.02	−0.96	0.708	−0.204

.

4.3. Linearized Aircraft Model

Applying the Jacobian linearization approach or small disturbance theory [96], at a steady-state condition (trimmed flight), the linearized state–space model of the aircraft was obtained. Since the aim of the proposed work is to track a reference pitch angle and yaw angle, the linearized model was decoupled into longitudinal and lateral models. The same servomotor was used for abdominal pitch and yaw motion control.

The aircraft trim data corresponding to steady cruise flight condition was estimated for an airspeed of $V_0=10$ m/s and height $H_0=100$ m for an undeflected abdomen. For steady cruise, the aircraft flies with constant translational and angular velocity components in the body frame, $(\dot{u},\dot{v},\dot{w},\dot{p},\dot{q},\dot{r}=0)$. Furthermore, $(v,p,r,\varphi ,\psi =0)$. Only the elevator ${\delta }_e$ and thrust $T_n$ were used as control inputs to trim the aircraft. MATLAB’s trim algorithm that uses sequential quadratic programming was used to obtain the trim solution. The resulting trimmed steady cruise flight condition parameters are presented in

Table 3. Aircraft initial state and control inputs at steady cruise trim condition.

State (${\mathbf{x}}_0$)		Control Input (${\mathbf{u}}_0$)
$u_0$ = 10 (m/s)	$p_0$ = 0 (${}^{\circ }$/s)	${\delta }_{e0}$ =  0.228 (${}^{\circ }$)	${\varphi }_{T_0}$ = 0 (${}^{\circ }$)
$v_0$ = 0 (m/s)	$q_0$ = 0 (${}^{\circ }$/s)	${\delta }_{a0}$ = 0 (${}^{\circ }$)	${\theta }_{T_0}$ = 0 (${}^{\circ }$)
$w_0$ = −0.198 (m/s)	$r_0$ = 0 (${}^{\circ }$/s)	$T_{n_0}$ = 0.68 (N)	${\psi }_{T_0}$ = 0 (${}^{\circ }$)
$X_0$ = 0 (m)	${\varphi }_0$ = 0 (${}^{\circ }$)
$Y_0$ = 0 (m)	${\theta }_0$ = −1.13 (${}^{\circ }$)
$Z_0$ = −100 (m)	${\psi }_0$ = 0 (${}^{\circ }$)

. Simulation constraints on states and control inputs are given in Equation (67).

The nonlinear model was trimmed about steady cruise flight condition $(v,p,q,r,\varphi ,\psi =0)$ with Matlab’s trim algorithm that uses sequential quadratic programming.

$$\left\{\begin{array}{c}\left|\alpha \right|\le 20{}^{\circ },\left|\beta \right|\le 30{}^{\circ },\left|p\right|\le 30{}^{\circ }/s,\left|q\right|\le 30{}^{\circ }/s,\left|r\right|\le 30{}^{\circ }/s\hfill \\ {\delta }_e,{\delta }_a,\in [-20,20]{(}^{\circ }),{\theta }_T,{\psi }_T\in [-60,60]{(}^{\circ })\hfill \end{array}\right.$$

(67)

The linearization algorithm in MATLAB was also used to linearize the trimmed aircraft. The longitudinal part of the linearized model consists of the following states: ${\mathbf{x}}_{long}=\left[uwq\theta \right]{}^{\prime }$, and control inputs: ${\mathbf{u}}_{long}=\left[{\delta }_e{\theta }_T\right]{}^{\prime }$. The lateral–directional part consists of the following states: ${\mathbf{x}}_{lat}=\left[vpr\varphi \psi \right]{}^{\prime }$, and control inputs: ${\mathbf{u}}_{long}=\left[{\delta }_a{\varphi }_T\right]{}^{\prime }$. The longitudinal state and control matrices, as well as the lateral state and control matrices are given in Equations (68) and (69), respectively.

$$A_{long}=\left[\begin{array}{cccc}-0.39& 0.18& 0.13& -9.81\\ -3.16& -26.62& 6.82& 0.26\\ -7.80& -67.54& -8.74& 0.64\\ 0& 0& 1& 0\end{array}\right],B_{long}=\left[\begin{array}{cc}0& 0\\ -2.2& 1.65\\ -6.43& -190.03\\ 0& 0\end{array}\right]$$

(68)

$$A_{lat}=\left[\begin{array}{ccccc}-0.25& -0.32& -9.91& 12.22& 0\\ -14.09& -85.94& 12.47& 0& 0\\ 0.45& -0.014& -0.49& -23.33& 0\\ 0& 1& -0.02& 0& 0\\ 0& 0& 1& 0& 0\end{array}\right],B_{lat}=\left[\begin{array}{cc}-0.033& -7.11\\ -18.71& 0\\ 0.32& -137.27\\ 0& 0\\ 0& 0\end{array}\right]$$

(69)

4.4. Effect of Abdomen Mass on Longitudinal Stability

Static stability in aircraft design provides insight into the flying qualities of an aircraft. It is the initial tendency of an aircraft to return to its equilibrium state when disturbed [96]. Where the main focus is longitudinal motion, static stability theory can be applied to cruise conditions, as well as to quasi- steady maneuvers. A way to determine the static longitudinal stability of an aircraft is through the change in pitching moment with lift, $\frac{d{C}_m}{d{C}_L}$. This derivative must be negative for a stable aircraft. In addition, the pitching moment at zero angle of attack must be positive. The same condition applies using the change in pitching moment with angle of attack as an indicator in some cases, i.e., $\frac{d{C}_m}{d\alpha }<0$ is indicative of a stable aircraft [96,135]. Another measure of static stability is the static margin ($SM$), which will be used to analyze the static stability of the aircraft in consideration in this study due to its simplistic and intuitive nature. The $SM$ of an aircraft is given by [136]

$$SM=\frac{X_{NP}-X_{c{g}_m}}{c_{ref}},$$

(70)

where $X_{NP}$ is the neutral point ($NP$), and $X_{c{g}_m}$ is the center of gravity ($cg$) position of the whole aircraft along $x_B$ axis.

Generally, if $cg$ is ahead of the $NP$, then $SM$ is positive and the aircraft is stable. The larger the $SM$, the more stable and less maneuverable the aircraft is. However, if the $cg$ is behind the $NP$, resulting in a negative $SM$, the aircraft is unstable in pitch, making it difficult or impossible to fly manually [136]. The $SM$ value is usually chosen based on the desired performance of the aircraft in terms of handling qualities and usually ranges from $(10\%\le SM\le 5\%)$, for piloted aircraft. The neutral point of the DISWA was estimated as ($X_{NP}=$ 0.046 m) in AVL [137,138,139].

The eigenvalues of the state matrix of a linearized model also provide insight into the dynamic stability of the system. The effect of abdominal mass of the aircraft of interest was investigated. Five abdominal mass values were considered to be a percentage of central body mass of the aircraft, namely 6%, 11%, 16%, 21% and 26%, such that they resulted in negative and positive $SM$ values. Each time, the aircraft was trimmed and linearized about steady cruise flight condition for an airspeed of $V_0=10$ m/s and height $H_0=100$ m with an undeflected abdomen. The obtained eigenvalues of the various $\mathbf{A}$ matrices are summarized in

Table 4. Effect of abdominal mass as a percentage of central body mass on longitudinal stability.

Percentage of Central Body Mass (%)	Static Margin (%)	Eigenvalues	Dynamic Stability
6	37.3	−20.80 ± 28.21i	Stable
		−0.22 ± 0.95i
11	24.6	−18.81 ± 23.02i	Stable
		−0.20 ± 0.69
16	9.5	−18.00 ± 20.39	Stable
		−0.16 ± 0.41
21	−2.7	−17.52 ± 18.91i	Unstable
		−0.34
		0.06
26	−13.9	−17.24 ± 17.77i	Unstable
		−0.55
		0.385

.

4.5. Control Design

For control in this study, the abdomen has been constrained to 2 DoF, relative to the central body: the abdominal pitch and yaw motions. In this section, the LQI controllers for lateral and longitudinal control are presented. For the pitch or yaw angle controller, it is desired that the reference value of the pitch angle ${\theta }_{ref}$ or yaw angle ${\psi }_{ref}$ be tracked with an overshoot of less than 4%, a steady-state error of less than 1%, and a settling time of less than 4 s. Commanding the yaw angle instead of the roll angle achieves turning by generating a sideslip angle, which then generates a lateral force to turn the aircraft. From a practical standpoint, it is desirable that all states and control effectors remain within physical limits. Additionally, the aircraft was assumed to be cruising at a constant altitude and airspeed.

4.5.1. Longitudinal Control Simulation Study

The eigenvalues of the ${\mathbf{A}}_{long}$ matrix of the pre-augmented state–space model shown in Equation (68) indicate that the system is stable since all real parts of the eigenvalues are negative

$$eig\left({\mathbf{A}}_{long}\right)=\left[\begin{array}{c}-17.73\pm 19.58\mathrm{i}\\ -0.15\pm 0.24\mathrm{i}\end{array}\right]$$

(71)

For the longitudinal controller, the variable being controlled is the pitch angle; hence, (${\mathbf{y}}_{long}=e_{\theta }$), furthermore, it is assumed that changes in $\theta$ do not affect the speed of the aircraft. Therefore, the corresponding axial velocity u components are ignored. In addition, thrust is assumed to be constant, therefore, the effects of changes in thrust are ignored. The $z_B$ axis velocity w was replaced with the angle of attack $\alpha$ using the relationship ($\alpha =\frac{w}{u_0}$). The longitudinal augmented states considered are ${\mathbf{z}}_{long}=\left[\alpha q\theta e_{\theta }\right]{}^{\prime }$, and the longitudinal control inputs: ${\mathbf{u}}_{long}=\left[{\delta }_e{\theta }_T\right]{}^{\prime }$. The resulting augmented matrices are

$$\begin{array}{l}{\widehat{\mathbf{A}}}_{long}=\left[\begin{array}{cccc}-2.66& 0.68& 0.025& 0\\ -67.54& -8.74& 0.64& 0\\ 0& 1& 0& 0\\ 0& 0& -1& 0\end{array}\right], {\widehat{\mathbf{B}}}_{long}=\left[\begin{array}{cc}0& 0\\ -2.20& 1.65\\ -6.43& -190.03\\ 0& 0\end{array}\right]\\ {\widehat{\mathbf{C}}}_{long}=\left[\begin{array}{cccc}0& 0& 0& 1\end{array}\right], {\widehat{\mathbf{D}}}_{long}=\left[\begin{array}{cc}0& 0\end{array}\right]\end{array}$$

(72)

The obtained longitudinal linearized model was checked to be controllable and observable using MATLAB functions ctrb $({\widehat{\mathbf{A}}}_{long},\widehat{\mathbf{B}}{}_{long})$ and obsv $({\widehat{\mathbf{A}}}_{long}$, ${\widehat{\mathbf{C}}}_{long})$, respectively. The observability and controllability matrices were fully ranked $(rank\left({\mathcal{O}}_{long}\right)=rank\left({\mathcal{C}}_{long}\right)=4)$. Therefore, the longitudinal model is controllable and observable.

Three pitch control (PC) strategies were investigated based on the selection of weighting matrices ${\mathbf{R}}_{long}$ to reflect the following:

• PC 1 —Prioritized use of elevator for pitch angle control, interpreted as cheap elevator cost and expensive abdominal pitch cost.
• PC 2—Prioritized use of abdominal pitch for pitch angle control, interpreted as expensive elevator cost and cheap abdominal pitch cost.
• PC 3—Relatively evenly prioritized use of elevator and abdominal pitch for pitch angle control, interpreted as cheap elevator cost, cheap abdominal pitch cost.

The same state-weighting matrices ${\mathbf{Q}}_{long}=diag\left[464500\right]$, which were selected using Bryson’s rule as a guide [140], were used for all cases considered. Then, the optimal gains were calculated using MATLAB’s lqr (A,B,Q,R) function. A summary of the cases considered for pitch angle control are presented in

Table 5. Simulation parameters for pitch control (PC) cases.

Pitch Control (PC) Case	${\mathbf{R}}_{\mathit{long}}$	${\mathbf{K}}_{\mathit{long}}$
PC 1	diag [$2\times {10}^2$ $8\times {10}^7$]	$\left[\begin{array}{cccc}-87.68& 10.29& 93.47& -158.11\\ 0& 0& 0& 0\end{array}\right]$
PC 2	diag [$8\times {10}^7$ $2\times {10}^{-2}$]	$\left[\begin{array}{cccc}0& 0& 0& 0\\ 0.03& -17.30& -75.53& 158.06\end{array}\right]$
PC 3	diag [$2\times {10}^{-2}$ $2\times {10}^{-2}$]	$\left[\begin{array}{cccc}-12.78& -0.69& -1.98& 4.22\\ 0.69& -17.29& -75.38& 158.06\end{array}\right]$

. As the control objective was to track a reference pitch angle while keeping all states and control efforts within bounds, the performance characteristics for all cases are compared in terms of settling time, overshoot and steady-state error. The selected gains resulted in the closed-loop step response shown in Figure 4 for the states, while the control efforts are shown in Figure 5. The performance characteristics of the control strategies are summarized in

Table 6. Comparison of performance characteristics for the different pitch control (PC) cases.

Characteristic	PC 1	PC 2	PC 3
Settling time (s)	3.08	2.97	2.96
Overshoot (%)	3.84	3.82	3.81
Steady-state error (%)	0	0	0

.

4.5.2. Lateral–Directional Simulation Study

Estimating the eigenvalues of the ${\mathbf{A}}_{lat}$ matrix of the pre-augmented state–space model in Equation (69) reveals that the system is unstable due to the existence of a positive real part in the characteristic roots.

$$eig\left({\mathbf{A}}_{lat}\right)=\left[\begin{array}{c}1.02\pm 3.67\mathrm{i}\\ -86.05\\ -2.68\\ 0\end{array}\right]$$

(73)

For the yaw angle controller, the variable being controlled is the yaw angle; hence, (${\mathbf{y}}_{lat}=e_{\psi }$). The $y_B$ axis velocity v was replaced with the sideslip angle, $\beta$, using the relationship $\beta =\frac{v}{u_0}$. The lateral augmented states considered were ${\mathbf{z}}_{lat}=\left[\beta pr\varphi \psi e_{\psi }\right]{}^{\prime }$, and the control inputs: ${\mathbf{u}}_{lat}=\left[{\delta }_a{\psi }_T\right]{}^{\prime }$. The resulting augmented matrices are

$$\begin{array}{l}{\widehat{\mathbf{A}}}_{lat}=\left[\begin{array}{cccccc}-0.025& -0.032& -0.99& 1.22& 0& 0\\ -14.09& -85.94& 12.47& 0& 0& 0\\ 0.45& -0.014& -0.49& -23.33& 0& 0\\ 0& 1& -0.02& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& -1& 0\end{array}\right], {\widehat{\mathbf{B}}}_{lat}=\left[\begin{array}{cc}-0.032& -7.11\\ -18.71& 0\\ 0.32& -137.27\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right]\\ {\widehat{\mathbf{C}}}_{lat}=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 1\end{array}\right], {\widehat{\mathbf{D}}}_{lat}=\left[\begin{array}{cc}0& 0\end{array}\right]\end{array}$$

(74)

Checking the lateral–directional linearized model controllability and observability using MATLAB functions ctrb $({\widehat{\mathbf{A}}}_{lat},{\widehat{\mathbf{B}}}_{lat})$ and obsv $({\widehat{\mathbf{A}}}_{lat}$, ${\widehat{\mathbf{C}}}_{lat})$, respectively resulted in fully ranked matrices $(rank\left({\mathcal{O}}_{lat}\right)=rank\left({\mathcal{C}}_{lat}\right)=6)$. Therefore, the lateral–directional model is controllable and observable.

Again, three yaw control (YC) strategies were investigated based on the selection of weighting matrices ${\mathbf{R}}_{lat}$ to reflect the following:

• YC 1—Prioritized use of aileron for yaw angle control, interpreted as cheap aileron cost, expensive abdominal yaw cost.
• YC 2—Prioritized use of abdominal yaw for yaw angle control, interpreted as expensive aileron cost, cheap abdominal yaw cost.
• YC 3—Relatively evenly prioritized use of elevator and abdominal pitch for yaw angle control, interpreted as cheap aileron cost, cheap abdominal yaw cost.

The same state-weighing matrix ${\mathbf{Q}}_{lat}=diag\left[66633500\right]$ was used for all cases considered. The choice of the state-weighting matrix was also influenced by the need to minimize the sideslip when the aircraft is turning. The estimated optimal gains, as well as other parameters for the yaw controllers, are presented in

Table 7. Simulation parameters for yaw control (YC) cases.

Roll Control (YC) Case	${\mathit{R}}_{\mathit{lat}}$	${\mathit{K}}_{\mathit{lat}}$
YC 1	diag [18 $\times {10}^7$]	$\left[\begin{array}{cccccc}-0.68& -0.92& 6.30& -33.16& 18.01& -22.35\\ 0& 0& 0& 0& 0& 0\end{array}\right]$
YC 2	diag [$8\times {10}^7$ 2]	$\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ -0.86& -0.08& -1.72& 0.27& -8.68& 15.81\end{array}\right]$
YC 3	diag [10 10]	$\left[\begin{array}{cccccc}-0.31& -0.09& 0.01& -2.59& -0.03& -0.38\\ -0.06& -0.03& -0.8& -0.05& -3.55& 7.06\end{array}\right]$

. The state trajectory and control effort responses are shown in Figure 6 and Figure 7 respectively. Additionally, a comparison of performance characteristics is presented in

Table 8. Comparison of performance characteristics for the different yaw control (YC) cases.

Characteristic	YC 1	YC 2	YC 3
Settling time (s)	3.33	6.25	2.29
Overshoot (%)	5.80	0.89	2.95
Steady-state error (%)	0	2.00	0.04

.

5. Discussion

5.1. Model Validation: Single-Body vs. Multibody

As shown in Table 2, the trim results for the DISWA using the single-body model and multibody model are similar at various flight speeds. The only difference is the motor torque ${\tau }_{ay}$, acting on the abdomen which represents the amount of torque needed to keep the abdominal pitch inclination angle at the specified values and can be treated as a disturbance torque for control purposes.

5.2. Effect of Abdominal Mass on Longitudinal Stability

From the results presented in Table 4, the longitudinal stability of the aircraft decreases with increasing mass of the abdomen. This is associated with the static margin of the aircraft. A statically stable aircraft can be dynamically unstable; however, the reverse is not the case. The aircraft was stable up until 16%; however, as it became statically unstable from 21%, it also became dynamically unstable. Additionally, one can observe two complex conjugate pairs of eigenvalues for the dynamically stable system (6–16%). The highly damped pairs are the short-period modes, and the lightly damped modes represent the phugoid modes. For dynamically unstable systems in general, control systems must be designed that can stabilize the aircraft during flight.

5.3. Pitch Attitude Control

The longitudinal control input plots in Figure 5 show which longitudinal control effectors were prioritized for use in pitch-up control. For PC 1, the elevator was used exclusively, indicated by a downward elevator deflection in Figure 5a, and no abdominal deflection in Figure 5b. PC 2 used abdominal pitch deflection exclusively as seen in Figure 5a, and no elevator deflection in Figure 5b. A downward elevator and upward abdominal pitch deflection shown in Figure 5a,b, respectively, show that the elevator and abdomen were simultaneously used for control in PC 3.

The effects of the control effector choice on longitudinal states over time according to the controller used are shown in Figure 4. The angle of attack $\alpha$ time history in Figure 4a reveals that exclusive use of a positive elevator deflection (PC 1) caused an initial decrease in angle of attack as expected. Similar results have been obtained in [141,142] for tailless aircraft. The deflection of the elevator causes a change in the camber of the airfoil of the wing and consequently changes the lift coefficient. A positive elevator deflection for a reflex airfoil, as in the case of the aircraft wing in this study, decreases the wing lift and pitching moment, therefore resulting in a decrease in angle of attack. However, the exclusive deflection of the abdomen (PC 2) shifts the combined $cg$ forward, and consequently causes a change in lift and drag reference offset, as well as gravitational moments, therefore generating control moments. The resulting effect is the initial increase in angle of attack. In the case of PC 3, these effects are combined and reduced. Regarding the time histories of the pitch rate q (Figure 4b) for the three pitch controllers, an increase in pitch angle resulted in an increase in pitch rate, as expected.

The simulation results presented in Figure 4, Figure 5 and Table 6 showed that the PC 3 controller had better tracking performance for pitch angle than the PC 1 and PC 2 controllers, with a settling time of 2.96 s and the least overshoot of 3.81%. The PC 1 controller was the poorest performer of the three controllers with the highest overshoot of 3.84% and the longest settling time of 3.08 s. All three controllers had zero steady-state error and fall within the set performance requirements. The corresponding control efforts for elevator deflection ${\delta }_e$ and abdominal pitch angle are well within the physical bounds with maximum absolute deflections being 0.37${}^{\circ }$ and, 0.04${}^{\circ }$ respectively. The states also fall within the simulation constraints set for $\alpha$ and q, with the maximum absolute values of 5.00${}^{\circ }$ and 27.14${}^{\circ }$/s, respectively. For parallel comparison, the results indicate that a dragonfly using its abdomen in conjunction with its wings would obtain improved controllability and control performance compared to just using the wings.

5.4. Yaw Attitude Control

Figure 7 shows the lateral–directional control input plots for yaw angle control. YC 1 uses ailerons exclusively and zero abdominal yaw deflection, abdominal yaw deflection exclusively and zero aileron deflection for YC 2, and both aileron and abdominal yaw deflection for YC 3. The transient response for the lateral–directional dynamic response in terms of sideslip angle $\beta$, roll angle $\varphi$ and roll rate p are presented in Figure 6a–c. Please note that prior to implementation of a controller, the aircraft was laterally unstable (Equation (73), which is not unusual for tailless aircraft [141,143]. Similar behavior can be observed based on the effect of yaw control effector choice on the transient of these states. The exclusive use of the aileron (YC 1) produces the most stable response in $\beta$, $\varphi$ and p, while using the abdominal yaw alone for control (YC 2) caused the aircraft to be more unstable as a result of dynamic effects of the oscillatory movement of the abdomen as seen in Figure 7b. Simultaneously using both ailerons and abdominal yaw (YC 3), on the other hand, reduced the instability experienced with the exclusive use of abdominal yaw in YC 2. The main remark regarding the state evolution of roll angle (Figure 6b) and roll rate (Figure 6c) is the presence of a strong coupling between lateral and directional dynamics, caused by aileron and abdominal yaw deflection.

A yaw controller generates a sideslip to induce a turn and a significant problem associated with turning of a rudderless aircraft is adverse yaw resulting from the differential lift and drag on the left and right wing. The sideslip angle $\beta$ is a key performance indicator for aircraft turn coordination, i.e., with zero sideslip angle $\beta =0$. The time history for sideslip angle $\beta$ in Figure 6a shows that the exclusive use of ailerons (YC 1) and exclusive use of abdominal yaw (YC 2) resulted in negative sideslip. However, simultaneously using ailerons and abdominal yaw achieved better turn coordination over time. Regarding the time histories of the yaw rate r (Figure 6d), for the three yaw controllers, an increase in yaw angle resulted in an increase in yaw rate as expected.

From the simulation results presented in Figure 6, Figure 7 and Table 8, it can be observed that YC 3 has a better response in terms of settling time of 2.29 s compared to that of YC 1 and YC 2 with settling times of 3.33 s and 6.25 s, respectively. In terms of overshoot, YC 2 has the least overshoot 0.89%, compared to YC 1 and YC 3 with overshoots of 5.80% and 2.95%, respectively. YC 2 gives the highest steady-state error of 2%, with the steady-state errors of YC 1 and YC 3 being 0% and 0.04%, respectively. In addition, the corresponding control efforts for aileron deflection ${\delta }_a$ and abdominal yaw angle fall within the set limits with maximum absolute deflections being 0.056${}^{\circ }$ and, 0.013${}^{\circ }$ respectively. The states also fall within the simulation constraints set for $\beta$, p and r, with maximum absolute values of 22.24${}^{\circ }$, 6.53${}^{\circ }$/s and 26.56${}^{\circ }$/s, respectively. Of all the controllers, only YC 3 meets the controller performance requirements. For insects, the lack of a dedicated yaw-controlling surface similar to the aerodynamic rudder of many aircraft causes inherent instability on the lateral axis; however, “ruddering” movements of the abdomen have been observed [144]. The results attest to the potential use of a yaw-like abdominal airframe-based control mechanism that can be beneficial for the reduction of sideslip during turns (see Figure 6a) as sideslip is a critical control state of the entire aircraft.

The results obtained for yaw control are consistent with the observation in the analysis by Eturk in [36] that used internal moving mass for turn trim and trajectory control. With the internal moving mass, the aircraft was able to achieve trimmed turning flight with no sideslip, which was not possible using conventional aerodynamic control surfaces. Unlike that study, however, reaction torques are apparent, due to the size and mass of the appendage.

The application of an externally articulated mass specifically and solely for lateral control should be designed systematically due to the tendency of the articulated mass in this configuration to reduce the lateral–directional flight quality of a tailless aircraft.

6. Conclusions

A tensor-based multibody flight dynamics simulation model of a dragonfly-inspired straight-flying wing aircraft that considers abdominal articulation was presented in this study. The longitudinal and lateral flight altering effects of the mechanism have been demonstrated dynamically in simulation. An effective pitch and yaw angle tracking controller based on optimal control theory was implemented that was only possible with this new formulation. Allocation of control effort between the moving mass and aerodynamic actuators was achieved using modern control techniques. This control concept offers an alternative to conventional aerodynamic control surfaces for unmanned aerial vehicles. In addition, the results demonstrate the benefit of abdominal articulation of insects or tailless aircraft for longitudinal and lateral attitude control, especially for minimizing sideslip angle during turns.

The aerodynamic effects of the abdomen have been neglected in this study. Future work will extend this model to create a more comprehensive model that incorporates appendage aerodynamics. Although testing was done in the range in which the aircraft had adequate aerodynamic control, it is worth considering the operational range of extremely small aircraft. Furthermore, only one value of actuated mass was considered for control design. Simulation results show the actuated mass ratio to the mass of the aircraft has effects on the control sensitivity of the mechanism. Therefore, control sensitivity may be increased using a heavier mass. Finding an optimal mass for such a mechanism should be considered for future analysis. Additionally, the LQI control design procedure followed in this study can be extended for robustness to include uncertainties and actuator dynamics. Control allocation strategies can also be added to address issues such as actuator input saturation and rate limits, actuator redundancy, and fault tolerance of actuators and effectors.