Modeling, Simulation, and Cruise Characteristics of Wingtip-Jointed Composite Aircraft

Abstract

In this paper, multibody dynamic modeling and a simulation method for the wingtip-jointing process of a new-concept composite aircraft system are investigated. When the wingtips of two aircraft are jointed, the resultant wingtip-jointed aircraft is regarded as variable-geometry multiple rigid bodies, and a seven-degree-of-freedom non-linear dynamic model is established by mathematical derivation. The slip-meshing method is adopted to analyze the unsteady aerodynamic influence. We also present specific aerodynamic database acquisition methods under the quasi-steady assumption. Based on this, the simulation results indicate that the longitudinal and lateral movements are highly jointed and complex. A new composite aircraft system is investigated, in order to meet the balance requirement. With the lift–drag ratio (K) considered, the piecewise cubic Hermite interpolation (PCHIP) method, with a sufficient sample size, was utilized to help the cruise strategy optimization analysis under fixed altitude and speed conditions. Meanwhile, distribution of cruise characteristics with different sampling values of composite flight characteristic parameters were also analyzed. The research results can be used as a reference for new-concept composite aircraft model establishment, simulation, and multibody dynamic characteristic investigation.

1. Introduction

As aeronautical technologies develop, new-concept aircraft with better flight performance are created. Wingtip-jointed composite aircraft, representing a specific sort of new-concept aircraft, can form a composite aircraft system with wingtips jointed together when flying. As their aerodynamic performance is optimized by reduced induced drag force, such composite aircraft have the capacity for better range and endurance.

The concept of wingtip-docking technology was originally proposed by Dr. Richard Vogt, a German scientist, during World War II [1]. The United States Air Force (USAF) carried out the Tom-Tom project [2], the FICON project [3], and several other projects and experiments [4] to validate this theory, utilizing varieties of conceptual wingtip-docking models. Professor Edgar of Draper Laboratory and Cornell University demonstrated the aerodynamic modification created by the combination of two separate wings using the vortex lattice method [5], and wind tunnel testing [6]. The FICON project conducted by the USAF designed a strategic bombing unit named a flying aircraft carrier with a B-29 bomber carrying two F-84 D/E fighters. After the verification of wind tunnel testing, ground effect experiments and several successful hookups were performed; one fighter was torn away from the right wingtip of the other due to insufficient collaboration among the aircrafts [7]. The F-84 immediately flipped over onto the wing of the B-29 and both crashed, with the loss of all onboard [8].

Since then, relevant research into this technology has focused on the dynamics of multibody systems with rotational hinge constraints. A multibody system refers to a mechanical system composed of multiple bodies interconnected in a certain way. Multibody dynamics is a new branch of the discipline based on classical mechanics.

To deal with a multibody system, in principle, traditional classical mechanics methods are typically used, including the vector mechanics method represented by the Newton–Euler equation and the analytical mechanics method represented by the Lagrange principle. Based on this, the Roberson–Wittenburg method mainly uses graph theory and mathematical tools to describe the structure of a multibody system, using the relative displacement between adjacent rigid bodies as generalized co-ordinates to derive universal dynamics suitable for arbitrary multi-rigid body systems. Fisher’s augmented body concept was used to create a physical interpretation of the coefficient matrix A of the equation [9].

Bryan Cannon et al. analyzed the modeling and control method of an aeroelastic morphing UAV [10]. Animesh Chakravarthy of Florida University analyzed the modeling method of a specific variable-sweep aircraft [11]. Yue Ting and Wang Lixin have conducted dynamic characteristics analyses and flight control designs for Z-wing morphing aircraft [12,13], oblique wing aircraft [14,15], and variable-sweep morphing vehicles [16]. Seigler et al. [17,18], Valasek et al. [19], S. Li [20], and B. Yan [21] have developed different dynamic models for other similar multibody aircraft systems.

For a wingtip-jointed composite aircraft system, Brandon Troub and Carlos Montalvo investigated a traditional Newton–Euler dynamic modeling method which adds constraints to a single aircraft. At the same time, the gain scheduling method was used to complete the control characteristics study, and the composite trajectory was discussed [22]. However, for a composite system composed of multiple aircraft based on the aerodynamic data and trim strategy of a single aircraft, a multi-aircraft composite method giving a multi-aircraft trim strategy was not established. The existing research studies mainly finish the composite aircraft study based on the aerodynamic data and trim strategy of a single aircraft, while the dynamic model establishment, aerodynamic database, and composite trimming strategy simulation of a composite system will help the aircraft cruise simulation results be more reasonable.

When a multibody system consisting of two aircraft is dynamically analyzed, the force and the movement (especially for the relative roll motion) of the composite aircraft are highly coupled between the two composite aircrafts. A non-linear multibody dynamic model should be established with the force and movement on each degree of freedom (DOF) analyzed. To guarantee flying stability, a new trimming strategy for the composite system may have to be arranged, where simulation can be conducted to obtain the associated dynamic characteristics.

To present our study regarding the above problem, this paper is organized as follows: Section 2 presents the derivation of the full-parameter dynamic equations and kinematic equations of the wingtip-jointed composite system, including the establishment of aerodynamic force and momentum models under the quasi-steady assumption. Section 3 presents the simulation method and analysis of the intrinsic dynamics responses with two different trimming strategies, in order to judge the simulation results of the strategy, which can help the composite system to meet the cruise requirement. With the lift-to-drag ratio mainly considered, the optimal cruise strategy under fixed altitude and speed is investigated using the PCHIP method, as detailed in Section 4, where optimized values of parameters among all the samples are correspondingly presented.

2. Dynamic Modeling of Wingtip-Jointed Composite Aircraft

2.1. Wingtip-Jointed Composite Aircraft

Composite aircraft consist of two fixed-wing aircraft with a column hinge mechanism emplaced separately, which can permit relative roll movement as well as providing relative translation movement, relative pitch movement, and relative yaw movement constraints, as shown in Figure 1.

The wingtip hinge mechanism is mainly composed of the hinge column restraint mechanism, as shown in Figure 2. This mechanism restricts three relative translational degrees of freedom (DOFs) and two relative rotational DOFs, including the relative translational motion, pitch relative rotation, and yaw relative rotation in the fore–after, up–down, and left–right directions, where a relative roll motion between the two aircraft is allowed.

Multibody system dynamics is a branch of science based on classical mechanics. A system of rigid bodies is defined as a multibody system of multiple bodies. Research on wingtip-jointed composite aircraft concerns the dynamic problem of a rootless tree system with rotational hinge constraints. The number of equations can be determined according to the degrees of freedom (DOFs) and constraint forms.

In general, a conventional single aircraft can be regarded as one rigid body that has 6 DOFs, including 3 translational DOFs and 3 rotational DOFs, so the multibody system has 12 DOFs originally, while the internal hinge joint structure regulates 5 relative DOFs. Therefore, as shown in Figure 3, this multibody system has 7 DOFs, including 3 translational DOFs, 1 pitch rotational DOF, 1 yaw rotational DOF, and 2 separate rotational DOFs. The dynamic and kinematic equations to be derived contain overall equations and inner equations of the composite system.

2.2. Derivation of Dynamic and Kinematic Equations of the Composite System

In the flying process of wingtip-jointed composite aircraft, the relative movement between each aircraft can cause a shift of the mass center of the multibody system. As shown in Figure 4, composite body co-ordinate system ($o_{fb}{x}_{fb}{y}_{fb}{z}_{fb}$) was chosen, where the origin of this co-ordinate system is located at the wingtip joint. A composite air co-ordinate system ($O{x}_{fa}{y}_{fa}{z}_{fa}$) and path coordinate system ($O{x}_{fk}{y}_{fk}{z}_{fk}$) are defined in the corresponding method. The ground co-ordinate system is described as $o_g{x}_g{y}_g{z}_g$.

According to multibody dynamics theory, the velocities of the two aircraft in the multibody system can be described as:

$$\overrightarrow{v_i}=\overrightarrow{v_o}+{\dot{\overrightarrow{r}}}_{OCi}+{\overrightarrow{\omega }}_i\times {\overrightarrow{r}}_{OCi}, i=1,2.$$

(1)

The momentum $\overrightarrow{p}$ and the moment of momentum $\overrightarrow{L}$ can be described as:

$$\begin{array}{l}\overrightarrow{p}={\displaystyle \sum _i^n{m}_i\left({\overrightarrow{v}}_i+{\dot{\overrightarrow{r}}}_{OCi}+{\overrightarrow{\omega }}_i\times {\overrightarrow{r}}_{OCi}\right)}, i=1,2.\\ \overrightarrow{L}={\displaystyle \sum _i^n\left({\overrightarrow{r}}_{OCi}\times m{\overrightarrow{v}}_i+{\overleftrightarrow{I}}_{Ci}\cdot \left(\overrightarrow{\omega }+{\overrightarrow{\omega }}_{ir}\right)\right)}, i=1,2.\end{array}$$

(2)

The augmented body method [13] introduced the definition of the augmented body mass M, the augmented body rotation tensor $\overleftrightarrow{I_c}$, and the augment body static moment ${\stackrel{\rightharpoonup }{S}}^{*}$. The relative equations are defined, respectively, as follows:

$$M=\sum _{i=1}^k\left(m_i\right),$$

(3)

$${\overleftrightarrow{I}}_C^{*}={\displaystyle {\iiint }_V\rho \left(r^2\overleftrightarrow{1}-\overrightarrow{r}\overrightarrow{r}\right)}dV=\left[\begin{array}{ccc}{I}_x^{*}& -I_{xy}^{*}& -I_{zx}^{*}\\ -I_{xy}^{*}& I_y^{*}& -I_{yz}^{*}\\ -I_{zx}^{*}& -I_{yz}^{*}& I_z^{*}\end{array}\right],$$

(4)

$${\overrightarrow{S}}^{*}={\displaystyle \int \overrightarrow{r}\times dm}=S_x\overrightarrow{i}+S_y\overrightarrow{j}+S_z\overrightarrow{k}.$$

(5)

The numbers of augmented body parameters in different cases are shown in

Table 1. Augmented body parameters in different cases.

		Case	Single Aircraft	Double Composite Aircraft
	Number
Parameters				${\mathit{\phi }}_{12}=-20°$	${\mathit{\phi }}_{12}=-10°$	${\mathit{\phi }}_{12}=0°$	${\mathit{\phi }}_{12}=10°$	${\mathit{\phi }}_{12}$ = 20°
$m$				1500	3000
$I_x$				1420	72,840
$I_y$			4060	10,920	9680	8130	9680	10,920
$I_z$			478	11,780	10,840	9570	10,840	11,780
$\left|\overrightarrow{S}\right|$			0	2390	1210	0	−1210	−2390

. The variable double composite aircraft numbers change obviously, compared with those for the single aircraft case. The relative roll angle, ${\phi }_{12}$, can determine the corresponding numerical value of the augmented body parameters.

Based on Equations (2)–(5), the equations of $\overrightarrow{p}$ and $\overrightarrow{L}$ become

$$\begin{array}{l}\overrightarrow{p}=M\cdot \overrightarrow{v}+{\dot{\overrightarrow{S}}}^{*}+\overrightarrow{\omega }\times {\overrightarrow{S}}^{*}\\ \overrightarrow{L}={\overrightarrow{S}}^{*}\times \overrightarrow{v}+{\overleftrightarrow{I}}_O^{*}\cdot \overrightarrow{\omega }\end{array}$$

(6)

The origin of the composite body co-ordinate system is not always located at the mass center of the whole multibody system, according to the momentum theorem:

$$\begin{array}{l}\overrightarrow{F}=\dot{\overrightarrow{p}}\\ \overrightarrow{M}=\dot{\overrightarrow{L}}+\overrightarrow{r}\times m\dot{\overrightarrow{v}}=\dot{\overrightarrow{L}}+\overrightarrow{v}\times \dot{\overrightarrow{S}}\end{array}$$

(7)

Equation (7) can be further derived as:

$$\{\begin{array}{l}\overrightarrow{F}=m(\dot{\overrightarrow{v}}+\overrightarrow{\omega }\times \overrightarrow{v})+\frac{{\delta }^2\overrightarrow{S}}{\delta t^2}+2\overrightarrow{\omega }\times \frac{\delta \overrightarrow{S}}{\delta t}+\frac{\delta \overrightarrow{\omega }}{\delta t}\times \overrightarrow{S}+\overrightarrow{\omega }\times (\overrightarrow{\omega }\times \overrightarrow{S})\\ \overrightarrow{M}=\overrightarrow{S}\times \frac{\delta \overrightarrow{v}}{\delta t}+\overrightarrow{\omega }\times \left(\overrightarrow{S}\times \overrightarrow{v}\right)+{\overleftrightarrow{I}}_O^{}\cdot \frac{\delta \overrightarrow{\omega }}{\delta t}+\frac{\delta {\overleftrightarrow{I}}_O^{}}{\delta t}\cdot \overrightarrow{\omega }+\overrightarrow{\omega }\times ({\overleftrightarrow{I}}_O^{}\cdot \overrightarrow{\omega })\end{array}.$$

(8)

In Equation (8), the augmented body rotational tensor and static moment are associated with the connection form and motion of multiple rigid body system, describing the shift of the mass center of the whole multibody system. According to the augmented body method, the scalar form of the non-linear dynamic equation is derived as:

$$\begin{array}{l}\{\begin{array}{l}{F}_x=m({\dot{v}}_x+{\omega }_y{v}_z-{\omega }_z{v}_y)+2{\omega }_y{\dot{S}}_z-2{\omega }_z{\dot{S}}_y+{\dot{\omega }}_y{S}_z-{\dot{\omega }}_z{S}_y+{\omega }_y\left({\omega }_x{S}_y-{\omega }_y{S}_x\right)-{\omega }_z\left({\omega }_z{S}_x-{\omega }_x{S}_z\right)+{\ddot{S}}_x\\ F_y=m({\dot{v}}_y+{\omega }_z{v}_x-{\omega }_x{v}_z)+2{\omega }_z{\dot{S}}_x-2{\omega }_x{\dot{S}}_z+{\dot{\omega }}_z{S}_x-{\dot{\omega }}_x{S}_z+{\omega }_z\left({\omega }_y{S}_z-{\omega }_z{S}_y\right)-{\omega }_x\left({\omega }_x{S}_y-{\omega }_y{S}_x\right)+{\ddot{S}}_y\\ F_z=m({\dot{v}}_z+{\omega }_x{v}_y-{\omega }_y{v}_x)+2{\omega }_x{\dot{S}}_y-2{\omega }_y{\dot{S}}_x+{\dot{\omega }}_x{S}_y-{\dot{\omega }}_y{S}_x+{\omega }_x\left({\omega }_z{S}_x-{\omega }_x{S}_z\right)-{\omega }_y\left({\omega }_y{S}_z-{\omega }_z{S}_y\right)+{\ddot{S}}_z\end{array}\\ \{\begin{array}{l}{M}_x=S_y{\dot{v}}_z-S_z{\dot{v}}_y+{\omega }_y\left(S_x{v}_y-S_y{v}_x\right)-{\omega }_z\left(S_z{v}_x-S_x{v}_z\right)+{\dot{I}}_x{\omega }_x-{\dot{I}}_{xy}{\omega }_y-{\dot{I}}_{zx}{\omega }_z+I_x{\omega }_x-I_{xy}{\omega }_y-I_{zx}{\omega }_z+{\omega }_y\left(-I_{zx}{\omega }_x-I_{yz}{\omega }_y+I_z{\omega }_z\right)-{\omega }_z\left(-I_{xy}{\omega }_x+I_y{\omega }_y-I_{yz}{\omega }_z\right)\\ M_y=S_z{\dot{v}}_x-S_x{\dot{v}}_z+{\omega }_z\left(S_y{v}_z-S_z{v}_y\right)-{\omega }_x\left(S_x{v}_y-S_y{v}_x\right)-{\dot{I}}_{xy}{\omega }_x+{\dot{I}}_y{\omega }_y-{\dot{I}}_{yz}{\omega }_z-I_{xy}{\omega }_x+I_y{\omega }_y-I_{yz}{\omega }_z+{\omega }_z\left(I_x{\omega }_x-I_{xy}{\omega }_y-I_{zx}{\omega }_z\right)-{\omega }_x\left(-I_{zx}{\omega }_x-I_{yz}{\omega }_y+I_z{\omega }_z\right)\\ M_z=S_x{\dot{v}}_y-S_y{\dot{v}}_x+{\omega }_x\left(S_z{v}_x-S_x{v}_z\right)-{\omega }_y\left(S_y{v}_z-S_z{v}_y\right)-{\dot{I}}_{zx}{\omega }_x-{\dot{I}}_{yz}{\omega }_y+{\dot{I}}_z{\omega }_z-I_{zx}{\omega }_x-I_{yz}{\omega }_y+I_z{\omega }_z+{\omega }_x\left(-I_{xy}{\omega }_x+I_y{\omega }_y-I_{yz}{\omega }_z\right)+{\omega }_y\left(I_x{\omega }_x-I_{xy}{\omega }_y-I_{zx}{\omega }_z\right)\end{array}\end{array}$$

(9)

Meanwhile, the overall kinematic equations of the composite system can be derived as

$$\begin{gathered} \left[\begin{array}{c}{\dot{x}}_{\mathrm{e}}\\ {\dot{y}}_e\\ {\dot{z}}_e\end{array}\right]=L_{ef}(\phi ,\theta ,\psi )\left[\begin{array}{c}u\\ v\\ w\end{array}\right], \\ \begin{array}{c}\dot{\phi }=p+q\sin \phi \tan \theta +r\cos \phi \tan \theta \\ \dot{\theta }=q\cos \phi -r\sin \phi \\ \dot{\psi }=q\sin \phi \sec \theta +r\cos \phi \sec \theta \end{array}. \end{gathered}$$

(10)

In the composite flight process, internal relative rotational movement exists on the additional 7th DOF, compared with a regular single aircraft’s flight process. As shown in Figure 5, ${\phi }_{12}$ is defined as the relative roll angle between the aircraft on the right side (subscript 1) and the one on the left side (subscript 2).

According to the moment of momentum theorem, the relationship between the motion and moment of momentum of the composite system can be expressed as:

$${\dot{\phi }}_{12}=-\frac{l_1}{I_{O_{f}1}}+\frac{l_2}{I_{O_{f}2}},$$

(11)

which can be derived as:

$$I_{O_{f}1}{I}_{O_{f}2}{\dot{\phi }}_{12}=l_2{I}_{O_{f}1}-l_1{I}_{O_{f}2}.$$

(12)

The geometric movement state and motion relationship can be derived as:

$$(\begin{array}{c}{\phi }_{12}={\phi }_2-{\phi }_1\\ V_{x2}=V_{x1}\\ V_{y2}=V_{y1}-{\dot{\phi }}_{12}l\sin {\phi }_{12}\\ V_{z2}=V_{z1}-{\dot{\phi }}_{12}l\cos {\phi }_{12}\end{array}$$

(13)

Based on the second-order Lagrange equation

$$\frac{d}{dt}[\frac{\partial L}{\partial {\dot{\phi }}_{12}}]-\frac{\partial L}{\partial {\phi }_{12}}=0,$$

(14)

where the Lagrange function is described as

$$L={\displaystyle \sum _{k=1}^n\{\frac{1}{2}{m}_k(V_{kx}^2+V_{ky}^2+V_{kz}^2)+}\frac{1}{2}{I}_{kx}{\phi }_k^2+mg{z}_k\}.$$

(15)

Thus, Equation (14) can be expanded as:

$${\dot{\phi }}_{12}=\frac{I_x\dot{p}+(m{l}^2+\frac{1}{2}{I}_x){\dot{p}}_{12}+ml({\ddot{z}}_{1e}\cos {\phi }_{12}-{\ddot{y}}_{1e}\sin {\phi }_{12})+\frac{1}{2}mgl}{ml({\dot{z}}_{1e}\sin {\phi }_{12}+{\dot{y}}_{1e}\cos {\phi }_{12})},$$

(16)

where $Z_{1e}$ and $y_{1e}$ represent the center of gravity of the aircraft on the right side. The zero-order variables in the composite body co-ordinate system are ($x_c,y_c,z_c,{\phi }_c,{\theta }_c,{\psi }_{12},{\phi }_{12}$). The final forms for the multibody dynamic equation and the kinematic equation are derived as Equations (17) and (18), respectively:

$$\begin{array}{l}(\begin{array}{l}X-mg\sin \theta +T\cos {\phi }_T=m({\dot{v}}_x+{\omega }_y{v}_z-{\omega }_z{v}_y)+2{\omega }_y{\dot{S}}_z-2{\omega }_z{\dot{S}}_y+{\dot{\omega }}_y{S}_z-{\dot{\omega }}_z{S}_y+{\omega }_y\left({\omega }_x{S}_y-{\omega }_y{S}_x\right)-{\omega }_z\left({\omega }_z{S}_x-{\omega }_x{S}_z\right)+{\ddot{S}}_x\\ Y+mg\cos \theta \sin \phi =m({\dot{v}}_y+{\omega }_z{v}_x-{\omega }_x{v}_z)+2{\omega }_z{\dot{S}}_x-2{\omega }_x{\dot{S}}_z+{\dot{\omega }}_z{S}_x-{\dot{\omega }}_x{S}_z+{\omega }_z\left({\omega }_y{S}_z-{\omega }_z{S}_y\right)-{\omega }_x\left({\omega }_x{S}_y-{\omega }_y{S}_x\right)+{\ddot{S}}_y\\ Z+mg\cos \theta \cos \phi -T\sin {\phi }_T=m({\dot{v}}_z+{\omega }_x{v}_y-{\omega }_y{v}_x)+2{\omega }_x{\dot{S}}_y-2{\omega }_y{\dot{S}}_x+{\dot{\omega }}_x{S}_y-{\dot{\omega }}_y{S}_x+{\omega }_x\left({\omega }_z{S}_x-{\omega }_x{S}_z\right)-{\omega }_y\left({\omega }_y{S}_z-{\omega }_z{S}_y\right)+{\ddot{S}}_z\end{array}\\ \{\begin{array}{l}l=S_y{\dot{v}}_z-S_z{\dot{v}}_y+{\omega }_y\left(S_x{v}_y-S_y{v}_x\right)-{\omega }_z\left(S_z{v}_x-S_x{v}_z\right)+{\dot{I}}_x{\omega }_x-{\dot{I}}_{xy}{\omega }_y-{\dot{I}}_{zx}{\omega }_z+I_x{\dot{\omega }}_x-I_{xy}{\dot{\omega }}_y-I_{zx}{\dot{\omega }}_z+{\omega }_y\left(-I_{zx}{\omega }_x-I_{yz}{\omega }_y+I_z{\omega }_z\right)-{\omega }_z\left(-I_{xy}{\omega }_x+I_y{\omega }_y-I_{yz}{\omega }_z\right)\\ m=S_z{\dot{v}}_x-S_x{\dot{v}}_z+{\omega }_z\left(S_y{v}_z-S_z{v}_y\right)-{\omega }_x\left(S_x{v}_y-S_y{v}_x\right)-{\dot{I}}_{xy}{\omega }_x+{\dot{I}}_y{\omega }_y-{\dot{I}}_{yz}{\omega }_z-I_{xy}{\dot{\omega }}_x+I_y{\dot{\omega }}_y-I_{yz}{\dot{\omega }}_z+{\omega }_z\left(I_x{\omega }_x-I_{xy}{\omega }_y-I_{zx}{\omega }_z\right)-{\omega }_x\left(-I_{zx}{\omega }_x-I_{yz}{\omega }_y+I_z{\omega }_z\right)\\ n=S_x{\dot{v}}_y-S_y{\dot{v}}_x+{\omega }_x\left(S_z{v}_x-S_x{v}_z\right)-{\omega }_y\left(S_y{v}_z-S_z{v}_y\right)-{\dot{I}}_{zx}{\omega }_x-{\dot{I}}_{yz}{\omega }_y+{\dot{I}}_z{\omega }_z-I_{zx}{\dot{\omega }}_x-I_{yz}{\dot{\omega }}_y+I_z{\dot{\omega }}_z+{\omega }_x\left(-I_{xy}{\omega }_x+I_y{\omega }_y-I_{yz}{\omega }_z\right)+{\omega }_y\left(I_x{\omega }_x-I_{xy}{\omega }_y-I_{zx}{\omega }_z\right)\end{array}\\ I_{O_{f}1}{I}_{O_{f}2}{\dot{\phi }}_{12}=l_2{I}_{O_{f}1}-l_1{I}_{O_{f}2}\end{array}$$

(17)

$$\begin{gathered} \left[\begin{array}{c}{\dot{x}}_{\mathrm{e}}\\ {\dot{y}}_e\\ {\dot{z}}_e\end{array}\right]=L_{ef}(\phi ,\theta ,\psi )\left[\begin{array}{c}u\\ v\\ w\end{array}\right], \\ \begin{array}{c}\dot{\phi }=p+q\sin \phi \tan \theta +r\cos \phi \tan \theta \\ \dot{\theta }=q\cos \phi -r\sin \phi \\ \dot{\psi }=q\sin \phi \sec \theta +r\cos \phi \sec \theta \end{array}, \\ {\dot{\phi }}_{12}=\frac{I_x\dot{p}+(m{l}^2+\frac{1}{2}{I}_x){\dot{p}}_{12}+ml({\ddot{z}}_{1e}\cos {\phi }_{12}-{\ddot{y}}_{1e}\sin {\phi }_{12})+\frac{1}{2}mgl}{ml({\dot{z}}_{1e}\sin {\phi }_{12}+{\dot{y}}_{1e}\cos {\phi }_{12})}. \end{gathered}$$

(18)

Equations (17) and (18) show that there are additional terms in the multibody dynamic equations, caused by the mass center shift in the composite flight process. Thus, there is a significant difference in the non-linear equations between composite aircraft with jointed wingtips and single conventional configuration aircraft. In Equation (17), $S_x$, $S_y,$ and $S_z$ represent the additional forces produced by mass center shift in composite flight, while $\dot{S_z}$ and $\ddot{S_z}$ represent the additional forces produced by the velocity and acceleration, respectively, of the aircraft’s mass center shift.

In Equation (17), the aerodynamic variables are computed as follows: If the relative roll angle of the wings can be fixed constantly—which means that ${\overleftrightarrow{I}}_C$ and ${\stackrel{\rightharpoonup }{S}}^{*}$ are constant—and their all-order derivatives are 0, then Equations (17) and (18) become similar to the six DOF non-linear equations for a single (conventional configuration) aircraft.

2.3. Aerodynamic Database Establishment

The relative rolling movement of wings and fuselages produces an unsteady aerodynamic force in the aircraft system, especially lateral and internal forces. In general, the unsteady aerodynamic influence on the drag force, lift force, and pitch moment of a conventional rigid-body aircraft can be ignored. In this study, the aerodynamic drag force, lift force, and pitch moment were obtained under the quasi-steady assumption [12], while the unsteady aerodynamic effect was not considered; the drag force, lift force, and pitch moment of the wingtip-jointed composite aircraft under a certain configuration were nearly equal to those of the corresponding static configuration. The aerodynamic force and moment of the wingtip-jointed composite aircraft system can be expressed as:

$$\{\begin{array}{l}D=qS(C_{D_0}+C_{D\alpha }\alpha +C_{DV}\cdot \Delta V/V+C_{D{\phi }_{12}}{\phi }_{12}+C_{D\delta }\cdot \delta \\ L=qS(C_{L_0}+C_{L\alpha }\alpha +C_{LV}\cdot \Delta V/V+C_{Lq}q/(2V/\overline{c})+C_{L\dot{\alpha }}\dot{\alpha }/(2V/\overline{c})+C_{L{\phi }_{12}}{\phi }_{12}+C_{L\delta }\cdot \delta \\ \mathrm{C}=qSb{(\mathrm{C}}_{C\beta }\beta +{\mathrm{C}}_{Cp}p/2V+C_{Cr}r/2V+C_{C{p}_{12}}{p}_{12}+C_{C{\phi }_{12}}{\phi }_{12}+C_{C\delta }\cdot \delta )\\ l=qSb\left({\mathrm{C}}_{l\beta }\beta +{\mathrm{C}}_{lp}pb/2V+{\mathrm{C}}_{lr}rb/2V+C_{l{\phi }_{12}}{\phi }_{12}+{\mathrm{C}}_{l{p}_{12}}{p}_{12}b/2V+C_{l\delta }\cdot \delta \right)\\ m=qS\overline{c}(C_{m0}+C_{m\alpha }\alpha +C_{mV}\cdot \Delta V/V+C_{mq}q/(2V/\overline{c})+C_{m\dot{\alpha }}\dot{\alpha }/(2V/\overline{c})+C_{m{\phi }_{12}}{\phi }_{12}+C_{m\delta }\cdot \delta )\\ n=qSb\left({\mathrm{C}}_{n\beta }\beta +{\mathrm{C}}_{np}pb/2V+{\mathrm{C}}_{nr}rb/2V+C_{n{p}_{12}}{p}_{12}+C_{n{\phi }_{12}}{\phi }_{12}+C_{n\delta }\cdot \delta \right)\\ l_1=-l_2=qSb\left(C_{l{1}_0}+{\mathrm{C}}_{l1\beta }\beta +C_{l1{\phi }_{12}}{\phi }_{12}+{\mathrm{C}}_{l1p_{12}}{p}_{12}b/2V+C_{l1\delta }\cdot \delta \right)\end{array},$$

(19)

where $\delta$ is a six-dimensional array consisting of the control surface deflection, which can be expressed as $\delta =\left({\delta }_{a_1},{\delta }_{a_2},{\delta }_{e_1},{\delta }_{e_2},{\delta }_{r_1},{\delta }_{r_2}\right)$.

Establishment of the aerodynamic database included the acquisition of aerodynamic derivatives of relative roll motion, aerodynamic derivatives of composite control surface deflection, and general aerodynamic parameters. In Figure 6, the field of calculus is a cylindrical flow field with 10 times forward fuselage length, 20 times body length backward, and 10 times span diameter. The Y axis points back along the axis of the fuselage. The unstructured grid of the wingtip-jointed aircraft system is shown in Figure 7, while detailed symbols and the corresponding data acquisition methods are listed in

Table 2. Aerodynamic database acquisition method.

Aerodynamic Parameters	Symbol	Data Acquisition Method
Aerodynamic derivatives of relative roll motion	$C_{D\phi 12}$, $C_{L\phi 12}$, $C_{m\phi 12}$, etc.	Solver: Pressure-based Turbulence model: Spalart–Allmaras Solution algorithm: simple Spatial discretization flow formation: second order upwind method
Aerodynamic derivatives of composite control surface deflection	$C_{D\delta }$, $C_{L\delta }$, $C_{m\delta }$, $C_{C\delta }$, $C_{l\delta }$, $C_{n\delta }$, $C_{l1\delta }$, $C_{l2\delta }$, etc.
General aerodynamic parameters	Basic data $C_{D0}$, $C_{L0}$, $C_{m0}$
	Static derivative $C_{D\alpha }$, $C_{L\alpha }$, $C_{m\alpha }$, etc.
	Dynamic derivatives $C_{lp12}$, $C_{np12}$, $C_{cp12}$, $C_{l1p12}$ $C_{Lq}$, $C_{L\dot{\alpha }}$, $C_{mq}$, $C_{m\dot{\alpha }}$, etc.	Nastran MSC program calculation

.

3. Multibody Dynamic Equation of the Composite Aircraft System

According to the established dynamic model of the wingtip-jointed composite aircraft system, the simulation platform was built based on MATLAB/Simulink with a multibody system state module, a multibody dynamics solver module, a force model calculation module, an automatic control module, and a record and report module, as shown in Figure 8.

The multibody system state module mainly takes charge of the corresponding parameters including mass, static moment, and rotation tensor of the augmented body under different real-time relative roll angles (${\phi }_{12}$). In the force model calculation module, the aerodynamic force and moment are calculated based on the multibody system state and the real-time motion parameters. The multibody dynamics solver module mainly takes charge of accepting the input parameters from the multibody system state module and force model calculation module. The automatic control module mainly conducts the design work of control strategy based on the real-time errors of the aircraft system’s position, posture, and velocity. The record and report module mainly saves the data and display in the scope at every iteration. The iterative format was selected as ode45 and the step size was 0.01 s.

The dynamic simulation platform of the composite aircraft system was built based on the equations with total parameter use, such that the dynamic response of the single- and composite-trimming strategies could be carried out. The single-trimming strategy refers to each control surface deflection of the composite aircraft needing to meet the composite requirements. The overall motion dynamic equation of the aircraft system is balanced, and there is no need to satisfy the control input scheme of the internal relative motion dynamic equation. The composite-trimming strategy means that the control input of each control surface deflection of the composite aircraft system needs to meet the balance of the overall motion and internal relative motion dynamic equations of the composite aircraft system. Detailed parameters and corresponding numbers are listed in

Table 3. Aerodynamic database acquisition method.

Motion Parameters	Number	Rudder Surface Parameters	Single-Trimming Strategy		Composite-Trimming Strategy
			Right Aircraft	Left Aircraft	Right Aircraft	Left Aircraft
Height	5000 m	Throttle angular displacement	58.03%	58.03%	58.70%	58.70%
Speed	95.91 m/s	Elevator deflection angle	3.68°	3.68°	4.07°	4.07°
Attack angle	2.27°	Aileron deflection angle	0°	0°	−3.80°	3.80°
Relative roll angle	0°	Rudder surface angle	0°	0°	0°	0°

.

According to the simulation results of the compound flight dynamics response with the single-trimming strategy, the system dynamics simulation of the composite aircraft is terminated at 4.59 s when the relative roll angle reaches 180 degrees; that is, the two composite aircraft finally collide in the form of relative roll motion at 4.59 s, through the wing tip hinge mechanism.

In contrast to the single aircraft flight dynamics response, the relative roll moment exists in composite flight, which results in an internal relative roll motion. It is necessary to deflect the aileron, elevator, or rudder, in order to balance the relative roll moment in the system internal dynamic equation. Differential dual elevators cause redundant pitching torque, which makes the balancing problem of the composite aircraft more complex. The angle of differential rudder causes an internal relative yaw moment, which leads to higher requirements for the strength of the wingtip hinge structure and the docking accuracy of the mechanism. Therefore, the composite-trimming strategy mainly uses the deflection angle of the differential aileron surface to balance the internal relative rolling moment $(l_2-l_1)$. This can be obtained from the aerodynamic database for composite motion characteristic parameters and control parameters of rudder components. The simulation results are shown in Figure 9. The divergence velocity of the Mach number, angle of attack, pitch angle, and relative roll angle were all less than in the single-trimming strategy simulation results. Meanwhile, the variation amplitudes of the zero-order parameters of the whole system and the system including the Mach number, angle of attack, pitch angle, and relative roll angle were all less than 0.3% within 10 s of simulation time, such that the aircraft system can be regarded as meeting the cruise requirements for straight and level flight.

4. Lift Drag Ratio Analysis of Composite System Cruise Process

During cruising, the fuel consumption rate of the engine is related to $T_R$ (the thrust required) by the engine and $W$ (the weight of the aircraft system) [23].

$$T_R=\frac{W}{K}.$$

(20)

The specific mathematical relationship shown in Equation (20) proves that a larger value of $K$ leads to a smaller value of $T_R$ and provides better cruising performance. Therefore, the aircraft system will fly longer and further with the same volume of fuel.

In order to analyze the distribution law of the lift–drag ratio with the relative roll angle, the full-trim strategy that meets the requirements of compound flight cruise was adopted.

Polynomial approximation is a commonly used tool to approximate functions. In the interval of seeking the minimum point of a function, we can use the value of the objective function at several points to construct a polynomial as the approximate expression of the objective function, then use the minimum point of the polynomial as an approximation for the original objective function. Iterative calculation is repeated until a satisfactory result is obtained; this method is called interpolation. Not only are the function values on the nodes equal, but the corresponding derivative values are also required to be equal; even the higher-order derivatives are equal. The interpolation polynomial satisfying this requirement is called a Hermite interpolation polynomial [24].

The initial conditions of the calculation example were set to a height of 5000 m and a Mach number of 0.282. The angle of attack varied from −4° to 10° with sampling interval of 1°. The ${\phi }_{12}$ sampling data include ±40°, ±20°, ±10°, and 0°, such that 105 sample data were obtained, according to the CFD modeling and composite flight dynamics model; these are shown as the discrete points in Figure 10.

Figure 10 also shows the lift-to-drag ratio curves for different relative roll angles obtained using the PCHIP method, where ${\alpha }_{Kmax}$ represents the corresponding attack angle for maximum K value. The curves prove ${\alpha }_{Kmax}\left(40°\right)$ > ${\alpha }_{Kmax}\left(-40°\right)$ > ${\alpha }_{Kmax}\left(20°\right)$ > ${\alpha }_{Kmax}\left(0°\right)$ > ${\alpha }_{Kmax}\left(-20°\right)$ = ${\alpha }_{Kmax}\left(10°\right)$ > ${\alpha }_{Kmax}\left(-10°\right)$. The data curves also demonstrate that the optimal relative roll angle for the maximum value of lift-to-drag ratio among all the different sampling relative roll angles was 0°, and the maximum $K$ value was 29.31 with an angle of attack of 2.27°. Compared with other relative roll angles, the lift-to-drag ratio of ${\phi }_{12}$ = 0° was always the largest; in the range of −2° to 5°, $K\left(40°\right)$ < $K\left(-40°\right)$ < $K\left(20°\right)$ < $K\left(-20°\right)$ < $K\left(10°\right)$ < $K\left(-10°\right)<K\left(0°\right)$, while the distribution rules in other ranges are complicated to summarize.

5. Conclusions

According to the structure of the aircraft system and the constraint form of the articulated hinge mechanism, the degrees of freedom (DOFs) for the compound flight were analyzed. The Rottenburg–Wittenburg method was used to derive the full format of the 7-degree-of-freedom composite body dynamics equations and motion equations. The simulation results under given conditions demonstrated that the single-aircraft trim strategy, which is similar to the single-aircraft flight trim strategy, was not able to keep the composite aircraft system stable for a long time. With the overall and internal dynamic equations and kinematics equations balanced, the new composite trimming strategy is adopted to meet the requirement for steady composite flight. Based on the composite trimming strategy, 105 sampled data points under a given angle-of-attack range and relative roll angle range were obtained, while the piecewise cubic Hermite interpolation method was used to obtain the lift–drag ratio distribution curves corresponding to seven different composite motion characteristic parameters (mainly the relative roll angle). The maximum lift–drag ratio and the corresponding angle of attack values under the seven composite motion characteristic parameter schemes were obtained, and the influences of different composite motion characteristic parameters on the lift–drag ratio distribution curves were compared. Among all of the different sampling relative roll angles, the optimal relative roll angle for the maximum value of lift-to-drag ratio was 0°, while the maximum $K$ value was 29.31 when the angle of attack was 2.27°. The optimal relative roll angle for the maximum value of lift-to-drag ratio and the distribution rules in different range are obtained, which provides reference for cruising characteristics analysis (especially for lift–drag ratio) of composite aircrafts.