A Comprehensive Framework for Coupled Nonlinear Aeroelasticity and Flight Dynamics of Highly Flexible Aircrafts

: A framework to model and analyze the coupled nonlinear aeroelasticity and ﬂight dynamics of highly ﬂexible aircrafts is presented. The methodology is based on the dynamics of 3D co-rotational beams. The coupling of axial, bending and torsional effects is added to the stiffness and mass matrices of Euler–Bernoulli beam to capture the most relevant characteristics of a real wing structure. The ﬁnite-state aerodynamic model is coupled with the structural model to simulate the unsteady aerodynamics. A scheme of mixed end-point and mid-point time-marching algorithms is proposed and applied into the implicit predictor–corrector integration, where the end-point algorithm is used in the predictor step for efﬁciency and mid-point algorithm in corrector step for accuracy. The ground, body and airﬂow axes for ﬂight dynamics are re-deﬁned by the global and elemental ones for structural dynamics, followed by the redeﬁnitions of local Euler angles and airﬂow angles of each element. The framework can be used for quick analyses of ﬂexible aircrafts in conceptual and preliminary design phases, including linear and nonlinear trim, aerodynamic load estimation, stability assessment, time-domain simulations and ﬂight performance evaluations. The results show the payload mass and its distributions will signiﬁcantly affect the trim state and longitudinal stability of highly ﬂexible aircrafts. predictor–corrector integration for the ﬁrst time, by which the computation efﬁciency, accuracy and stability are improved. Two representative numerical benchmarks, especially a highly ﬂexible ﬂying-wing aircraft, are presented and analyzed with the trimmed states calculated and dynamic responses after ﬂap perturbations computed numerically. The results show that with the increase of the concentrated payload, large deformations occur in the highly ﬂexible wings, resulting in signiﬁcant changes in trim state and the reduction in phugoid mode stability. The effects of spanwise mass distribution on trim and longitudinal stability are further studied, indicating that, for highly ﬂexible aircrafts, the longitudinal performances will be affected by the spanwise mass distribution; span-loaded aircrafts with distributed payloads have much smaller deformations and higher longitudinal stabilities than the center-loaded ones with concentrated payloads. time-domain response simulations and further control evaluations. Comparisons with published literature show the validity and reliability of the present framework. The investigations on lateral and directional characteristics and symmetric and unsymmetrical gust responses with dynamic stall effects of large-scale ﬂexible aircrafts based on the present framework will be the focus of the future work.


Introduction
High-altitude long-endurance (HALE) aircrafts, which are being considered for aerial reconnaissance, long-time surveillance, environmental sensing and communication relays, have been the topic of modern flying vehicles [1]. Due to the energy and weight restrictions, HALE aircrafts, especially solar-powered aircrafts, notably feature very low wing-loading and ultra-light wing structure with high aspect-ratio. The long and slender wings inherently have relatively low structural stiffness, which gives rise to large geometrically nonlinear deformations during normal flights. The aerodynamics, inertia distributions and even the thrust locations and directions could be significantly affected, resulting in the changes of aeroelastic characteristics, trim conditions and the overall flight mechanics of the whole aircraft. The mishap of Helios Prototype in 2003 is a typical result brought by the complex interactions among multiple disciplines and the lack of adequate analysis tools [2].
Thus, a simple yet realistic model and analysis framework that can accurately assess the flight performances of a largely deformed highly flexible aircraft, and can fully account for the coupling between flight dynamics and both structural nonlinearity and unsteady aerodynamics, which play a key role in the vehicle's static and dynamic characteristics, are of paramount importance in the analysis and design of these highly flexible aircrafts [3][4][5]. Such models are typically built using nonlinear beam representations of the structure [6][7][8], with either unsteady strip theory [9][10][11][12][13] or unsteady vortex lattice aerodynamics [14,15].
Patil [16] pointed out that many kinds of analyses with different models and complexities can be applied to investigate the aeroelastic characteristics of a flexible aircraft. Firstly, with the linear stiffness and mass matrices, an eigenvalue analysis can be conducted to know the structural natural mode shapes and frequencies of the unloaded aircraft, which are the basics of further dynamic studies. Secondly, the aircraft can be trimmed by a nonlinear static aeroelastic approach to get the deformed positions under given aerodynamic loads, thrusts and gravity. Then, the aeroelastic equations can be linearized and eigenvalue analysis be conducted about this equilibrium state, which will give a more reasonable estimation of the stability of the deformed flying vehicle. However, in many cases, especially when the wing-loading and structural stiffness of the aircraft are low enough, such a linear stability analysis based on the hypothesis of small disturbance cannot give a precise prediction to the nonlinear dynamic responses of a flying aircraft under large disturbances. Thus, the final option to analyze the responses of a flexible aircraft is the direct solution of the nonlinear aeroelastic equations via time-marching numerical integration. This kind of analysis will give responses of an aircraft as a function of time, which can be considered as the realistic indications of its dynamic behaviors.
To investigate the nonlinear aeroelastic responses, trim solutions and flight stability of a complete highly flexible aircraft, Patil and Hodges [17] developed an analysis package named NATASHA (Nonlinear Aeroelastic Trim and Stability of HALE Aircraft). In their framework, the 2D Peters' finite state unsteady aerodynamic theory is used to calculate the aerodynamic loads. The geometrically exact fully intrinsic equations are used to model the highly flexible structure of the aircraft. With all the variables defined in the local reference system, the computational costs of coordinate transformation and assembly of global tangential stiffness matrix are reduced, with which the solution efficiency is improved. Their studies on a typical flexible flying-wing aircraft show that there are significant differences in flight dynamic characteristics between the deformed and undeformed aircraft.
A nonlinear aeroelastic simulation toolbox named UM/NAST (The University of Michigan's Nonlinear Aeroelastic Simulation Toolbox) was developed by Cesnik [18] and perfected by Shearer [19,20] and Su [13], aimed at the investigations of dynamic aeroelastic responses coupled with flight dynamics of highly flexible aircrafts. Their analyses are based on the strain-based nonlinear beams and Peters' 2D aerodynamic model. Two simplified models that can take the effects of skin wrinkling and aerodynamic stall into consideration are added into their framework.
A simulation framework called SHARP (Simulation of High-Aspect-Ratio Planes) was developed by Palacios and co-workers [21,22] to study the dynamic characteristics of highly flexible aircraft with large deformations. The displacement-based nonlinear dynamic formulations coupled with an unsteady vortex lattice aerodynamic model are employed to capture the large deformations and nonlinear responses of the flexible aircraft. With the computation of nonlinear static trim parameters and dynamic responses to maneuver and gust loads, this framework can be used to conduct stability analyses, dynamic response simulations, flutter suppressions, and gust load alleviations.
In addition, there are several more powerful nonlinear aeroelastic computer codes on the subject of computational aeroelasticity with high fidelity models and fluid structure interaction with reduced order models (ROMs) such as ARMA-ROM, Volterra-ROM, and POD-ROM. Many important works on this subject were done by Farhat et al. [23][24][25][26] at Stanford University (Mech and Aerospace Engr Dept). A linearized synchronous and staggered fluid-structure time-integration interaction algorithms [27][28][29] were introduced and applied to predict the transient aeroelastic response of a reduced-order model of a complete F-16 aircraft configuration at a given Mach number.
A rather comprehensive review, in which the up-to-date models, methodologies, applications, important observations and current challenges in the research of nonlinear aeroelasticity and flight dynamics of highly flexible aircraft are presented and discussed in detail, was recently carried out by Afonso et al. [30]. Readers are referred to this review for the state-of-the-art and further details on the design, analysis and performance predictions of this kind of highly flexible aircraft.
The geometrically exact intrinsic beam proposed by Hodges and strain-based nonlinear beam proposed by Cesnik are frequently used in the modeling of highly flexible wings. However, the displacement-based models, which are more intuitive and can be easily transported to the existing finite element programs, have rarely been utilized. As for the coupling between flight dynamics and aeroelasticity, many studies adopt the mean axes or its ramifications [31,32], which assume that the elastic deformations are small, and essentially still being the 6-DOF formulation. Even in those approaches in which large deformations are considered, the Euler angles and airflow angles are not defined locally, resulting in the costs of coordinate transformations of displacements and velocities between global and local frames.
The main contribution of this paper, which perfects the foundational work of forerunners, is to propose a new and more comprehensive framework to model and analyze the coupled nonlinear aeroelasticity and flight dynamics of highly flexible aircrafts, aimed at the requirements of time-domain simulations and performance evaluations in the preliminary design phase of HALE aircraft.
The starting point of this framework is the dynamics of 3D co-rotational nonlinear beam model, which is low-order and capable of capturing the nonlinear elastic deformations and rigid-body motions in a computationally effective formulation, targeted for the low computational costs and quick analyses in the preliminary design. The original equations of nodal inertial force vector of the element modified to enhance its versatility. The coupling between axial, bending and torsional effects that often arise in wings, especially when using composite materials, is added to the traditional Euler-Bernoulli beams in order to capture some of the most relevant characteristics of real wing structure.
The mixed time-marching algorithm is proposed and applied into the implicit predictor-corrector integration scheme for the first time, in which the end-point algorithm is used in the predictor step for efficiency and mid-point algorithm with numerical damping in corrector step for energy conserving and computation accuracy and stability. By this means, the time step of nonlinear numerical integration can be extended to as long as 0.1 s (note that this time step is not a universal reference and it depends on the frequencies of the structure and the external forces, and is a function of wing size, mass, stiffness, etc.) and the Newton-Raphson correction will converge in no more than five iterations with a tolerance of displacement of 1 × 10 −6 .
To couple the flight dynamics with nonlinear aeroelasticity in a more intuitive and straightforward way, via the spatial locations and orientations of every element, the ground, body and airflow axes for flight dynamics are re-defined analogously to the global and elemental reference frames for structural dynamics at the local element level, followed by the redefinitions of local Euler angles (roll, pitch and yaw angles) and airflow angles (angles of attack and sideslip) of each element. With all the variables defined in local reference systems, the computational costs of coordinate transformations can be reduced, which improves the solution efficiency.
The Peters' finite-state subsonic inflow model is coupled with the structural model to simulate the unsteady aerodynamics of the aircraft in flight, and it will be perfected by the ONERA dynamic stall model [33,34] in the following work.
As for the nonlinear trim of a complete flexible aircraft, the Newton-Raphson method is implemented with the forming of Jacobian matrix, which relates the resultant force and moment of the whole aircraft with the trim parameters.
The presented framework can be used for the quick analyses of a complete flexible aircraft in its conceptual and preliminary design phases, including linear and nonlinear trim, aerodynamic load estimation, stability assessment, as well as linear and nonlinear time-domain simulations and flight performance evaluations. With the proposed framework, numerical studies based on two representative benchmarks were conducted and carefully compared with published literature. With excellent agreements reached, the accuracy, stability, efficiency and reliability of the presented framework and its corresponding computing codes are validated.

Co-rotation Based Geometrically Nonlinear Structural Model
In the co-rotational formulation, beams are assumed to be linearly elastic, undergoing large translations and rotations but small strains. The rigid-body motions are separated from the strain-producing deformations at the local element level. It is assumed that the internal element behavior is linear, whereas nonlinearity is introduced via the co-rotational technique.
The kinematics and statics of 3D co-rotational beams were presented in the previous work of the authors [35,36], where the local and global degrees of freedom, nodal and elemental reference systems, tangential transformation matrix T, material stiffness matrix, geometric stiffness matrix K tσ and static tangent stiffness matrix K t,stat are defined.

End-Point Algorithm
The discretized dynamic equilibrium equation at the end of n + 1 time step is given by: where q i,n+1 , q mas,n+1 and q e,n+1 are, respectively, the nodal internal, inertial and external force vectors at time n + 1. According to Crisfield [37], q mas can be expressed as: where l is the length of element; M is the consistent mass matrix; I ρ is the tensor of mass moment of inertia; w i is the vector of angular velocity at node i in the body-attached elemental frame U e ; and S(·) is the 3 × 3 skew-symmetric matrix.p is the vector of nodal accelerations, whose translational components are defined in global reference system U g and rotational ones are given in the body-attached frame U e . Furthermore, is the diagonal matrix composed of U e and I 3×3 , where I 3×3 is the 3 × 3 identity matrix. The superscript ( * ) is used only to distinguish the 3 × 3 elemental frame U e and the 12 × 12 diagonal matrix. Equation (2) implies that the initial elemental reference frame U 0 should be in the same direction as the global one. It is no longer applicable when the two frames point in different directions. Thus, to make it more generic and adaptive, the inertial force vector q mas in the present framework is modified as: and U * e is modified as: where and U 0 is the initial elemental reference frame expressed in the global reference system. To solve such a system of algebraic non-linear equations with the Newton-Raphson method, it is necessary to compute the variations of vectors q i,n+1 and q mas,n+1 to take into account the contributions of internal and inertial forces to the generalized tangent stiffness matrix (details on the derivation of δq i,n+1 , which leads to the static tangent stiffness matrix K t,stat , are given in Refs. [38,39]).
The variation of q mas,n+1 can be expressed as: where δp is the variation of the global generalized displacement p. The full expressions of the matrices K mas1 , K mas2 and K mas3 are given as: where β, γ and ∆t are, respectively, the time integration parameters of the Newmark method [40] and time step, which are given in Section 2.4.1.
where U i,n is the nodal reference frame of node i at time n. H −1 (·) is the matrix that relates the non-additive pseudo-vector changes to the additive ones ∆ψ i . Matrices A, L(r 2 ) and L(r 3 ) depend on the definition of the co-rotational element frame. They are all given in Refs. [38,39].

Mid-Point Algorithm
The dynamic equilibrium equation of the mid-point algorithm is given by: where q e,m , q i,m and q mas,m are vectors of mid-point external, internal and inertial nodal forces, respectively. The mid-point inertial force vector in present framework is defined as (note it is different from the one given in [37]): whereṗ n andṗ n+1 are, respectively, the vectors of nodal velocities at step n and n + 1, whose translational components are defined in the global reference system U g and the rotational ones in the body-attached frame U e . The nodal external and internal forces acting during the time step are represented by their mid-point values: where q il,n and q il,n+1 are the local nodal internal forces at time n and n + 1, respectively, and the matrices T n and T n+1 are the corresponding tangential transformation matrices. The inertial contribution to the effective stiffness matrix is obtained from the variation of Equation (12): The full expressions of the matrices K mas1,m and K mas2,m are given as: where M r = l 6 I ρ . The mid-point static tangential stiffness matrix K tstat,m is given by: which comes from the variation of Equation (14). K l is the local linear stiffness matrix, and the geometric stiffness matrix K tσ takes the identical form to that of the statics [36]. The effective tangential stiffness matrix of the mid-point algorithm is then expressed as the sum of the static term K tstat,m and the dynamic term K tmas,m : In some cases, the approximate conservation of energy of the mid-point algorithm is not sufficient to ensure the numerical stability, which undergoes an uncontrolled growth of energy and finally is unable to reach convergence. A numerical damping correction is introduced as applied to the mid-point algorithm, in which, in place of Equation (14), the mid-point internal force q i,m is redefined as: where positive parameter ξ controls the numerical damping.

Mixed Time-Marching Integration Algorithm
The mixed time-marching algorithm is applied in the implicit predictor-corrector time integration scheme for the first time, in which the end-point algorithm is used in the predictor step for efficiency and mid-point algorithm with numerical damping in corrector step for energy conserving and computation accuracy and stability.

Predictor Step with End-Point Algorithm
Assuming all of the required information is available at step n, we can adopt the standard Newmark interpolation to get the vectors of velocity and acceleration at time n + 1, and then the corresponding internal and inertial forces. Here, we use the HHT-α time integration scheme, thus the end-point dynamic equilibrium equation (Equation (1)) is now expressed as: The values at time n and n + 1 are then substituted into Equation (21), which finally leads to a set of equations of the form: is the equivalent dynamic tangent stiffness matrix, which includes contributions from both the internal and inertial terms. (24) is the equivalent incremental loads, and ∆p is the predicted incremental displacement from time step n to n + 1.

Corrector
Step with Mid-Point Algorithm Having solved Equation (22) for ∆p, the vectors of displacement, velocity and acceleration at step n + 1 can be obtained by the Newmark interpolation relations; the corresponding internal and inertial force vectors are to be obtained successively.
The values at the mid-point can be interpolated and then substituted into the mid-point dynamic equilibrium equation (Equation (11)), which will, in general, lead to a residual g m that is not zero. Therefore, a standard process of Newton-Raphson (or modified Newton-Raphson or quasi-Newton) iteration is to be applied, which gives: is the equivalent dynamic tangent stiffness matrix of mid-point algorithm. Note it takes a different form to that previously given for the predictor step. By solving Equation (25), the improvement δp to p n+1 can be obtained, so that: The more detailed descriptions of 3D formulations for a spatial beam and the numerical implementations of the rotational updates, as well as the overall solution strategy, can be found in the co-rotational nonlinear finite element literature [37,39,[41][42][43].

Element Stiffness and Mass Matrices
Different from the traditional Euler-Bernoulli beams, the coupling among axial, bending and torsional effects is added into the element stiffness and mass matrices in the present framework. These effects are included by considering the flap-torsion (FT), lag-torsion (LT) and axial-torsion (AT) cross-stiffness properties in the stiffness matrix, and dynamic coupling terms between torsion and bending in the mass matrix.
As illustrated in Figure 1, there is a basic reference coordinate system O xyz , whose origin is chosen arbitrarily on the cross-section. The x axis lies parallel to the beam axis, and the y and z axes correspond, respectively, to the lagwise movement direction (or flapwise bending axis) and the flapwise one (or lagwise bending axis), which in general will not be the principal bending directions.
The local reference system for the stiffness matrix T xyz is parallel to the basic one with origin at the tension center (T) of the cross-section, and the x axis will be the neutral axis of the beam. It can be noted that T is generally different from the shear center (C). The local reference system for mass matrix G xy z is obtained by translating the basic system to the center of gravity (G) of the cross-section, and then rotating with an angle α m to coincide with the principal axes of inertia {y , z }.
The element stiffness matrix in the basic coordinate system is given as: where The element stiffness matrix in local axes [44] can be expressed as: in which Descriptions of the symbols in Equation (29) are given in Table 1. The element mass matrix in the basic coordinate system is: where where c(α m ) = cos(α m ), s(α m ) sin(α m ), and (y c , z c ) is the coordinates of shear center (C) in the local reference system for the mass matrix G xy z . The element mass matrix in local axes can be expressed as: in which Descriptions of the symbols in Equation (31) are given in Table 2.
Given the element stiffness and mass matrices for all the elements, the assembling of the global stiffness and mass matrices can be done in the usual way: where and U ei is the elemental coordinate system of the ith element expressed in the global reference system.

Coupling of Aeroelasticity and Flight Dynamics
In the traditional flight dynamics, the aircraft is modeled as a rigid body with six degrees-of-freedom. However, due to the different velocities and orientation angles at different locations of the flexible aircraft, it is no longer proper to represent the aircraft as a rigid body. The Euler angles (roll, pitch and yaw angles) and airflow angles (angles of attack and sideslip) for a flexible aircraft should be redefined at the local element level.
In the present framework, given the spatial locations and orientations of every element, the ground, body and airflow axes for flight dynamics are re-defined analogously to the global and elemental reference frames for structural dynamics at the local element level.
As illustrated in Figure 2, the ground axes U g for flight dynamics is defined to be identical to the global reference frame for aeroelasticity. For each element, there is a body-attached element frame U e , which translates and rotates continuously with the varying generalized displacements of the element. The element frame U e in aeroelasticity is then taken as the body axes U b for flight dynamics. Note that, according to the most usual practice in flight dynamics, the x axis of U b is directed to "forward", y axis to "rightward", and z axis to "downward", while, in aeroelasticity, the x axis of the element frame is preferred to lie along the axial line of the beam, and y and z axes along the principal axes of inertia. Thus, the body axes U b of a wing element is defined as: where U ex , U ey and U ez are the base vectors of U e .

Undeformed shape
Deformed shape With ground axes U g and elemental body axes U b defined, the Euler angles (roll, pitch and yaw angles) of each element can be defined in the usual way. The translational velocity and acceleration of the element in body axes can be expressed as: where V g and a g are, respectively, the inertial velocity and acceleration vectors expressed in the ground or global axes. The angular velocity and acceleration in body axes can be expressed as: where w x , w y , w z and a x , a y , a z are, respectively, the components of the elemental angular velocity and angular acceleration, which are defined in the element frame U e , in light of the co-rotational approach.
With the translational velocity in body axes, the angles of attack and sideslip of an element can be calculated as: where V bx , V by and V bz are the components of V b and V is the norm of V b .

Peters' Finite State Inflow Theory for Aerodynamic Modeling
The unsteady aerodynamic loads used in the current framework are based on the Peters' 2D finite state inflow theory [45]. The model calculates aerodynamic loads on a thin airfoil section undergoing large motions in an incompressible flow. The aerodynamic loads per unit span calculated about the aerodynamic center (A.C.) are given as: where ρ is the air density, b is the semichord, d is the distance of the mid-chord in front of the reference axis (here, it is the shear center (C)), and δ is the trailing-edge flap deflection angle. c lα is the lift curve slope and c lδ and c mδ are the lift and moment slopes due to flap deflection, respectively. c l0 , c m0 and c d0 are the lift, moment and drag coefficients for zero angle of attack (AOA), respectively. Furthermore, λ 0 is the inflow parameter, which accounts for induced flow due to free vorticity.
where λ is the vector of inflow states [45], which is given by: where A, b and c are constants given as: where N is the number of finite states. To transfer the loads from A.C. to the wing reference axis (here, it is shear center (C)), one may use l ra = l ac ; m ra = m ac + ( Furthermore, the aerodynamic loads are rotated to the body coordinate system, and then the global coordinate system, which yields: where l is the local span, which is here the length of the corresponding element. Finally, the elemental force and moment are allocated equally to the two nodes of each element.

Nonlinear Trim of the Flexible Aircraft
The nonlinear trim of the flexible aircraft is performed for zero resultant force and moment of the whole aircraft. A cost function is defined as: where F R and M R are the resultant force and moment of the whole aircraft, which can be expressed as: where F ai , M ai , F ti and F gi are, respectively, the aerodynamic force and moment, thrust and gravity acting on the ith node of the aircraft, and d i is the distance vector from the moment reference point (e.g., center of gravity or the constraint point of the aircraft) to the ith node.
The cost function f is minimized over the solution space using the body angle of attack α, elevator deflection angle δ and thrust T. Newton-Raphson iteration is employed to find the optimal trim parameter S, i.e., where S = [α, δ, T] T and ∂ f ∂S is the Jacobian matrix that relates the resultant force and moment of the whole aircraft with the trim parameters, which is computed numerically trough finite differences as: It is worth pointing out that, theoretically, Equation (49) is at risk of falling into a local minimum, but this situation is never seen to happen in practical applications.

Numerical Validation
Validations of the structural static solution under given loads and static aeroelastic responses with different flight conditions are presented in the previous work of the authors [35,36] and are not repeated here. To validate the implementation of the present framework and the corresponding computing codes, the dynamic responses of two representative numerical benchmarks, namely the cantilevered wing proposed by Goland [46] and the highly flexible flying wing proposed by Patil et al. [17], were analyzed and compared carefully to the results taken from Murua [22], Patil [17] and Su [13].

Cantilevered Goland Wing with Inertial Coupling
The cantilevered Goland wing was analyzed as a method-verification study for the present time-marching aeroelastic approach and the coupling of bending and torsional inertia. The properties of the Goland wing are given in Table 3. With the effects of gravity neglected, the wing starts from rest with a very small AOA of 0.05 degree. The flutter speed was estimated in the time domain with different freestream velocities evaluated until the flutter onset was found.  Figures 3 and 4, which, respectively, correspond to the damped, neutral and divergent states.
At the air speed of 131.4 m/s, the amplitude of oscillation decays and finally approaches to a steady state. As the air speed goes up to 136.4 m/s, a periodic motion with unchanged amplitude occurs, indicating the flutter speed is 136.4 m/s. The flutter frequency was calculated to be 69.2 rad/s, with a time-domain peak-finding technique. Above this speed, the wing oscillation tends to grow.   Table 4 compares the flutter speed and frequency obtained in the present framework to the results found in the literature. It can be seen that the flutter frequencies calculated with different models are almost the same, but the flutter speeds calculated with strip theory are approximately 18% lower than those with unsteady vortex lattice method (UVLM). This is because the two-dimensional strip theory cannot simulate the loss of lift near the tip of a wing with finite span, while the three-dimensional vortex lattice method can. Different from high aspect-ratio wings, the Goland wing has an aspect-ratio of only 3.33, in which the 3D effect cannot be neglected. Using the strip theory to calculate the flutter speed of 3D finite-span wings may result in underestimation. To use 2D theory for a 3D wing, a spanwise distribution coefficient was applied to correct for the tip loss of lift: where τ is a user-defined parameter that controls the spanwise lift distribution and s is the distance between wingtip and the spanwise location of interest. Despite the difference from 3D UVLM, excellent agreement with other models based on strip theory is reached, from which the sufficient accuracy of the established framework and its corresponding computing codes can be observed.

Flying Wing with Coupled Flight Dynamics and Aeroelasticity
To validate the implementation of the coupling of flight dynamics with nonlinear aeroelasticity and the corresponding computing codes, the trim and dynamic responses of a representative highly flexible flying-wing aircraft with a span of 72.8 m, a constant chord of 2.44 m and a payload varying from 0 to 227 kg on the central pod (note that the payload is not considered the entire mass at the point where the central pod locates; apart from the payload, there is also the structural mass of the central pod itself, which is 27.23 kg), as shown in Figure 5 and Table 5, were analyzed and compared carefully to the results taken from Patil [17] and Su [13], respectively.  At the flight speed of 12.2 m/s at sea level, the highly flexible aircraft with the varying concentrated payload was trimmed for zero resultant force and moment. Figure 6 shows the variations of the trimmed AOA, flap deflection angle (assuming constant throughout the span) and thrust per motor with the varying payload. It can be seen that the trim AOA increases with the payload, as the heavier is the payload, the higher is the lift required, and thus the higher are the AOAs. When the concentrated payload increases, the aerodynamic loads get higher, which leads to large bending deformations in the highly flexible wings, and the aircraft turns into a "U" shape, as illustrated in Figure 7.  As the aircraft deforms to a "U" shape, the aerodynamic center moves backward, which brings about a pitch-down moment, and the flap deflection decreases to reduce the additional pitch-down moment of the flap. It can also be seen that the trimmed thrust does not change significantly with the payload. This phenomenon can be attributed to the simplicity of the drag model, in which the profile drag and frictional drag of the aircraft are constant and the induced drag that is proportional to the lift is not considered. Trim   Figures 8 and 9 show, respectively, the trimmed configurations, AOAs and flap deflection angles obtained with different structural models and different payload distributions. It can be seen that, for the concentrated payload located at the central pod, results from rigid, linear and nonlinear structural models are quite different. The rigid model keeps its initial shape, while the AOA increases and the flap deflection decreases linearly with the payload. The difference in trimmed AOAs between linear and rigid structural models is brought by the aeroelastic deformation, especially the torsional deformations. Under the aerodynamic loads, the wings bend upward and twist positively, i.e., nose-up. With the positive torsion, a lower AOA is sufficient to provide the required lift, thus the linear model shows a lower AOA compared to the rigid one. On the other hand, as a follower force, the local lift is always perpendicular to the airfoil, with the directions changing continuously with the deformed configurations of the wing, so that it no longer acts in the vertical direction. Large bending deflections and the associated loss of vertical lift lead to the requirement of a higher AOA. This is the main reason the AOAs from the nonlinear model are higher than those from the rigid one when the payload exceeds 160 kg.

Effects of Structural Model and Mass Distribution on
Furthermore, it can also be seen in Figure 7 that the bending deflection at wingtip of the nonlinear structural model is about 53% larger than that of the linear one, which indicates higher loss of lift, and results in a higher AOA, as illustrated in Figure 8.
The flap is mainly used to balance the pitch moment. Since there are no deformations, the flap deflection of rigid model decreases linearly with the payload as a result of the increasing pitch-down moment brought by the increasing AOA. The nonlinear structural model has larger deformations than the linear one, resulting in a larger back-shift of the aerodynamic center and larger pitch-down moment, and the flap deflections are much lower than those given by the linear structural model, as shown in Figure 9.
In Figure 7, we can see that, if the payload is concentrated on the central pod, the flexible aircraft will perform a highly deformed "U" shape configuration, leading to the dramatic changes in its trim and further in flight dynamic performances. If we now spread the payload equally to the three pods, there will be an almost evenly distributed mass that is balanced by the lift, and the deformations will be nearly negligible even with the heaviest payload. The variations of trim AOA and flap deflection angle of a highly flexible aircraft will behave almost as a rigid one, despite the slight differences that are mainly brought by the torsional deformations, as shown in Figures 8 and 9.
Given the trim parameters and the deformed configurations, the aerodynamic characteristics, static and dynamic stabilities of the aircraft at each trimmed condition can be assessed accordingly, which are the topics of static aeroelasticity, and related studies were done in the previous works of the authors [35,36]. In addition, the comparison in Figure 6 shows excellent agreement with results from Patil [17] and Su [13], indicating the static aeroelastic characteristics of the flying wing calculated here are sufficiently accurate.

Stability Analysis
Time-domain simulations of the nonlinear dynamic responses of the highly flexible aircraft with three different payload masses (122, 152 and 227 kg) were performed. The aircraft was initially at its trimmed state, and perturbation was introduced by adding a pre-defined flap deflection change: the flap deflection angle was ramped linearly up to 5 • between 1 and 2 s, then ramped linearly back to its trimmed angle between 2 and 3 s, and finally kept constant.
The simulation was conducted by the direct solution of the complete nonlinear dynamic equations by implicit predictor-corrector schemes with the mixed time-marching algorithms and a time step of up to 0.1 s. Figures 10-12 show the variations of altitude, airspeed and AOA of the mid-span (the center) of the aircraft, respectively.
It can be seen in Figures 10 and 11 that, with a payload of 122 kg at the central pod, the phugoid mode of the highly flexible aircraft is neutrally stable, as the altitude and airspeed of the vehicle after perturbation oscillate with equal amplitudes. When the payload increases to 152 or 227 kg, the phugoid mode becomes unstable with increasing amplitudes of oscillations: the heavier is the payload, the higher is the frequency and the larger is the amplitude of the oscillation. For the differences between the two results given by Patil [17] and Su [13], Su attributed them to the differences of damping of the models used, which the authors think was not explained precisely. In the opinion of the authors of the present paper, it is more likely to be the difference of payload masses that leads to the different results.  S u ( P a y l o a d = 1 5 2 k g ) P a t i l ( P a y l o a d = 2 2 7 k g ) P r e s e n t ( P a y l o a d = 1 2 2 k g ) P r e s e n t ( P a y l o a d = 1 5 2 k g ) P r e s e n t ( P a y l o a d = 2 2 7 k g ) Figure 11. Airspeed of flight with different payload.
In Figures 10 and 11, the exchange between potential energy and kinetic energy of the highly flexible aircraft can be seen clearly from the out-of-phase variations between the altitudes and airspeeds. Besides, the aircraft altitude loss indicates a loss of energy due to the unstable phugoid mode and constant thrust. T i m e ( s ) S u ( P a y l o a d = 1 5 2 k g ) P a t i l ( P a y l o a d = 2 2 7 k g ) P r e s e n t ( P a y l o a d = 1 5 2 k g ) P r e s e n t ( P a y l o a d = 2 2 7 k g ) Figure 12. AOA at mid-span with different payload.
As shown in Figure 12, within several cycles of oscillations, the mid-span (also the whole aircraft) starts to experience very high AOAs in every cycle, which is the inevitable result of the excited unstable phugoid mode.
It is to be noted that, since stall is not yet included in the present framework, the results at high AOAs cannot be considered to be the real motion of the aircraft. Once the dynamic stall is included, one would see the different responses at the highest altitudes, where the lift will be insufficient to balance the weight and the aircraft drops with the increase in velocity, resulting in instantaneous higher AOAs.
In addition, it can be clearly seen that excellent agreements are made with the results supplied by Patil [17] for 227 kg payload mass and Su [13] for 152 kg payload mass, which demonstrates the validity and accuracy of the presented framework and its corresponding computing codes.

Effects of Mass Distribution on Longitudinal Stability
It is known from the results and analyses presented in Section 3.2.1 that the spanwise distribution of mass will affect the trimmed state of the aircraft, which will further influence the flight dynamic characteristics. It can be inferred that the mass distribution in spanwise will alter the longitudinal stability, changing the unstable phugoid mode into a stable one.
In this experiment, part of the 152 kg payload remained in the central pod, and the remaining mass was distributed equally to the left and right ones. Figure 13 shows the variations of the altitude of the mid-span of the highly flexible aircraft with, respectively, 100%, 95%, 90% and 80% of the total payload mass on the central pod and the rest on the side ones.
It can be clearly seen that, in the case of the same total payload, with the reduction of the mass at central pod, the amplitude of oscillation decreases dramatically; the frequency also decreases to some extent. When the mass on the central pod reduces to 80% of the total payload, the amplitude of oscillation after perturbation even begins to decay.
With the varying distributions of payload mass, other parameters (the airspeeds and AOAs) have the same trend as the altitude of the aircraft, which are not explain in details here. Generally, for a rigid straight wing without sweep and dihedral, the effects of spanwise mass distributions on longitudinal stability and dynamic responses are negligible. However, it clearly tells that the longitudinal stability and dynamic responses of a flexible wing will be significantly affected by the spanwise mass distributions. With the same mass of payload, a span-loaded aircraft is more stable than a center-loaded one.

Concluding Remarks
A new and comprehensive analytical framework for the coupled flight dynamics and nonlinear aeroelasticity of highly flexible aircrafts is presented. The theoretical methodology is based on the statics and dynamics of a 3D co-rotational geometrically nonlinear beam model for elastic deformations and rigid-body motions, with the original equations for nodal inertial force vector of element modified to enhance its versatility. The coupling of axial, bending and torsional effects is considered by adding the coupling terms into the element stiffness and mass matrices of traditional Euler-Bernoulli beam in order to capture the most relevant characteristics of real wing structure. The Peters' finite-state unsteady aerodynamic model is coupled with the structural model to simulate the unsteady aerodynamics in flight. The flight dynamics are coupled with nonlinear aeroelasticity in a more intuitive and straightforward way by re-defining the ground, body and airflow axes analogously to the global and elemental reference frames for structural dynamics at the local element level, followed by the definitions of local Euler angles and airflow angles of each element. The scheme of mixed end-point and mid-point time-marching algorithm is proposed and applied into the implicit predictor-corrector integration for the first time, by which the computation efficiency, accuracy and stability are improved.
Two representative numerical benchmarks, especially a highly flexible flying-wing aircraft, are presented and analyzed with the trimmed states calculated and dynamic responses after flap perturbations computed numerically. The results show that with the increase of the concentrated payload, large deformations occur in the highly flexible wings, resulting in significant changes in trim state and the reduction in phugoid mode stability. The effects of spanwise mass distribution on trim and longitudinal stability are further studied, indicating that, for highly flexible aircrafts, the longitudinal performances will be affected by the spanwise mass distribution; span-loaded aircrafts with distributed payloads have much smaller deformations and higher longitudinal stabilities than the center-loaded ones with concentrated payloads.
The present framework accounts for realistic design elements including arbitrary structural shapes and properties, concentrated masses, multiple thrusts and control surfaces and unsteady aerodynamics. It is adaptive and easy to accommodate with other models, and can be used for complete aircraft analyses in the conceptual and preliminary design phases, including the quick linear and nonlinear trim, aerodynamic load estimation, stability assessment, as well as linear and nonlinear time-domain response simulations and further control evaluations. Comparisons with published literature show the validity and reliability of the present framework. The investigations on lateral and directional characteristics and symmetric and unsymmetrical gust responses with dynamic stall effects of large-scale flexible aircrafts based on the present framework will be the focus of the future work.

Conflicts of Interest:
The authors declare no conflict of interest.