Torbjørn Cunis
email: torbjoern.cunis@onera.fr

Jean-Philippe Condomines

Laurent Burlion
email: laurent.burlion@onera.fr

Full-envelope Equations of Motion of the Generic Transport Model based on Piece-wise Polynomial Aerodynamic Coefficients

I. Introduction

Defined as any deviation from the desired flight-path [START_REF]Airplane Flying Handbook[END_REF], in-flight loss of control (LOC-I) includes upset situations such as stall, high and inverted bank angle, as well as post-stall spirals and rotations. As such, LOC-I has remained the foremost cause of fatal accidents for the last decades and still imposes the highest risk to aviation safety [START_REF]Unstable Approaches: Risk Mitigation Policies, Procedures and Best Practices[END_REF]. Today's autopilots are not capable of recovery from such situations. Normally, pilots are able to pull back the aircraft by reducing the angle of attack, but if the stall occurs suddenly due to vertical gusts, pilots often don't have enough time to react. Especially for operators of unmanned aerial vehicles (UAV) who lack of awareness for the current flight situation, recovering their drone in case of stall is hard. Flying into clouds, where vertical gusts are more likely but visual inspection of the flight condition is impossible, prospects of recovery are worse. This can lead to catastrophic consequences, establishing the necessity of autopilots which are capable of upset recovery.

With their unstable and highly non-linear characterizations, LOC-I situations require extensive control effort and adequate approaches. Non-linear behaviour of aircrafts in the post-stall flight regime has been investigated analytically [START_REF] Caroll | Bifurcation Analysis of Nonlinear Aircraft Dynamics[END_REF][START_REF] Jahnke | Application of Dynamical Systems Theory to Nonlinear Aircraft Dynamics[END_REF][START_REF] Goman | Application of bifurcation methods to nonlinear flight dynamics problems[END_REF][START_REF] Kwatny | Nonlinear Analysis of Aircraft Loss of Control[END_REF][START_REF] Engelbrecht | A Multi-mode Upset Recovery Flight Control System for Large Transport Aircraft[END_REF] and researchers developed control laws for upset recovery [START_REF] Burcham | Development and Flight Test of an Emergency Flight Control System Using Only Engine Thrust on an MD-11 Transport Airplane[END_REF][START_REF] Burcham | Manual Throttles-Only Control Effectiveness for Emergency Flight Control of Transport Aircraft[END_REF][START_REF] Urnes Sr | Flight Control for Multi-engine UAV Aircraft using Propulsion Control[END_REF][START_REF] Chang | Aircraft Trim Recovery from Highly Nonlinear Upset Conditions[END_REF][START_REF] Xargay | L 1 adaptive controller for multi-input multi-output systems in the presence of nonlinear unmatched uncertainties[END_REF][START_REF] Stepanyan | Stall Recovery Guidance Algorithms Based on Constrained Control Approaches[END_REF][START_REF] Engelbrecht | Automatic Flight Envelope Recovery for Large Transport Aircraft[END_REF][START_REF] Crespo | Analysis of Control Strategies for Aircraft Flight Upset Recovery[END_REF][START_REF] Richards | Vehicle Upset Detection and Recovery for Onboard Guidance and Control[END_REF][START_REF] Gregory | Flight Test of L 1 Adaptive Control Law: Offset Landings and Large Flight Envelope Modeling Work[END_REF]. For the recovery approaches jean-philippe.condomines@enac.fr found in literature as well as proposed by the authors [START_REF] Cunis | Nonlinear Analysis and Control Proposal for In-flight Loss-of-control[END_REF] are model-based, there is a need for reliable flight dynamics data.

Due to its rich and freely available data, NASA's generic transport model (GTM) is well-recognized in aerospace engineering community and widely used in literature [START_REF] Kwatny | Nonlinear Analysis of Aircraft Loss of Control[END_REF][START_REF] Engelbrecht | A Multi-mode Upset Recovery Flight Control System for Large Transport Aircraft[END_REF][START_REF] Engelbrecht | Automatic Flight Envelope Recovery for Large Transport Aircraft[END_REF][START_REF] Crespo | Analysis of Control Strategies for Aircraft Flight Upset Recovery[END_REF][START_REF] Richards | Vehicle Upset Detection and Recovery for Onboard Guidance and Control[END_REF][START_REF] Chakraborty | Nonlinear region of attraction analysis for flight control verification and validation[END_REF][START_REF] Gill | Upset Dynamics of an Airliner Model: A Nonlinear Bifurcation Analysis[END_REF][START_REF] Allen | Safe Set Protection and Restoration for Unimpaired and Impaired Aircraft[END_REF][START_REF] Cunningham | A Piloted Evaluation of Damage Accommodating Flight Control Using a Remotely Piloted Vehicle[END_REF][START_REF] Stepanyan | Loss-of-Control Mitigation via Predictive Cuing[END_REF]. Representing a 5.5 % down-scaled, typical aerial transport vehicle, the GTM provides an unmanned aerial vehicle [START_REF] Jordan | AirSTAR: A UAV Platform for Flight Dynamics and Control System Testing[END_REF] with exhaustive, full-envelope aerodynamic data from wind-tunnel studies [START_REF] Foster | Dynamics Modeling and Simulation of Large Transport Airplanes in Upset Conditions[END_REF][START_REF] Murch | Recent NASA Research on Aerodynamic Modeling of Post-Stall and Spin Dynamics of Large Transport Airplanes[END_REF][START_REF] Frink | Computational Aerodynamic Modeling Tools for Aircraft Loss of Control[END_REF]. An open-source, 6-degree-of-freedom aerodynamic model is available today for use in MAT-LAB/Simulink [START_REF]Flight Dynamics Simulation of a Generic Transport Model[END_REF].

Polynomial fitting of the aerodynamic coefficients provide a constructive method to define and evaluate models based on analytical computation due to their continuous and differentiable nature. However, none of the polynomial models published recently [START_REF] Kwatny | Nonlinear Analysis of Aircraft Loss of Control[END_REF][START_REF] Chakraborty | Nonlinear region of attraction analysis for flight control verification and validation[END_REF] represent the aerodynamic coefficients well in the entire region of the envelope. To overcome this problem, fitting of the GTM aerodynamic coefficients by piece-wise defined polynomials accounts for both the pre-stall and the poststall behaviour of the coefficients. On the other hand, a piece-wise defined model can bring the disadvantage of discontinuity and hence needs to be treated with care.

In this paper, we propose a novel approach for fitting aerodynamic coefficients, namely piece-wise polynomial identification, and demonstrate a trim condition analysis for the full-envelope aerodynamic model of the GTM based on piece-wise fittings.

II. Preliminaries

In this paper, we will mainly refer to the axis systems of the international standard [START_REF]Flight dynamics -Concepts, quantities and symbols -Part 1: Aircraft motion relative to the air[END_REF]: the body axis system-x f -axis aligned with fuselage, z f -axis points vertically down, y f -axis completes the setup-; the airpath axis system-x a -axis aligned with aircraft velocity vector V, y a -axis lies in the x a -y f -plane, z a -axis points down completing the setup; and the normal earth-fixed axis system -z g -axis points towards the center of gravity, x g -axis and y g -axis are earth-fixed completing the setup.

The orientation of the body axes with respect to the normal earth-fixed system is given by the attitude angles Φ, Θ, Ψ and to the air-path system by angle of attack α and side-slip β; the orientation of the air-path axes to the normal earth-fixed system is given by azimuth χ, inclination γ, and bank µ. (Fig. 1.)

x g x a x f z g z f mg L F D V Θ γ α

III. Equations of Motion

The aerodynamic forces and moments in the body axis system are subject to the aerodynamic coefficients,

  X A Y A Z A   f = 1 2 SV 2   C X C Y C Z   ;   L A M A N A   f = 1 2 SV 2   bC l c a C m bC n   +   X A Y A Z A   f × x cg -x ref cg ; (1)
where S, b, and c a are wing area, aerodynamic span, and aerodynamic mean chord; x cg and x ref cg are the position and reference position center of gravity, respectively. Likewise, the thrust forces and moments are given to

  X F Y F Z F   f =   F 0 0   ;   L F M F N F   f =   0 l t F 0   ; (2)
assuming the engines to be aligned with the x f -axis, symmetric to the x f -y f -plane, and deviated from the origin only by the vertical offset l t . The gravity force is finally given in the normal earth-fixed axis system as

  X G Y G Z G   g =   0 0 mg   . (3)
We have the resulting effective lift force, drag force, and side force by rotation of Eqs. (1) to (3) into the air-path axis system,

L eff = X A f sin α -Z A f cos α + X F f sin α -mg cos γ cos µ; (4)
D eff = -X A f + Z A f sin α + X F f cos α cos β -Y A f sin β + mg sin γ; (5)
Q eff = -X A f + Z A f sin α + X F f cos α sin β + Y A f cos β + mg cos γ sin µ. (6)
Here, the effective drag force-negative x a -axis-leads to a change in airspeed; the normal component of the effective lift force, opposing the gravity vector,-negative z g -axis-changes the inclination; and the effective side force-positive y a -axis-contributes to the side-slip angle:

V = - 1 m D eff ; (7)
γ = 1 mV L eff cos µ; (8)
β = 1 mV Q eff . (9)
The horizontal component of L eff in case of a coordinated turn acts as centripetal force pulling the aircraft radially inwards, leading to a change of the azimuth χ:

χ = 1 mV L eff sin µ. ( 10
)
For a symmetric plane (I xy = I yz = 0), the resulting moments in body axis are derived from Eqs. ( 1) and ( 2) and the conservation of angular momentum to [START_REF] Brockhaus | Flugregelung[END_REF]

L f = L A f + L F f -q r (I z -I y ) + p q I zx ; (11)
M f = M A f + M F f -p r (I x -I z ) -p 2 -r 2 I zx ; (12) N f = N A f + N F f -p q (I y -I x ) -q r I zx ; (13)
with the inertias I x , I y , I z , I zx . The changes of angular body rates are then given as

ṗ = 1 I x I z -I 2 zx (I z L f + I zx N f ) ; (14)
q = 1 I y M f ; (15)
ṙ = 1 I x I z -I 2 zx (I zx L f + I x N f ) . (16)
Conveniently, the normalized body rates p, q, r are used rather, with

  p q r  = 1 2V   b p c a q b r   ⇐⇒   p q r   = 2V   b -1 p c -1 a q b -1 r  ( 17
)
and the time derivatives

  ṗ q ṙ  = 1 V   1 2   b ṗ c a q b ṙ   - V   p q r    . ( 18
)
Rotating the body rates into the normal earth-fixed axis, we conclude with the change of attitude:

Φ = p + q sin Φ tan Θ + r cos Φ tan Θ; (19) Θ = q cos Φ -r sin Φ; (20) Ψ = q sin Φ cos -1 Θ + r cos Φ cos -1 Θ. (21)
IV. Piece-wise Polynomial Identification For the full-envelope, the aerodynamic coefficients of the generic transport model are given by angle of attack, side-slip, surface deflections, and normalized body rates,

(α i , β i , η i , qi , C X,i , C Z,i , C m,i ) 1≤i≤k ; ( 22
)
(α i , β i , ζ i , ξ i , pi , ri , C Y,i , C l,i , C n,i ) 1≤i≤k . ( 23
)
Inspired by the AERODAS model [START_REF] Spera | Models of Lift and Drag Coefficients of Stalled and Unstalled Airfoils in Wind Turbines and Wind Tunnels[END_REF], it seems appropriate to fit the aerodynamic coefficients by piecewise defined polynomials to address the issue of accurate fitting over the full range of the angle of attack:

C (α, β, . . . ) = C pre (α, β, . . . ) if α ≤ α 0 C post (α, β, . . . ) else (24) with C = C X C Y C Z C l C m C n T .
In order to ensure continuity of the aerodynamic coefficients-and thus the resulting equations of motion-over the entire domain of α, we have the additional constraint

C pre (α 0 , •) ≡ C post (α 0 , •) . ( 25
)
The sub-functions are chosen to be sum of polynomials

C i X = C i Xα (α) + C i Xβ (α, β) + C i Xη (α, β, η) + C i Xq (α, q) ; (26)
C i Y = C i Yα (α) + C i Yβ (α, β) + C i Yζ (α, β, ζ) + C i Yp (α, p) + C i Yξ (α, β, ξ) + C i Yr (α, r) (27)
with

C X = C X C Z C m T , C Y = C Y C l C n T
for i ∈ {pre, post}. Optimal polynomials are subject to the cost functionals

i i=1 C pre * (α i , •) -C * ,i 2 + k i=i +1 C post * (α i , •) -C * ,i 2 (28)
with the data of the GTM given as

(α i , . . . , C α,i , C β,i ) i (α i , . . . , C Xη,i , C Xq,i ) i (α i , . . . , C Yζ,i , C Yξ,i , C Y p,i , C Yr,i ) i
and the pre-stall boundary α = α i selected from the GTM data. Continuity of the single terms,

C pre * (α 0 , Ξ) = C post * (α 0 , Ξ) (29)
for all C * and inputs Ξ, then implies Eq. 25.

We reduce this problem to a linear least-square optimization problem [START_REF] Lawson | Solving Least Squares Problems, ser. Classics in Applied Mathematics[END_REF], where the solution q contains the coefficients of the polynomial; the constraints of continuity are transformed into matrix equalities in q.

The coefficients with respect to angle of attack are fitted by piece-wise defined curves with the limit α 0 initially been found solving

C pre Xα (α 0 ) = C post Xα (α 0 ) , (30)
where C pre Xα , C post Xα has been obtained with out constraints. The coefficients with respect to the normalized rates are fitted by piece-wise defined surfaces. Finally, the coefficients with respect to side-slip angle and elevator deflections can be approximated by simple polynomial surfaces and hyper-surfaces, respectively.

A. Piece-wise curve fitting

Introducing the vector of monomials P n (α) up to degree n, i.e.

P n (α) = 1 α • • • α n T , (31)
and denoting the length of P n by r [n], we can write a polynomial as scalar product

C α (α) = P n (α) , q (32)
with the vector of coefficients

q T = b 1 • • • b r[n] .
The optimal polynomial sub-functions C pre α , C post α with coefficients q pre , q post minimizing the costs in Eq. ( 28) are subject to the linear least-square problem:

1 find q minizing K * q -κ 2 2 under the constraint A * q = 0, (33)
where • 2 2 is the L 2 -norm, q is the extended vector of coefficients, q = q pre q post T ,

K is the data monomials matrix

K =           P n (α 1 ) . . . 0 P n (α i ) P n (α i +1 ) 0 . . . P n (α k )           , (35)
and the constraint of continuity is written as

P n (α 0 ) T -P n (α 0 ) T =A * q pre q post = 0 (36)
with α 0 given a priori.

B. Piece-wise surface fitting

Fitting a piece-wise defined polynomial surface, where the polynomial sub-functions are given by C pre β = P n (α, β) , q pre ; C post β = P n (α, β) , q post ; with

P n (α, β) = 1 α β • • • α n α n-1 β • • • β n T , (37)
is similar to the piece-wise curve fitting discussed before. However, as the pre-stall boundary is extended to a curve rather than a point, the constraint of continuity alters to ∀β ∈ R. P n (α 0 , β) , q pre = P n (α 0 , β) , q post . (38)

Separation of the assigned parameter α ≡ α 0 yields P n (α 0 , β) , q i = Λ T 0 P n (β) , q i = P n (β) , Λ 0 q i (39) with i ∈ {pre, post},

P n (β) = 1 β • • • β n T , and
Λ 0 =      1 α 0 α n 0 1 • • • . . . α 0 1      . ( 40
)
As Eq. (39) again resembles polynomial curves with coefficients Λ 0 q i , we have that P n (β) , Λ 0 q pre = P n (β) , Λ 0 q post for all β ∈ R if and only if Λ 0 q pre = Λ 0 q post . Hence, the constraint of continuity for piece-wise defined surfaces is written as Λ 0 -Λ 0 q pre q post = 0 (41) which is used for the constrained least-square problem of Eq. (33).

V. Full-envelope Trim Analysis

The system of equations of motion of section III is in a trim condition if and only if the airspeed and air-path are constant, i.e. V = γ = 0; the side force vanishes, β = 0; the body rates remain unchanged, ṗ = q = ṙ = 0; and the attitude is constant, Φ = Θ = 0. In order to allow a coordinated turn, we relax azimuth and heading to be χ = Ψ; that is, the heading is changed by a bank turn only (cos µ = 0) and only if the azimuth changes in the same manner. We now have airspeed V , inclination γ, and bank µ as continuation parameters, leaving angle of attack α, side-slip β, the normalized rates p, q, r, the surfaces ζ, η, ξ and thrust F as free variables. The attitude angles Φ and Θ are fully determined by the aforementioned angles.

VI. Conclusion

Fig. 1 :

 1 Fig. 1: Axis systems with angles and vectors.

https://mathworks.com/help/optim/ug/lsqlin.html