1. 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 .
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.
2. Theoretical Modeling
2.1. 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
, material stiffness matrix, geometric stiffness matrix
and static tangent stiffness matrix
are defined.
2.2. End-Point Algorithm
The discretized dynamic equilibrium equation at the end of
time step is given by:
where
,
and
are, respectively, the nodal internal, inertial and external force vectors at time
.
According to Crisfield [
37],
can be expressed as:
where
l is the length of element;
is the consistent mass matrix;
is the tensor of mass moment of inertia;
is the vector of angular velocity at node
i in the body-attached elemental frame
; and
is the
skew-symmetric matrix.
is the vector of nodal accelerations, whose translational components are defined in global reference system
and rotational ones are given in the body-attached frame
. Furthermore,
is the diagonal matrix composed of
and
, where
is the
identity matrix. The superscript (
) is used only to distinguish the
elemental frame
and the
diagonal matrix.
Equation (
2) implies that the initial elemental reference frame
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
in the present framework is modified as:
and
is modified as:
where
and
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
and
to take into account the contributions of internal and inertial forces to the generalized tangent stiffness matrix (details on the derivation of
, which leads to the static tangent stiffness matrix
, are given in Refs. [
38,
39]).
The variation of
can be expressed as:
where
is the variation of the global generalized displacement
. The full expressions of the matrices
,
and
are given as:
where
,
and
are, respectively, the time integration parameters of the Newmark method [
40] and time step, which are given in
Section 2.4.1.
where
is the nodal reference frame of node
i at time
n.
is the matrix that relates the non-additive pseudo-vector changes to the additive ones
. Matrices
,
and
depend on the definition of the co-rotational element frame. They are all given in Refs. [
38,
39].
2.3. Mid-Point Algorithm
The dynamic equilibrium equation of the mid-point algorithm is given by:
where
,
and
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
and
are, respectively, the vectors of nodal velocities at step
n and
, whose translational components are defined in the global reference system
and the rotational ones in the body-attached frame
.
The nodal external and internal forces acting during the time step are represented by their mid-point values:
where
and
are the local nodal internal forces at time
n and
, respectively, and the matrices
and
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
and
are given as:
where
.
The mid-point static tangential stiffness matrix
is given by:
which comes from the variation of Equation (
14).
is the local linear stiffness matrix, and the geometric stiffness matrix
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
and the dynamic term
:
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
is redefined as:
where positive parameter
controls the numerical damping.
2.4. 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.
2.4.1. 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
, 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
are then substituted into Equation (
21), which finally leads to a set of equations of the form:
where
is the equivalent dynamic tangent stiffness matrix, which includes contributions from both the internal and inertial terms.
is the equivalent incremental loads, and
is the predicted incremental displacement from time step
n to
.
2.4.2. Corrector Step with Mid-Point Algorithm
Having solved Equation (
22) for
, the vectors of displacement, velocity and acceleration at step
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
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:
where
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
to
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].
2.5. 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
, 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 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 is obtained by translating the basic system to the center of gravity (G) of the cross-section, and then rotating with an angle to coincide with the principal axes of inertia .
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
,
, and
is the coordinates of shear center (C) in the local reference system for the mass matrix
. 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
is the elemental coordinate system of the
ith element expressed in the global reference system.
2.6. 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
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
, which translates and rotates continuously with the varying generalized displacements of the element. The element frame
in aeroelasticity is then taken as the body axes
for flight dynamics. Note that, according to the most usual practice in flight dynamics, the
x axis of
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
of a wing element is defined as:
where
,
and
are the base vectors of
.
With ground axes
and elemental body axes
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
and
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
and
are, respectively, the components of the elemental angular velocity and angular acceleration, which are defined in the element frame
, 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
,
and
are the components of
and
V is the norm of
.
2.7. 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.
is the lift curve slope and
and
are the lift and moment slopes due to flap deflection, respectively.
,
and
are the lift, moment and drag coefficients for zero angle of attack (AOA), respectively. Furthermore,
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
,
and
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
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.
2.8. 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
and
are the resultant force and moment of the whole aircraft, which can be expressed as:
where
,
,
and
are, respectively, the aerodynamic force and moment, thrust and gravity acting on the
ith node of the aircraft, and
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
and
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.
4. 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.