text:
David Quero
email: david.queromartin@dlr.de

Christoph Kaiser

Pierre Vuillemin

Charles Poussot-Vassal

A STATE-SPACE MODEL FOR LOADS ANALYSIS BASED ON TANGENTIAL INTERPOLATION

Keywords: Aeroservoelasticity, reduced order model, Loewner framework, tangential interpolation, gust load alleviation, H ∞ controller synthesis

In this work an approach for the generation of a generalized state-space aeroservoelastic model based tangential interpolation, also known as Loewner rational interpolation, is presented. The resulting differential algebraic system (DAE) system is reduced to a set of ordinary differential equations (ODE) by residualization of the non-proper part of the transfer function matrix. The generalized state-space is of minimal order and allows for the application of the force summation method (FSM) for the aircraft loads recovery, which shows a superior convergence when compared to the mode displacement method (MDM) for an increasing number of generalized coordinates for the cut loads recovery. Compared to the classical rational function approximation (RFA) approach, the presented method provides a minimal order realization with exact interpolation of the unsteady aerodynamic forces in tangential directions, avoiding any selection of poles (lag states). After a demonstration of the tangential interpolation techniques on the transcendental Theodorsen and Sears functions, the new approach is applied to the generation of an aeroservoelastic model for loads evaluation of the NASA Common Research model under atmospheric disturbances, showing an excellent agreement with the reference model in the frequency domain. Applications include the aerodynamic transfer function matrices generated by either potential flow or linearized computational fluid dynamics (CFD) solvers. The resulting aeroservoelastic model of minimal order is used for the design of an H ∞ -optimal controller for gust loads alleviation (GLA).

INTRODUCTION

One of the major tasks for aircraft certification is the consideration of dynamic load cases, in particular those caused by gust and continuous turbulence encounters, considered by the regulatory agencies [START_REF]Certification specifications and acceptable means of compliance for large aeroplanes[END_REF][START_REF]Federal aviation regulations part 25: Airworthiness standards: Transport category airplanes[END_REF]. In order to reduce the internal loads caused by these atmospheric disturbances without compromising the structural weight different load alleviation schemes are applied. To that aim, an accurate model suited for the design of advance control laws covering the complete flight envelope is of most relevance. Thus the need arises for a precise but at the same time efficient model description in order to develop suitable control laws for an efficient load reduction when encountering atmospheric disturbances. Moreover, modern transport aircraft are becoming more flexible and as a result the eigenfrequency values corresponding to the IFASD-2019-066 first flexible modes decrease and tend to interact with the rigid-body behaviour of the structure, requiring thus a model description containing both the rigid and flexible effects for an effective load alleviation.

Modern approaches of control theory require a description of the model in the time domain of minimum size. However, the aerodynamic models describing the flow over the aircraft are provided within the frequency domain in tabular form as they cannot usually be obtained in an explicit form and if they do, they contain transcendental functions which cannot be transformed in the time domain by a direct application of an inverse Laplace transform. Thus there is a need of converting the frequency domain description into time domain models.

The method of Roger [START_REF] Roger | Airplane math modeling methods for active control design[END_REF] and Abel [START_REF] Abel | An analytical technique for predicting the characteristics of a flexible wing equipped with an active flutter-suppression system and comparison with wind-tunnel data[END_REF] represents the three-dimensional unsteady subsonic aerodynamic flow by means of a least-squares technique with a series of poles which represent the aerodynamic lags due to the presence of the wake. All elements of the transfer function matrix share the same poles, which in turn have to be chosen in advance. Note that the least squares fit is ill-conditioned when increasing the number of real poles and prone to numerical instabilities [START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF][START_REF] Gustavsen | A robust approach for system identification in the frequency domain[END_REF]. For the representation of a gust excitation additional issues related to the time delay behaviour are encountered. In order to take into account the poles as parameters to the least-squares fit they can be further considered in an optimization problem [START_REF] Dunn | An analytical technique for approximating unsteady aerodynamics in the time domain[END_REF][START_REF] Tiffany | Nonlinear programming extensions to rational function approximations of unsteady aerodynamics[END_REF], ensuring a minimum error of the least-squares for a fixed number of poles. Karpel presented the minimumstate method [START_REF] Karpel | Design for active flutter suppression and gust alleviation using statespace aeroelastic modeling[END_REF], where an iterative least-squares process is applied to the aeroelastic system, reducing the total number of augmented states. This iterative process leads to a significantly increase of the computational effort which may be of up to three orders of magnitude [START_REF] Eversman | Consistent rational-function approximation for unsteady aerodynamics[END_REF].

In an effort to reduce the size of the aeroservoelastic systems formulated in the time domain, modern techniques of the control theory have been commonly applied to aeroelastic systems where the loads recovery is done by means of the mode displacement method (MDM) [START_REF] Pototzky | New and existing techniques for dynamic loads analyses of flexible airplanes[END_REF][START_REF] Karpel | Structural dynamic loads in response to impulsive excitation[END_REF][START_REF] Haghighat | Model-predictive gust load alleviation controller for a highly flexible aircraft[END_REF][START_REF] Liu | LQG based model predictive control for gust load alleviation[END_REF], even though the force summation method (FSM) method has a superior convergence of the cut loads values with respect to the number of structural modes considered. In this work the FSM method is considered for the recovery of the cut loads acting over the airframe, both in the time and in the frequency domains.

In general all the rational function approximation (RFA) methods introduce an error in the description of the aerodynamic loads due to the least squares fit. In order to avoid this error, different techniques based on rational interpolation may be used. More recently, the eigensystem realization algorithm (ERA) algorithm [START_REF] Juang | An eigensystem realization algorithm for modal parameter identification and model reduction[END_REF] has been widely used for system identification and model reduction [START_REF] Silva | Development of reduced-order models for aeroelastic analysis and flutter prediction using the CFL3Dv6.0 code[END_REF][START_REF] Gaitonde | Reduced order state-space models from the pulse responses of a linearized CFD scheme[END_REF] in the time domain, which presents two main limitations. Firstly, its application is limited to systems with a proper transfer function matrix and thus only the MDM method can be considered for the loads recovery. Secondly, if the system under consideration has a large number of inputs and outputs then a big computational effort is required to compute the singular value decomposition (SVD) of the dense Hankel matrix [START_REF] Kramer | Tangential interpolation-based eigensystem realization algorithm for MIMO systems[END_REF].

For the consideration of linear time invariant (LTI) multi-input multi-output (MIMO) systems with non-proper transfer function matrix a more general setting, namely that of the linear descriptor systems, is needed [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF]. The approach presented in this work is based on the Loewner framework for descriptor systems, which represents an extension of the rational interpolation at infinite frequency to a set of frequencies available at particular values [START_REF] Antoulas | Approximation of large-scale dynamical systems[END_REF].

In this work an approach to generate an aeroservoelastic state-space model suited for loads computation and control design which overcomes the above limitations and listed in

Formulation

Method Limitations [START_REF] Roger | Airplane math modeling methods for active control design[END_REF] [START_REF] Roger | Airplane math modeling methods for active control design[END_REF] Rational functional approximation (RFA) with fixed real poles -Use of real poles, lack of complex poles to represent oscillatory fluid modes typical in the transonic regime. -Increasing the number of poles may result in an ill-conditioned problem.

-Not appropriate for input delays (gust) [START_REF] Karpel | Design for active flutter suppression and gust alleviation using statespace aeroelastic modeling[END_REF] [START_REF] Karpel | Design for active flutter suppression and gust alleviation using statespace aeroelastic modeling[END_REF] Rational functional approximation (RFA) and nonlinear optimization -Use of real poles, lack of complex poles to represent oscillatory fluid modes typical in the transonic regime.

-Very high computational time. -Not appropriate for input delays (gust) [START_REF] Silva | Development of reduced-order models for aeroelastic analysis and flutter prediction using the CFL3Dv6.0 code[END_REF] [START_REF] Silva | Development of reduced-order models for aeroelastic analysis and flutter prediction using the CFL3Dv6.0 code[END_REF] Eigensystem realization algorithm (ERA) is presented. The proposed approach shall overcome all the current limitations, avoiding any selection of poles. Also, no iterative procedure is required for the precise formulation of the aeroservoelastic model in the time domain.

Regarding the design of proper controllers for load reduction, several methods have been applied for the design of GLA techniques [START_REF] Zhao | Gust load alleviation on a large transport airplane[END_REF][START_REF] Fonte | Design and validation of active gust load alleviation systems for aircraft[END_REF]. In the present work a feedback controller strategy based on measurements over the aircraft structure has been considered. The aeroservoelastic model presented in this work allows for the systematic consideration of the flexible effects into the controller design, which have been shown to play a role due to the interaction with the rigid-body behaviour of the aircraft [START_REF] Yamashiro | Reduction of flight control system/structural mode interaction[END_REF][START_REF] Kron | Robust 2-dof H-infinity controller for highly flexible aircraft: Design methodology and numerical results[END_REF]. Aouf et al. [START_REF] Aouf | H 2 and H ∞ optimal gust load alleviation for a flexible aircraft[END_REF] compared the effectivity of different controllers designed with H 2 , weighted-H 2 and H ∞ techniques, showing the superiority of both the weighted-H 2 and H ∞ controllers when compared to a regular H 2 design. Typically the available sensor measurements for the controller action are either structural displacements, velocities or accelerations. The aim of this work concerns the generation of a proper aeroservoelastic model and thus it is assumed that the wing root bending moment can be directly measured, leading to the design of an optimal H ∞ controller which makes use of the outer ailerons for the load reduction strategy. Additional sensors and actuators can be readily added to the developed model [START_REF] Fonte | Gust load alleviation for a regional aircraft through a static output feedback[END_REF].

The paper is organized as follows. In Section 2 the Loewner realization in connection with the theory of descriptor linear systems is presented and an application to the bidimensional unsteady incompressible flow over an airfoil containing transcendental functions is shown. A generalization to three-dimensional aeroservoelastic systems is then presented in Section 3. In Section 4 applications to the NASA common research model (CRM) / flutter reduced order assessment (FERMAT) model with aerodynamic models considering the doublet lattice (DLM)

IFASD-2019-066 potential method and high-fidelity CFD methods are presented. Additionally, an H ∞ -optimal controller for a GLA strategy is developed. Finally conclusions and future work are pointed out in Section 5.

GENERALIZED REALIZATION PROBLEM

Tangential interpolation

A descriptor linear system Σ is described in the time domain by a set of differential and algebraic equations (DAE), also denoted as generalized state-space system:

Σ : E ẋ (t) = Ax (t) + Bu (t) (1) y (t) = Cx (t) + Du (t) ,
where x (t) ∈ R n is the state vector, u (t) ∈ R nu the input vector, E, A ∈ R n×n , B ∈ R n×nu , C ∈ R ny×n and D ∈ R ny×nu are constant with E possibly singular. In this work regular systems are considered, that is, det (sE -A) = 0 except for a finite number of eigenvalues (which may have an infinite value), denoted by the set σ (E, A). The resolvent set ρ (E, A) is given by ρ

(E, A) = C \ σ (E, A).
The transfer function matrix H (s) of the system Σ is:

H (s) = C (sE -A) -1 B + D,
and is said to be proper if lim s→∞ H (s) is bounded (< ∞ for each of the its components) and strictly proper if lim s→∞ H (s) = 0. If the DAE system is split into the fast H ∞ (s) (including only the poles with infinite value) and slow subsystems H p (s) (including only the finite poles) and restricted system equivalence relations are applied so that H (s) = H p (s) + H ∞ (s), the non-proper part of the transfer function matrix H ∞ (s) can be written as a Neumann series expansion [START_REF] Sokolov | Contributions to the minimal realization problem for descriptor systems[END_REF]:

H ∞ (s) = C (sN -I) -1 B + D = D -CB - υ-1 j=1 s j CN j B,
where N is a nilpotent matrix with nilpotence index υ and coincides with the DAE index. This provides a direct relation between the so called DAE index υ and the non-proper part of the transfer function matrix H ∞ (s), with the DAE index ν one order less than that of the highest polynomial term.

The set (E, A, B, C, D) is called a realization of H (s). The realization is not unique and the one of the smallest possible order n is called minimal. The minimal realization used here is based on the rational interpolation for tangential data [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF][START_REF] Lefteriu | A new approach to modeling multiport systems from frequency-domain data[END_REF]. Tangential interpolation is a form of rational interpolation where the data is interpolated along particular directions. In particular, the data consist of the right interpolation data:

λ i , r i , w i | λ i ∈ C, r i ∈ C ny×1 , w i ∈ C nu×1 , for i = 1, ..., n r : IFASD-2019-066 Λ = diag [λ 1 , ..., λ nr ] ∈ C nr×nr R = [r 1 , ..., r nr ] ∈ C nu×nr W = [w 1 , ..., w nr ] ∈ C ny×nr ,
and the left interpolation data:

µ j , l j , v j | µ j ∈ C, l j ∈ C 1×nu , v j ∈ C 1×ny , for j = 1, ..., n l : M = diag [µ 1 , ..., µ n l ] ∈ C n l ×n l L = l T 1 , ..., l T n l T ∈ C n l ×ny V = v T 1 , ..., v T n l T ∈ C n l ×nu .
The transfer function matrix is evaluated at the values λ i , µ j and r i , l j are referred to as right and left tangential directions, while w i , v j are the right and left tangential data. The rational interpolation problem for tangential data aims at finding a realization (E, A, B, C, D) such that the associated transfer function H satisfies the right and left constraints,

H (λ i ) r i = w i , i = 1, ..., n r (2) l j H (µ j ) = v j , j = 1, ..., n l
which is achieved by means of the Loewner and shifted Loewner matrices. Next the structure of these two matrices is described.

The set Z of points in the complex plane Z = {z 1 , ...z nr+n l } and the corresponding values of the transfer function matrix is partitioned into the left and right data:

Z = {λ 1 , ..., λ nr } ∪ {µ 1 , ..., µ n l } , (3)
where the total number of sample points is n r + n l . The Loewner matrix is built as:

L =     v 1 r 1 -l 1 w 1 µ 1 -λ 1 • • • v 1 rn r -l 1 wn r µ 1 -λn r . . . . . . . . . vn l r 1 -ln l w 1 µn l -λ 1 • • • vn l rn r -lpwn r µn l -λn r     ,
The Loewner matrix can be also expressed in terms of the tangential controllability and observability matrices [START_REF] Lefteriu | A new approach to modeling multiport systems from frequency-domain data[END_REF]. If the directions r i , l j are selected generically (random in practice) the rank of the Loewner matrix L is equal to the rank of the underlying matrix E. The shifted Loewner matrix is the Loewner matrix corresponding to sH (s) and is built as:

L σ =     µ 1 v 1 r 1 -λ 1 l 1 w 1 µ 1 -λ 1 • • • µ 1 v 1 rn r -λn r l 1 wn r µ 1 -λn r . . . . . . . . . µn l vn l r 1 -λ 1 ln l w 1 µn l -λ 1 • • • µn l vn l rn r -λn r ln l wn r µn l -λn r     ,
Assuming that n r = n l and that det (L σ -sL) = 0 a minimal realization for the linear descriptor system is given by [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF][START_REF] Quero | A generalized statespace aeroservoelastic model based on tangential interpolation[END_REF]:

IFASD-2019-066 E = -L, A = -L σ , B = V, C = W, D = 0, (4)
and the associated transfer function H (s) = W (L σ -sL) -1 V satisfies the right and left interpolation conditions of Eq. 2. An additional SVD factorization of the Loewner matrix L allows for the further size reduction of the system, generating a reduced order model (ROM) of the descriptor system:

L = Y 1 Y 2 S 1 0 0 S 2 X * 1 X * 2 , (5)
where * denotes the conjugate transpose,

S 1 ∈ R r×r , S 2 ∈ R (n-r)×(n-r) and X 1 , X 2 , Y 1 , Y 2
are of appropriate dimensions. The reduced system of size r defined by

(E r , A r , B r , C r , D r )
is obtained by a Petrov-Galerkin projection and represents the best approximation to the full Loewner matrix in the Frobenius or 2-norm:

E r = -Y * 1 LX 1 , A r = -X * 1 L σ X 1 , B r = Y * 1 V, C r = WX 1 , D r = 0. (6)
Once the set (E r , A r , B r , C r , D r ) has been determined, this system may contain unstable poles.

Assuming that the aeroelastic system under consideration is stable, the system is then approximated by a stable one. This is done by means of an approximation in the Hardy space H ∞ [START_REF] Quero | A generalized statespace aeroservoelastic model based on tangential interpolation[END_REF][START_REF] Köhler | On the closest stable descriptor system in the respective spaces RH 2 and RH ∞[END_REF].

Application to transcendental aerodynamic transfer functions

The main task for the application of the FSM method for the cut loads recovery in an aeroservoelastic framework remains to properly describe the distributed unsteady aerodynamic forces in the time domain. In order to show the convenience of the tangential interpolation method described in Section 2.1 an application for the incompressible unsteady flow over a bidimensional profile with heave and pitch degrees of freedom is presented. Based on this application, a connection between the theory of DAE or descriptor systems and the unsteady aerodynamic flow is shown. In Section 2.2.2 including atmospheric disturbances is considered.

The classical unsteady model of Theodorsen [START_REF] Theodorsen | General theory of aerodynamic instability and the mechanism of flutter[END_REF] provides the lift and pitch moment coefficient (excluding the added-mass terms) to arbitrary input motions of the profile under the assumptions of inviscid flow, incompressibility and planar wake. As described by Brunton [START_REF] Brunton | Unsteady aerodynamic models for agile flight at low reynolds numbers[END_REF], several authors have constructed state-space models corresponding to the Theodorsen formulation. Note that excluding the models obtained by Peters [START_REF] Peters | Two-dimensional incompressible unsteady airfoil theory -an overview[END_REF], which requires up to eight states for an appropriate system representation, these state-space formulations do not include the added-mass terms. In this section a time domain model including also the added-mass terms is presented.

The total local lift c l (positive upwards) and local pitch moment at the quarter chord c m (positive nose up) can be described as function of the heave u z (positive down) and pitch u θ (positive nose up) motion in the frequency domain by means of an aerodynamic transfer function matrix:

c l (ω) c m (ω) = H t (ω) u z (ω) u θ (ω) ,
where the system output is y = c l c m T , the system input u = u z u θ T and the transfer function matrix is [START_REF] Hodges | Introduction to Structural Dynamics and Aeroelasticity[END_REF]:

IFASD-2019-066 H t (ω) =   -ω 2 πL ref U 2 ∞ + iω 2π U∞ C (k) ω 2 πaL 2 ref U 2 ∞ + iω πL ref U∞ (1 + (1 -2a) C (k)) + 2πC (k) ω 2 πL ref 4U 2 ∞ ω 2 1 8 -a 2 πL 2 ref 2U 2 ∞ -iω πL ref 2U∞   , (7)
with k the reduced frequency k = ωL ref /U ∞ , where L ref is a reference length (the half-chord) and U ∞ the true airspeed. The parameter a represents the pitch axis location with respect to the half-chord and ranges between -1 for the leading edge and 1 for the trailing edge. The Theodorsen function C (k) is given by

C (k) = H (2) 1 (k) H (2) 1 (k) + iH (2) 0 (k)
, with H 0 and H 1 the Hankel functions of first and second kind [START_REF] Theodorsen | General theory of aerodynamic instability and the mechanism of flutter[END_REF]. Interestingly, the pitch moment at the quarter chord position is independent of the Theodorsen function.

It is clear from Eq. 7 that the aerodynamic system is not proper, as the values lim s→∞ H (s) are not bounded. According to the description of Section 2.1 and to Eq. 7, the corresponding time domain description must be in the form of a descriptor system of index υ = 3. In Quero et al. [START_REF] Quero | A generalized statespace aeroservoelastic model based on tangential interpolation[END_REF] a strategy for the application of the Loewner realization to this aerodynamic system together with a reduction of the DAE index in order to obtain a regular ODE system including the time derivatives of the input has been presented.

In the case of bidimensional incompressible flow the DAE of index 3 can been reduced to a regular ODE provided the time derivatives of first and second order of the input signal are included. For compressible flow and due to the behaviour of the aerodynamic solvers at high frequencies this cannot be done in a straightforward manner and thus the dependence of the unsteady aerodynamic forces on the second time derivative resulting from the classical RFA fit [START_REF] Roger | Airplane math modeling methods for active control design[END_REF] for compressible flow may be seen as a low frequency residualization of its high frequency behaviour [START_REF] Fonte | Gust load alleviation for a regional aircraft through a static output feedback[END_REF].

Next the proposed Loewner realization is applied to two well-known transcendental aerodynamic transfer functions corresponding to the unsteady motion-induced and gust-induced aerodynamic incompressible flow, namely the Theodorsen function in Section 2.2.1 and the Sears function in Section 2.2.2.

Theodorsen function

In the case of bidimensional unsteady incompressible flow a general realization can be obtained which is then valid for all values of the parameters involved in Eq. 7, namely U ∞ , L ref and a. For a realization of the Theodorsen function the following strictly proper aerodynamic transfer function is defined:

H c (k) = C (k) - 1 2 ,
where C (k) represents the Theodorsen function which depends on the reduced frequency k.

Next the Loewner and shifted Loewner matrices are built and a realization is carried out: where the system (E cr , A cr , B cr , C cr , D cr ) is obtained after applying Eqs. 4 and 6 for a reduction on the number of states and a projection on the Hardy space H ∞ for the suppression of possible unstable poles on the aerodynamic transfer function H c . The resulting generalized state-space system is:

A c = E -1 cr A cr , B c = E -1 cr B cr , C c = E -1 cr C cr , D c = E -1 cr D cr = 0,
ẋc = A c x c + B c 2π U ∞ 0 U ∞ 1 L ref 1 2 -a 0 0   u u ü   c l c m = C c 0 x c +     0 π π U∞ πL ref U∞ 3 2 -a πL ref U 2 ∞ - πaL 2 ref U 2 ∞ 0 0 0 - πL ref 2U∞ - πL ref 4U 2 ∞ -1 8 -a 2 πL 2 ref 2U 2 ∞       u u ü   .
This generalized state-space form with a maximum time derivative term of order 2 for the input vector u corresponds to a DAE of index 3 [START_REF] Gerdin | Identification and estimation for models described by differentialalgebraic equations[END_REF]. Fig. 1 shows different realizations of the Theodorsen function C (k) in the complex plane for different values of the reduced frequency k, where a comparison with the Jones approximation [START_REF] Jones | The unsteady lift of a wing of finite aspect ratio[END_REF] involving two states is shown. For the Loewner realization the cases corresponding to a number of states equal to that of Jones (n = 2) and eight are depicted. It becomes clear that the Loewner realization for an increasing number of poles is in agreement with the Theodorsen function computed directly in the frequency domain (FD). The number of poles in the Loewner realization is increased by simply modifying the truncation after the SVD decomposition of Eq. 5. Note that unlike for the classical interpolation techniques [START_REF] Eversman | Modified exponential series approximation for the Theodorsen function[END_REF], the Loewner realization does not require the selection of any poles and no iterative optimization is involved. The matrices A c , B c , C c and D c obtained for the case where up to eight states are retained after the SVD decomposition of the Loewner matrix L are provided in Quero et al. [START_REF] Quero | A generalized statespace aeroservoelastic model based on tangential interpolation[END_REF].

Sears function

When subjected to a vertical atmospheric gust w g , the incremental lift coefficient (at the quarter chord position with no additional pitch moment) c l of the airfoil can be obtained in the frequency domain by means of the Sears function S (k), which depends on the reduced frequency k [START_REF] Filotas | Theory of airfoil response in a gusty atmosphere. Part I: Aerodynamic transfer function[END_REF]:

IFASD-2019-066 S (k) = 2 πk H (2) 0 (k) -iH (2) 1 (k)
Thus the aerodynamic transfer function matrix for the lift and pitch moment coefficients (c l and c m ) at the airfoil quarter chord position is:

c l c m = H s (ω) w g ,
where:

H s (ω) = 2πS (k) e -ik 0 .
The term S (k) e -ik is commonly referred to as the modified Sears function [START_REF] Giesing | Sears function and lifting surface theory for harmonic gust fields[END_REF] and its inverse Fourier transform provides the well-known Küssner function [START_REF] Garrick | On some reciprocal relations in the theory of nonstationary flows[END_REF]. A Loewner realization with 3 and 12 states of the Sears function S (k) is shown in Fig. 2, where a comparison with the three states approximation by Jones [START_REF] Jones | The unsteady lift of a wing of finite aspect ratio[END_REF] is also plotted. There is a very good agreement of the Loewner realization of 12 states with the reference data in the frequency domain (FD). The resulting state-space system in the time domain is:

ẋs = A s x s + 2π U ∞ B s w g c l c m = C s 0 x c . (8)
As in Section 2.2.1, the system (E sr , A sr , B sr , C sr , D sr ) is obtained after applying Eqs. 4 and 6 for a reduction on the number of states and a projection on the Hardy space H ∞ , this time on the modified Sears function S (k) e -ik . Thus, the matrices:

A s = E -1 sr A sr , B s = E -1 sr B sr , C s = E -1 sr C sr , D s = E -1 sr D sr = 0
correspond to the realization of the modified Sears function including the complex exponential factor e -ik . In Eq. 8 there is no need for an extended input vector including the time derivative of the gust velocity w g , as the magnitude of the Sears function tends to zero for increasing reduced frequencies.

AEROSERVOELASTIC SYSTEM FOR LOADS ANALYSIS

In this section the Loewner realization (tangential interpolation method) presented in Section 2.1 is applied to the more general case of an aircraft flying in compressible flow. Once the aerodynamic system is described in generalized state-space form in Section 3.1, the structural model is also considered and a coupled aeroservoelastic system for dynamic loads prediction including the FSM method for the cut loads recovery is presented in Section 3.2.

Loewner realization of the aerodynamic system

For the general case of three-dimensional unsteady compressible flow the transfer function matrix describing the aerodynamic forces cannot be explicitly obtained as for the bidimensional incompressible case, given by Eq. 7. A common technique for the computation of the subsonic compressible aerodynamic forces for loads application is the doublet lattice method (DLM), which provides the matrix of aerodynamic influence coefficients (AIC) at a set of particular frequencies. Due to the assumption of linear and inviscid flow, the DLM method is not appropriate to describe the transonic flow. In this region other methods such as the correction of the AIC matrices [START_REF] Quero | An Aeroelastic Reduced Order Model for Dynamic Response Prediction to Gust Encounters[END_REF] or linearized CFD solvers in the frequency domain can be used [START_REF] Kaiser | Aeroelastic gust load prediction based on time-linearized RANS solutions[END_REF]. An application including high-fidelity CFD data has been previously shown in Poussot-Vassal et al. [START_REF] Poussot-Vassal | Data-driven approximation of a high fidelity gust-oriented flexible aircraft dynamical model[END_REF].

In this work both DLM and high-fidelity linearized CFD methods are considered, see Section 4.

The input of the aerodynamic system is given by the vertical atmospheric disturbance w g at a particular spatial location (the aircraft nose) and by a set of generalized coordinates corresponding to the control surfaces u c and the structural modes u h of the aircraft structure (in vacuum) including the rigid-body and flexible modes, as the aeroservoelastic system is reduced by projection on the eigenvalues of the structure (modal truncation). Within the usual assumptions for aeroservoelastic modeling, the translational degrees of freedom u t of the aircraft structure do not cause any change in the aerodynamic forces and thus their time derivatives ut are considered as input to the aerodynamic system. The vector u rf contains the generalized coordinates corresponding to rotational rigid-body motions and to the flexible modes.

As output of the aerodynamic system the nodal forces and moments splined over the structural grid P ag by means of the transpose of the spline matrix H kg for the equivalence of work [START_REF] Beckert | Multivariate interpolation for fluid-structureinteraction problems using radial basis functions[END_REF] is chosen:

P ag = H T kg P ak .
The vector P ak represents the aerodynamic contribution to the nodal loads acting on the aerodynamic grid corresponding to the body surface and the spline matrix H kg transfers the structural

IFASD-2019-066
displacements to the aerodynamic grid over the body surface. With the aerodynamic distribution provided by P ag the cut loads distribution L F SM can be recovered over the load reference axis (LRA) by means of the FSM. The aerodynamic transfer function matrix of size N a × (N h + N c + 1) relating the input and output vectors is given in the frequency domain as:

P ag (ω) = H g (ω) H t (ω) H rf (ω) H c (ω)     w g (ω) iωu t (ω) u rf (ω) u c (ω)     , (9)
where H t (ω)) = T H H (ω) /iω with T H a matrix selecting the columns from H (ω) corresponding to the translational degrees of freedom u t . In this way the singularity caused by the translational motion at frequency zero can be avoided. The number of aerodynamic degrees of freedom is given by N a and N h and N c represent the number of structural modes and control surfaces respectively.

First a realization on the aerodynamic transfer function matrix as provided in Eq. 4 is done.

Once the realization ( Ê, Â, B, Ĉ, D) has been obtained the eigenvalues of the system are obtained by solving det(s Ê -Â) = 0. Note that at this stage there is no need for a SVD decomposition of the Loewner matrix as this realization is not the final one representing the aerodynamic system. Similarly, no stabilization by projection into the Hardy space H ∞ is enforced yet. The possible infinite poles due to the non-properness are excluded by setting a tolerance value S ∞ and the remaining finite poles are collected in the set S f :

S f = s ∩ {|s| ≤ S ∞ } | det(s Ê -Â) = 0 . (10)
The column of the aerodynamic transfer function matrix in Eq. 9 corresponding to the gust input H g (ω) is proper and no split into proper and non-proper part is required here. Similarly, no residualization is required when considering the gust disturbance as unique system input [START_REF] Poussot-Vassal | Data-driven approximation of a high fidelity gust-oriented flexible aircraft dynamical model[END_REF], as will be shown in Section 4. In order to split the proper and non-proper parts for the rest of the columns of the aerodynamic transfer function matrix H (ω) in a numerical stable way, a least-squares fit of a rational proper transfer function together with a polynomial part is done for each component of the transfer function matrix data:

H ij (ω) = Np k=1 R ij iω + p k + P 0,ij + P 1,ij iω -P 2,ij ω 2 , (11)
where p k , R ij can either be real or come in conjugate pairs, P 0,ij , P 1,ij , P 2,ij ∈ R and the sum runs over the number N p of finite poles p k given by Eq. 10 and which belong to the set S f defined in Eq. 10, p k ∈ S f . Note that compared to the classical least-squares fit of the RFA technique which uses real coefficients [3, 9, 10], Eq. 11 is able represent peaks of the aerodynamic transfer functions in the frequency domain which appear in the transonic flow regime [START_REF] Friedewald | Quasi-steady doublet-lattice correction for aerodynamic gust response prediction in attached and separated transonic flow[END_REF] by the presence of conjugate pair poles [START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF]. Also, the use of a number of poles N c as done here would lead to a prohibitive number of additional states if the classical approach of a least-squares fit were directly transformed into the time domain [START_REF] Smith | State-space model generation for flexible aircraft[END_REF].

In the next step the actual data after subtraction of the polynomial non-proper part instead of the strictly proper part of the fit given by the first term of Eq. 11 is used. Thus, there is no error introduced in the proper part of the aerodynamic transfer function matrix and the representation of the non-proper part by the polynomial part can be seen as a residualization of the high frequency behaviour of the aerodynamic transfer function matrix H (ω).

IFASD-2019-066 Next the proper transfer function matrix H a (ω) is obtained, H a (ω) = H (ω) -P 0 -P 1 iω + P 2 ω 2 . Similarly as done for the incompressible flow case, now the tangential interpolation method is applied to the proper transfer function matrix H a (ω) and the realization

(E a = I, A a = E -1 r A r , B a = E -1 r B r , C a = C r , D a = D r
) is obtained, where the SVD decomposition and the projection on the Hardy space H ∞ have been applied. For the sake of clarity, the above described steps are represented in Fig. 3.

The matrices B a and D a together with the input vector have to be properly completed to accommodate the non-proper polynomial part, see Eq. 12. The matrices corresponding to the polynomial part have been split in column blocks corresponding to the translational motion (P 0t , P 1t , P 2t ), rotational and flexible generalized coordinates (P 0rf , P 1rf , P 2rf ) and control surface deflections (P 0c , P 1c , P 2c ). The term corresponding to the third time derivative of the translational motion ... u t is explicitly neglected by setting P 2t = 0 and the submatrices in [B au1 B ah2 B ah1 B au2 0 0 0 0 0 0] have been defined for convenience in a subsequent reordering of terms.

ẋa = A a x a + B au1 B ah2 B ah1 B au2 0 0 0 0 0 0 u a (12)
P ag = C a x a + 0 P 0t P 0rf P 0c P 1t P 1rf P 1c 0 -P 2rf -P 2c u a ,

where

u a = [w g uT t u T rf u T c üT t uT rf uT c ... u T t üT rf üT c ] T .
Dropping the terms corresponding to zero submatrices and reordering, the generalized state-space model for the unsteady aerodynamic loads is:

ẋa = A a x a + B au1 B au2 0 0     w g u c uc üc     + B ah1 B ah2 u rf ut (13)
P ag = C a x a + 0 P 0c P 1c -P 2c     w g u c uc üc     + P 0rf D ah2 D ah3   u rf uh üh   , IFASD-2019-066
where u h = [u T t u T rf ] T and:

D ah2 = P 0t P 1rf D ah3 = P 1t -P 2rf
Note that some other approaches neglect the low frequency residualization of the non-proper part of the transfer function matrix by neglecting the high frequency content instead. For instance, Fonte [START_REF] Fonte | Design and validation of active gust load alleviation systems for aircraft[END_REF] considers the loads only up to a defined frequency value, introducing a lowpass filter with high bandwidth. In principle the time delay introduced by the casual filter would have then to be taken into consideration.

Aeroservoelastic system

In this section the interaction between the aerodynamic and structural parts are considered in order to obtain a complete aeroservoelastic formulation. As in Section 3.1, the transfer of the aerodynamic forces to the structural grid is done by the spline matrix H T kg . The aeroservoelastic problem can be formulated in the time domain as:

M hh üh + B hh uh + K hh u h =φ T gh     C a x a + 0 P 0c P 1c -P 2c     w g u c uc üc     + P 0rf D ah2 D ah3   u rf uh üh       ,
where M hh , B hh and K hh represent the generalized mass, damping and stiffness matrices respectively. The matrix φ gh contains a number N h (due to modal truncation) of the aircraft structural modes u h (in vacuum) including the rigid-body and flexible modes. As in Section 3.1, the vector u h = [u T t u T rf ] T represents the total number of generalized coordinates (including rigid-body and structural flexible modes) and u c the control surfaces. The subset u t ⊂ u h represents the translational rigid-body degrees of freedom and the subset u rf ⊂ u h the rotational rigid-body degrees of freedom together with the flexible modes.

In order to include the translational motion in the system output, the state vector x is augmented, x = [u T t u T rf uT h x T a ] T , and the following equation of first order in the time derivatives is obtained:

d dt         u t u rf uh x a         =     0 0 T t 0 0 0 T rf 0 0 A h1 A h2 A h3 0 B ah1 B ah2 T t A a         u t u rf uh x a     (14) +     0 0 0 0 0 0 0 0 0 α -1 φ T gh P 0c α -1 φ T gh P 1c -α -1 φ T gh P 2c B au1 B au2 0 0         w g u c uc üc     ,
with:

IFASD-2019-066 α = M hh -φ T gh D ah3 A h1 = α -1 -k hf T f + φ T gh P 0rf A h2 = α -1 -B hh + φ T gh D ah2 A h3 = α -1 φ T gh C a .
The matrix T rf selects the non translational degrees of freedom u rf from the generalized coordinates u h (that is, excluding the subset u t corresponding to the translational motion) and the matrix T t selects the translational degrees of freedom from the vector of generalized coordinates u h .

Once the aeroelastic system has been written in a generalized state-space form the loads can be recovered by the equilibrium of forces (FSM) [START_REF] Reschke | Integrated flight loads modelling and analysis for flexible transport aircraft[END_REF]:

L F SM = T cg -M gg φ gh üh + P ag ,
where the matrix M gg represents the physical mass matrix and the summation matrix T cg sums up the nodal forces and moments acting to the cut loads distribution at the LRA locations. Note that unlike for the MDM method, the FSM method requires not only the generalized aerodynamic forces φ T gh P ag but also the aerodynamic distribution P ag . In Eq. 15 the FSM method has been included in the output of the aeroservoelastic system:

    u t u rf P ag L F SM     =     I 0 0 0 0 I 0 0 0 P 0rf + D ah3 A h1 D ah2 + D ah3 A h2 C a + D ah3 A h3 0 T cg P 0rf + βA h1 T cg D ah2 + βA h2 T cg C a + βA h3         u t u rf uh x a     (15) +     0 0 0 0 0 0 0 0 0 P 0c + D ah3 α -1 φ T gh P 0c P 1c + D ah3 α -1 φ T gh P 1c -P 2c + D ah3 α -1 φ T gh P 2c 0 ΓP 0c ΓP 1c -ΓP 2c         w g u c uc üc     ,
where:

β = -T cg M gg φ gh + T cg D ah3 Γ = βα -1 φ T gh + T cg .
Eqs. 14 and 15 represent the generalized state-space and output equations for the aeroservoelastic model description in the time domain and Fig. 4 depicts the corresponding aeroservoelastic system including the cut loads obtained by the FSM. There, the matrix T h selects the subset u h from the vector [w g u T c u T h ] T .

APPLICATION TO THE NASA COMMON RESEARCH MODEL CONFIGURA-TION

In this section the structural FERMAT model as created by Klimmek [49] for the NASA CRM configuration [START_REF] Vassberg | Development of a common research model for applied CFD validation studies[END_REF] has been used. The mass case C2 corresponding to a 100 % of the fuel mass with a mass equal to the maximum takeoff weight MTOW = 260000 (kg) and a center of gravity position in x direction (see Fig. 5 for axes definition) with respect to the aircraft nose IFASD-2019-066

Figure 4: Aeroservoelastic system with the force summation method (FSM) for the cut loads recovery. position x 0 of x cg -x 0 = 33.714 (m) is considered. For the DLM case both open-loop and closed-loop configurations after the synthesis of an H ∞ -optimal controller are presented. Note that the cut loads presented in this section are incremental and do not include the loads acting at the 1g or level flight condition.

The aeroservoelastic system as described in Section 3.2 has been generated for the flightpoint corresponding to a Mach number of M ∞ = 0.86 and an altitude h = 9100 (m). The true airspeed is U ∞ = 260.892 (m/s) and the reference length for the reduced frequency definition L ref = 7 (m). In Section 4.1 the aerodynamic loads as obtained by the DLM method are used, whereas in Section 4.2 the high-fidelity linearized Reynolds-averaged Navier-Stokes (RANS) equations are considered.

DLM

As stated in Section 3.2, the spline from the structural model deformation to the aerodynamic grid is carried out by by the transpose of the spline matrix corresponding to a nearest-neighbour force mapping. Neglecting the aerodynamic modeling of the fuselage component there are a total of 1910 aerodynamic panels representing the complete configuration. The wing semi-span is b/2 = 29.38 (m). A view of the right wing component is provided in Fig. 5. Note that a strip can include panels from the wing component as well as from the aileron if they are topologically connected. The ailerons are taken into account in Section 4.1.2, where the effect of a gust loads alleviation (GLA) strategy is evaluated using the proposed aeroservoelastic model. In that case, the control surface modes corresponding to the ailerons are modeled with the transpiration approach as described in Section 3.1.

Regarding the dynamic structural model, two symmetric rigid-body modes corresponding to the heave and pitch (around the center of gravity) motion together with the first symmetric wing bending mode have been taken into account. The heave forms the vector u t and the pitch together with the first symmetric mode form the vector u rf in Eq. 15. The interpolation of the precomputed AIC matrices to the set of reduced frequencies selected for the Loewner realization is done using a special linear interpolation [START_REF] Msc | Aeroelastic Analysis User's Guide Version 68[END_REF]. The tolerance in Eq. 10 has been set to S ∞ = 10 (U ∞ /L ref ). complex plane due to a gust disturbance input for the different values of the reduced frequency k. In Fig. 7 the pitch motion around the center of gravity is taken as input. For comparison purposes the results obtained with six real poles applying the classical Roger's approach (RFA) is also shown. The data labeled as Loewner shows the present method and FD refers to the reference data in the frequency domain. It is clear from Figs. 6 and 7 that the present approach shows a better match compared to the Roger's approach. Because of the finer resolution for the low reduced frequency range in the set of reduced frequencies selected, this better match is more evident in this region which is of significant importance for aeroelastic applications.

Open-loop

Two different aeroservoelastic models are generated. Both contain two symmetric rigid-body modes corresponding to the heave and pitch motions but differ in the number of symmetric flexible modes considered. The first one includes the first flexible mode corresponding to the wing bending with a natural frequency of 1.0574 (Hz) resulting in a total of 96 states and the second one includes 16 symmetric flexible modes contained in a frequency interval which ranges up to a value of 16.3742 (Hz) with a total number of 231 states. These time-domain aeroservoelastic modes are validated in open-loop against the reference aeroelastic model available in the frequency domain. For the validation, the system is excited with a 1-cosine gust w g (t) in (m/s) given by (at the aircraft nose) [START_REF]Certification specifications and acceptable means of compliance for large aeroplanes[END_REF] of amplitude w 0 (m/s), translational speed U ∞ (m/s) and half length H (m). For this application the gust amplitude has been set to α eq = w 0 /U ∞ equivalent to both include two symmetric rigid-body modes. In both cases the aeroservoelastic model is able to reproduce the complete time history of the incremental cut loads distribution along the entire component.

Closed-loop

Once the aeroservoelastic model has been validated in open-loop, the implementation of a GLA strategy is shown. Here an H ∞ -optimal strategy for controller synthesis 53] including the two rigid-body modes and the first flexible mode has been considered. The name H ∞ makes reference to the Hardy space H ∞ where the optimization takes place. order of the controller is N K = 7 and all controller parameters are set to free for the optimization problem. Additionally, a loop-shaping by means of weighting performances has been applied [START_REF] Zhou | Robust and optimal control[END_REF]. The H ∞ controller synthesis is possible due to the generation of the complete aeroservoelastic model with a number of states in the order of hundreds, which is an acceptable size for the open-loop system [START_REF] Apkarian | Nonsmooth H ∞ synthesis[END_REF]. The objective is to minimize the incremental bending moment at the wing root position by means of a symmetric deflection of the outer ailerons.

As discussed in Section 3, the application of the FSM for the cut loads recovery results in a non-proper part of the corresponding transfer function components. In this work the weight performances applied to the H ∞ controller synthesis consider the transfer function up to a certain frequency value and thus neglect the high frequency behaviour for the optimization problem to be solved for the controller design. For an extension including the non-proper part of the transfer function caused by the time derivative of the aileron input see Poussot-Vassal et al. [START_REF] Poussot-Vassal | A structured approach for input derivative index-2 descriptor dynamical model control[END_REF].

In this case the first model of Section 4.1.1 including two symmetric rigid-body modes and one flexible mode has been considered. The selected flightpoint remains unchanged and a corresponding H ∞ -optimal controller is designed. The aeroelastic response to a gust input of amplitude α eq = 1 (deg) and gust gradient H = 350 (ft) has been considered next. For this specific gust input, the maximum outer aileron deflection produced by the controller is 1.93 (deg). In Fig. 12 the incremental wing root bending moment is shown as a function of time for both the open-loop (OL) and the closed-loop (H ∞ ) cases. By deploying the outer ailerons in a symmetric way, the controller is able to reduce the maximum value of the incremental wing root bending moment caused by the gust M X in 5.6 % for the maximum peak and in 10.4 % for the minimum peak at the expense of increasing the incremental torsional moment M Y in 3.8 % for its maximum peak value, see Fig. 13. Further, the proposed aeroservoelastic model enables the evaluation of the effects of the applied GLA strategy on the cut loads over the complete aircraft structure.

High-fidelity CFD

The aerodynamic half-model used for the high-fidelity computations is based on the publicly available grid from the 4th AIAA Drag Prediction Workshop [START_REF] Vassberg | Development of a common research model for applied CFD validation studies[END_REF] and includes the wing, horizontal tail plane and fuselage components. The hybrid unstructured and structured volume mesh contains around 3.7 million points, of which 1000014 are located over the component surface. The farfield boundary is a semisphere with a radius of 757 (m).

The aerodynamic shape corresponds to a cruise flight condition with a global lift coefficient of C L = 0.5 at a Mach number of M ∞ = 0.86, an altitude of h = 9100 (m) and a steady angle of attack α 0 = 1.641 (deg). The corresponding steady pressure coefficient distribution c p is shown in Fig. 14, where a recompression shock wave is present at the upper wing surface. However, the flow remains fully attached.

For the high-fidelity unsteady aerodynamic response the linear frequency domain (LFD) solver for the RANS equations of the DLR TAU-Code has been chosen [START_REF] Kaiser | Aeroelastic gust load prediction based on time-linearized RANS solutions[END_REF][START_REF] Thormann | Linear-frequency-domain predictions of dynamic-response data for viscous transonic flows[END_REF]. The aircraft is considered rigid and thus only the aerodynamic contribution without the aeroelastic response is considered, which is equivalent to setting all the terms to zero in the state-space vector except for x a in Eqs. 14 and 15.

The same gust input with H = 350 (ft) as in Section 4.1.2 is considered. When considering the high-fidelity modeling of the aerodynamic forces, a higher number of states are required in the Loewner realization when compared to the aerodynamic forces predicted by the DLM IFASD-2019-066

CONCLUSIONS

In this work an aeroservoelastic model for the prediction of the aeroelastic response together with the dynamic loads distribution over the complete aircraft structure has been presented. It represents an interpolation instead of an approximation of the reference model in the frequency domain and avoids any selection of poles as for the classical methods. Thus, the aerodynamic forces distribution caused by the gust disturbance and induced by the aircraft motion can be properly described in the time domain. Once the aerodynamic model has been obtained in the time domain the structural model is taken into account and an aeroservoelastic state-space model of first order in time which includes the effect of control surfaces deflection has been derived.

Aerodynamic forces obtained in the frequency domain from solvers with different fidelity levels have been considered for the NASA CRM / FERMAT configuration, showing that the number of states required for a proper representation of the aerodynamic model increases with the level of complexity. Additionally, the design of a H ∞ -optimal controller with the generated aeroservoelastic framework has been demonstrated, while at the same time the presented model is able to assess the cut loads distribution over the complete aircraft.

Related future work includes:

• Generation of a complete aeroservoelastic framework for high-fidelity aerodynamic methods, as presented in this work for the case of DLM solvers. • Application of the developed framework to flutter stability analysis.

• Inclusion of parametric generalized state-space aeroservoelastic models as an alternative to the classical gain scheduling approach [START_REF] Antoulas | On two-variable rational interpolation[END_REF]. • Extension to nonlinear generalized state-space formulation of a nonlinear aeroservoelastic system. In this case the Loewner framework in connection with a functional or Volterra series expansion theory can be followed [START_REF] Gosea | Model reduction of linear and nonlinear systems in the Loewner framework: A summary[END_REF].

Figure 1 :

 1 Figure 1: Theodorsen function in the complex plane and corresponding Jones and Loewner approximations.

Figure 2 :

 2 Figure 2: Sears function in the complex plane and corresponding Jones and Loewner approximations.

Figure 3 :

 3 Figure 3: Flow diagram for the Loewner realization of the aerodynamic transfer function matrix considering the proper and non-proper parts.

Fig. 6 Figure 5 :

 65 Fig.6shows the element of the aerodynamic transfer function matrix corresponding to the incremental lift coefficient over the strip element at the wing station y = 13.48 (m) in the

Figure 6 :

 6 Figure 6: Element of the transfer function matrix from the gust input to the lift coefficient for the strip at wing station y = 13.48 (m).

Figure 7 :

 7 Figure 7: Element of the transfer function matrix from the pitch input to the lift coefficient for the strip at wing station y = 13.48 (m).

Figure 8 :

 8 Figure 8: Incremental bending moment M X over the right wing component at different time instants for M ∞ = 0.86, h = 9100 (m), H = 50 (ft), α eq = 1 (deg). Results obtained in the frequency domain (FD) and with the generalized state-space model in the time domain (TD) by the FSM with one flexible mode. The number of states is 96.

Figure 9 :

 9 Figure 9: Incremental bending moment M X over the right wing component at different time instants for M ∞ = 0.86, h = 9100 (m), H = 50 (ft), α eq = 1 (deg). Results obtained in the frequency domain (FD) and with the generalized state-space model in the time domain (TD) by the FSM with 16 flexible modes. The number of states is 231.

Figure 10 :

 10 Figure 10: Incremental bending moment M X over the right wing component at different time instants for M ∞ = 0.86, h = 9100 (m), H = 350 (ft), α eq = 1 (deg). Results obtained in the frequency domain (FD) and with the generalized state-space model in the time domain (TD) by the FSM with one flexible mode. The number of states is 96.

Figure 11 :

 11 Figure 11: Incremental bending moment M X over the right wing component at different time instants for M ∞ = 0.86, h = 9100 (m), H = 350 (ft), α eq = 1 (deg). Results obtained in the frequency domain (FD) and with the generalized state-space model in the time domain (TD) by the FSM with 16 flexible modes. The number of states is 231.

Figure 12 :Figure 13 :

 1213 Figure 12: Time history of the incremental bending moment M X at the wing root with and without GLA strategy.

Figure 14 :

 14 Figure 14: Steady pressure coefficient distribution c p for M ∞ = 0.86, h = 9100 (m) and α 0 = 1.641 (deg). NASA CRM model.

Figure 15 :

 15 Figure 15: Element of the aerodynamic transfer function matrix from the gust input to the the vertical force for a node at position x = 30.1076 (m), y = 3.6521 (m), z = 3.9362 (m).

Figure 16 :

 16 Figure 16: Incremental bending moment M X over the right wing component at different time instants for M ∞ = 0.86, h = 9100 (m), α 0 = 1.641 (deg), H = 350 (ft), α eq = 1 (deg). Results obtained by LFD in the frequency domain (FD) and with the generalized state-space model in the time domain (TD) for a rigid aircraft. The number of states is 300.

Table

Table 1 :

 1 Summary of methods applied for the formulation of aerodynamic and aeroservoelastic models in the time domain and current limitations.

ACKNOWLEDGMENTS

This work has been funded within the frame of the Joint Technology Initiative JTI Clean Sky 2, AIRFRAME Integrated Technology Demonstrator platform AIRFRAME ITD (contract N. CSJU-CS2-GAM-AIR-2014-15-01 Annex 1, Issue B04, October 2nd, 2015) being part of the Horizon 2020 research and Innovation framework program of the European Commission.

COPYRIGHT STATEMENT

The authors confirm that they, and/or their company or organization, hold copyright on all of the original material included in this paper. The authors also confirm that they have obtained permission, from the copyright holder of any third party material included in this paper, to publish it as part of their paper. The authors confirm that they give permission, or have obtained permission from the copyright holder of this paper, for the publication and distribution of this paper as part of the IFASD-2019 proceedings or as individual off-prints from the proceedings.


Marie-Odile Soyer-Gobillard

Ph.D Françoise Paris

Ph.D Laura Gaspari

M.D Philippe Courtet

Ph. D Charles Sultan

Association between fetal DES-exposure and psychiatric disorders in adolescence/adulthood: evidence from a French cohort of 1002 prenatally exposed children

Keywords: Prenatal DES exposure, Psychiatric disorders, Suicide, Adolescence, Adulthood

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

From the 1950s to the late 1970s, millions of women took diethylstilbestrol (DES) during pregnancy to prevent miscarriages and premature births, to inhibit the milk inflow after childbirth, or to treat infertility, primary or secondary amenorrhea, dysmenorrhea and other genital disorders [START_REF] Smith | Diethylstilbestrol in the prevention and treatment of complications of pregnancy[END_REF]. Strikingly, DES continued to be used despite various alerts that were published as early as 1940 [START_REF] Smith | Diethylstilbestrol in the prevention and treatment of complications of pregnancy[END_REF] and the demonstration of Diekmann et al. [START_REF] Dieckman | Does the administration of diethylstilbestrol during pregnancy have therapeutic value?[END_REF] that DES was ineffective in preventing miscarriage or premature birth. In 1971, it emerged that in utero DES exposure was associated with somatic effects in adulthood, including female genital abnormalities [START_REF] Giusti | Diethylstilbestrol revisited: A review of the long term health effects[END_REF][START_REF] Newbold | Lessons learned from perinatal exposure to diethylstilbestrol[END_REF], vaginal cancer [START_REF] Herbst | Adenocarcinoma of the vagina. Association of maternal stilboestrol therapy with tumor appearance in young women[END_REF], and male urogenital disorders [START_REF] Gill | Association of diethylstilbestrol exposure in utero with cryptorchidism, testicular hypoplasia and semen abnormalities[END_REF]. In experimental studies, mice exposed to DES during gestation were found to be more aggressive than controls [START_REF] Palanza | Prenatal exposure to endocrine disrupting chemicals: effects on behavioral development[END_REF], and the offspring of rats treated by intraperitoneal administration of another synthetic estrogen, ethinylestradiol (EE), were shown to present altered emotional states like anxiety and depressive-like behaviors [START_REF] Dugard | Prenatal exposure to ethinylestradiol elicits behavioural abnormalities in the rat[END_REF][START_REF] Arabo | Administration of 17-alpha-ethinylestradiol during pregnancy elicits modifications of maternal behaviour and emotional alteration of the offspring in the rat[END_REF]. In humans, only three large epidemiological studies [START_REF] Vessey | A randomized double-blind controlled trial of the value of stilboestrol therapy in pregnancy: long term follow-up of mothers and their offspring[END_REF][START_REF] Verdoux | Serious psychiatric outcome of subjects prenatally exposed to diethylstilbestrol in the E3N cohort study[END_REF][START_REF] O'reilly | Diethylstilbestrol exposure in utero and depression in women[END_REF] and seven studies on a small number of exposed subjects (as summarized in the critical review of Kebir and Krebs) [START_REF] Kebir | Diethylstilbestrol and risk of psychiatric disorders: A critical review and new insights[END_REF] have pointed out a risk of psychological or psychiatric disorders during adolescence or post-adolescence after in utero DES exposure. However, the data are discrepant [START_REF] Verdoux | Serious psychiatric outcome of subjects prenatally exposed to diethylstilbestrol in the E3N cohort study[END_REF] regarding the impact of DES on the developing fetal brain and the consequences for adolescents and adults. The aim of this work was to determine whether prenatal exposure to DES can induce psychological and/or psychiatric disorders in teenagers and young adults.

Materials and methods: The HHORAGES cohort

The HHORAGES Association (Halt to Synthetic Hormones for Pregnancies), a patient support group, has assembled a cohort of 1280 women who took DES during pregnancies. A detailed questionnaire (Table 1) was prepared and validated by researchers and doctors and sent to these families. As psychiatric diseases generally appear in the post-adolescent period, after 18 years, the questionnaire concerned people born between 1946 and 1995. An authorization request was sent to the French "Commission Nationale de l'Informatique et des Libertés" (CNIL), which authorized the use of the questionnaire. We included items on patient sex, rank among siblings, specific exposure during gestation (DES alone or in association with EE), and somatic and psychiatric disorders.

Out of these 1280 families, 529 families responded to this extensive questionnaire, providing information on the family history; the mother's hormone treatment before and during pregnancy; medical records, including prescriptions; and the health problems of the children, including diagnoses and any medical treatments or hospitalizations. The psychiatric disorders reported in the questionnaires were classified as: eating disorders, behavioral disorders, obsessivecompulsive disorders, depression and bipolar disorders, schizophrenia, and suicides (attempts, death). When suicides or suicide attempts were reported, a second questionnaire validated by members of the Research Group on Suicide (CHU Lapeyronie Montpellier) was sent to the families.

Results

The results are presented in Table 2.

We collected completed questionnaires from 529 families, corresponding to 1182 children. Among these 1182 children, 180 firstborn children without DES exposure served as intrafamilial controls (Table 2). The analysis was conducted according to specific descriptive criteria: treatment by synthetic xenohormones (DES, EE), exposed and unexposed children, patient sex, birth order rank, and somatic and psychiatric disorders.

Of the 1182 children (552 sons + 630 daughters) born from 529 mothers, we noted that 180 of them (102 sons + 78 daughters) were born without any DES exposure, before a DES-exposed pregnancy. We named this group "Pre-DES" (Group 1). 740 children (315 boys + 405 daughters + 20 stillborn) were born after DES exposition, this synthetic estrogen (banned only in 1977 in France) having been prescribed alone or in combination with ethinylestradiol (EE), we called this group "DES exposed" (Group 2). The third group was composed of the 262 children (130 sons + 132 daughters) born from a mother that had taken DES for a previous pregnancy: we named this group "Post-DES" (Group 3) (Table 2). In total, 1002 children were exposed (740 children exposed under prescription plus 262 Post-DES children, i.e., born after a previous pregnancy treated with DES).

In the exposed and post-DES children (Group 2: n=740 -20 stillborns + Group 3: n=262, for a total of 982) we found that 603 children from Group 2 (250 sons + 353 daughters) and 16 from Group 3 were affected by psychiatric disorders (Table 2), in addition to 81 children with genital malformation and/or sterility but no apparent psychiatric disorders. In Group 2, 70 children were without any disorder and in Group 3, 242 children were without disorder.

For the 250 boys, psychiatric disorders were determined on the basis of family testimonies or medical records of antipsychotic or antidepressant prescription and/or psychiatric hospitalization, as follows: behavioral disorders, violence, aggressiveness, and obsessivecompulsive disorders (n=47), eating disorders (n=6), schizophrenia (n=112), and depression, bipolar disorder and anxiety (n=85), with schizophrenia being the most prevalent pathology. For the 353 girls, psychiatric disorders were as follows: behavioral disorders and obsessivecompulsive disorders (n=62), eating disorders (n=75), schizophrenia (n=53), and recurrent depressive or bipolar disorders (n=163), with affective disorders being the most prevalent.

Among the 262 post-DES cases, 16 of them (6.1%) presented psychiatric disorders and only four presented somatic disorders (4 girls). Among the 16 subjects presenting psychiatric disorders, one suffered behavioral disorders, six presented schizophrenia (3 boys + 3 girls), and nine presented bipolar disorders and depression (3 boys + 6 girls); among them, one suicide and five suicide attempt series were noted.

From the 1182 children of our cohort, 180 firstborns were unexposed (102 firstborn boys + 78 firstborn girls): none of these firstborns presented any disorders, had ever made a suicide attempt or had committed suicide. Conversely, among the 982 exposed and post-DES children (1002 minus the 20 stillborns), we noted 33 suicides (24 sons + 9 daughters) and 107 series of suicide attempts (42 sons + 65 daughters), a series comprising between two and 15 suicide attempts per person. Thus, if we assume that the mean number of suicide attempts was six per person, about 642 suicide attempts were made (Table 3). The psychological and/or psychiatric disorders are presented in Table 4 along with the prevalences for these disorders in the general population [START_REF] Nicoli | Ideation et conduites suicidaires en France: prevalence sur la vie et facteurs de risque dans l'étude ESEMeD[END_REF][START_REF] Lépine | Prévalence et comorbidité des troubles psychiatriques dans la population générale française: résultats de l'étude épidémiologique ESEMeD/MHEDEA 2000/(ESEMeD)[END_REF].

Discussion

Very few studies have investigated the impact of prenatal exposure to DES or synthetic hormones on psychiatric outcome. Up to now, most of the studies have been carried out on rodents and have shown the toxic effects of synthetic estrogens on the descendants, particularly in terms of behavioral disorders, as well as other effects like cancers [START_REF] Turusov | Occurrence of tumors in the descendants of CBA male mice prenatally treated with diethylstilbestrol[END_REF], with multigenerational carcinogenic effects in mice [START_REF] Walker | Intensity of multigenerational carcinogenesis from diethylstilbestrol in mice[END_REF]. Injection of 17-alpha-ethinyloestradiol (EE) in pregnant rats induces not only a high rate of abortion, but also anxiety and depression-type disorders in pups [START_REF] Dugard | Prenatal exposure to ethinylestradiol elicits behavioural abnormalities in the rat[END_REF][START_REF] Arabo | Administration of 17-alpha-ethinylestradiol during pregnancy elicits modifications of maternal behaviour and emotional alteration of the offspring in the rat[END_REF]. Based on these findings, it is hypothesized that estrogenic hormones induce neurodevelopmental disturbances in exposed human subjects and may potentially mediate an increased risk of behavioral and psychiatric disorders. Although selection biases (family association) may somewhat affect the impact of our data, the high prevalence of psychiatric disorders in the exposed subjects is quite striking, and it should be noted that the intra-familial controls were completely without disease. Our data therefore clearly demonstrate that DES exposure during pregnancy is associated with a strikingly high incidence of behavioral and/or psychiatric disorders.

The prevalence of eating disorders in our cohort was eight times higher in the exposed adolescents and young women, as previously shown in a study of Gustavson et al., who reported on 1711 women with an OR of 5.72 [START_REF] Gustavson | Increased risk of profound weight loss among women exposed to diethylstilbestrol in utero[END_REF]. Behavioral disorders (11.2%) were also higher in the exposed population. Depression was four times , 26.2%) more prevalent than in the general population, with the DES daughters being more affected than the sons. This observation confirms the epidemiological study from O'Reilly et al. [START_REF] O'reilly | Diethylstilbestrol exposure in utero and depression in women[END_REF]. This group reported an OR of 1.47 in a group of 1612 exposed DES daughters. Our data also agree with the findings of Vessey et al. [START_REF] Vessey | A randomized double-blind controlled trial of the value of stilboestrol therapy in pregnancy: long term follow-up of mothers and their offspring[END_REF], who showed that the offspring of the exposed group had twice as much anxiety as the offspring of the unexposed group. Several small case-control studies have since confirmed these findings [START_REF] Meyer-Bahlburg | Depression n adults with a history of prenatal DES exposure[END_REF][START_REF] Meyer-Bahlburg | A prenatal-hormone hypothesis for depression in adults with history of fetal DES exposure[END_REF][START_REF] Ehrhardt | Psychopathology in prenatally DES-exposed females: current and lifetime adjustment[END_REF][START_REF] Fried-Cassorla | Depression and diethylstilbestrol exposure in women[END_REF][START_REF] Pillard | Psychopathology and social functioning in men prenatally exposed to diethylstilbestrol (DES)[END_REF], all reporting rates of major depressive disorders (MDD) significantly higher than for the general population. Our data illustrate a higher risk of schizophrenia, as this disease was 17 times more prevalent than in the general population, with sons being more affected than DES daughters. To date, there have been few data regarding this association. In 1987, Katz et al. [START_REF] Katz | Psychosis and prenatal exposure to diethylstilbestrol[END_REF] for the first time reported four subjects in a very small cohort suffering from psychosis after DES exposure; in 1993, Pillard et al. [START_REF] Pillard | Psychopathology and social functioning in men prenatally exposed to diethylstilbestrol (DES)[END_REF] observed several subjects in another cohort showing "paranoid ideation" (the diagnostic term today being "paranoid schizophrenia"). More recently, Kebir and Krebs (unpublished data in 13) analyzed the data from a small HHORAGES cohort and noted the high frequency of psychotic disorders but claimed that this finding did not demonstrate any causality between in utero DES exposure and a specific disorder. Nevertheless, these authors also examined three epidemiological studies with worldwide cohorts and several smaller cohorts with the objective of determining the impact of DES on the development of psychiatric disorders, including schizophrenia [START_REF] Kebir | Diethylstilbestrol and risk of psychiatric disorders: A critical review and new insights[END_REF]. They found that out of all these studies, only one did not support the hypothesis of a causal link [START_REF] Verdoux | Serious psychiatric outcome of subjects prenatally exposed to diethylstilbestrol in the E3N cohort study[END_REF]. With regard to suicides, our work clearly demonstrates a drastically increased risk of suicide attempts (65.4% versus 0 in the unexposed controls and 0.25% in general population) and suicides (3.4% versus 0 in the unexposed controls and 0.02% in general population). It could be noted that, as in the general population, DES sons commit more suicides than DES daughters, and inversely for suicide attempts. Moreover, our data reveal that 50% of the sons who committed suicide suffered from schizophrenia. Psychiatric studies in general have shown that the percentage of suicides is generally higher in individuals with psychiatric disorders than in the general population. But to our knowledge, there is no information or specific studies concerning this association in the context of DES exposure.

Sixteen subjects in Group 3 (post-DES children) (Table 2) had diagnosed psychiatric disorders.

An explanation for this finding might be that DES , being very lipophilic synthetic estrogens, remains in the mothers' fat after estrogenic impregnation in a previous pregnancy and are then released through the placental barrier in the next gestation.

Conclusion

In conclusion, we demonstrate here that prenatal exposure to DES is associated with psychiatric disorders in adolescence and adulthood. These data have been supported by recent preliminary work by Réseau DES-France suggesting a higher risk of psychiatric disorders in DES daughters, based on an analysis of the data from 3436 DES daughters versus 3256 control women [START_REF]Distilbène, Etude sur trois générations. Résultats préliminaires[END_REF].

Moreover, in a recent review on the risk of psychiatric disorders in adolescents exposed in utero to DES, Kebir and Krebs [START_REF] Kebir | Diethylstilbestrol and risk of psychiatric disorders: A critical review and new insights[END_REF], in agreement with Abdomalecky et al. [START_REF] Abdolmaleky | Methylomics in psychiatry: modulation of gene-environment interactions may be through DNA methylation[END_REF], noted that the role of prenatal exposure to DES as an environmental factor for psychiatric disorders may involve a gene-environment interaction, as well as an epigenetic mechanism. This latter group reported for the first time that prenatal exposure to DES induced specific modifications in fetal Relational difficulties: in married life, in professional life?

Children, grandchildren: How many? Full term? Health condition? Malformations or observed disorders? Table 1. Questionnaire sent to families.

  : Surname, First Name, Date of Birth, Address, Home phone number, Professional phone number, Mobile number, E-mail address Family Situation, Number of children, children's first names and dates of birth. Professional situation. Pill taken: before first pregnancy, and between pregnancies Hormonal treatment before pregnancy (or pregnancies). (If a treatment was prescribed, please indicate the nature of the treatment and the time elapsed between the end of treatment and the beginning of the pregnancy) Miscarriage(s): Indicate the time when it occurred, in relation to the other full term pregnancies, if any (before the first one, or between two full term pregnancies? Professional Exposure to hormones, to chemicals (pesticides, etc...) Psychiatric or psychopathological family history (father, mother, .). Health problems after first child birth Pregnancy and childbirth (For more than 2 children, T.O.P.) 1st child, 2 nd child.. Treatment during pregnancy: (Yes/No) -Nature of medicines -Doses -What time during pregnancy (first and last day of treatment, expressed in weeks from the beginning of pregnancy)? What medical reason was put forward (possible miscarriage, comfort, etc...)? At what month did the delivery occur, from the start of the pregnancy? General condition of the child, Sex, Weight, Health problems at birth (mother and child) Existing documents (prescriptions, medical files, etc...) and testimoniesWould you be so kind as to send us the copies of the files and documents? If yes, in order to save time, please attach the copies of all relative documents in your possession.Health problems of your children(For more than 2 children, T. OP) 1st child, 2nd child First name, Sex, Birth Date, Rank in the sibling order Physical disorders: Which ones? What age? Sterility treatment? Surgery? Psychological disorders: Nature of first symptoms? What age? Subsequent aggravation? And at what age? Other data: Hospitalization, Violence, Suicide attempt(s), Medical Treatments, Diagnosis, AAH, Disability?

Table 1. Questionnaire sent to families

Acknowledgments

Mr René Alexandre, Dr Henri Pézerat and Prof. Jean Caston have worked tirelessly to demonstrate the harmful role of endocrine disruptors, particularly synthetic hormones, before passing away: we shall never forget their great help. This work could not continue without the daily ongoing support of the HHORAGES board, especially Geneviève Alchourroun,

DNA methylation that could be involved in postnatal neurodevelopment and the subsequent psychiatric disorders in adolescence (27).

Authors

Marie-Odile Soyer-Gobillard (mog66@orange.fr) Centre National de la Recherche Scientifique, Unité Mixte de Recherche 7628 Université P. et M. Curie (Paris6) Laboratoire Arago, F-66650 Banyuls sur mer and Association HHORAGES, (Halte aux HORmones Artificielles pour les GrossessES), Le Prieuré de Baillon, F-95270 Asnières sur Oise.www.hhorages.com Charles Sultan (c-sultan@chu-montpellier.fr) Françoise Paris (f.paris@chu-montpellier.fr) Laura Gaspari (dr.lauragaspari@gmail.com) Service d'Hormonologie (Développement et Reproduction), Hôpital Lapeyronie, CHU de Montpellier, Montpellier, France et Institut de Génétique Humaine, CNRS UPR 1142, Montpellier, France Service de Pédiatrie I, Unité d'endocrinologie et gynécologie pédiatrique, Hôpital Arnaud-de-Villeneuve, CHU de Montpellier, Montpellier Cedex 5, France Philippe Courtet (p-courtet@chu-montpellier.fr) Département d'Urgence et Post Urgence Psychiatrique, CHU Montpellier & Université Montpellier, France

Disclaimer

The authors declare they have no competitive financial interest concerning this paper; HHORAGES-France Association is financed exclusively by donations.

Mauricette Puillandre, Denise Hemmerdinger, Michel Datry, Yette Blanchet and Charles Zelwer. We are very grateful to Professor MO Krebs and Dr O Kébir for believing in our hypothesis nine years ago and for working so hard to scientifically demonstrate it. Table 2. Classification of the three groups of the HHORAGES cohort.

Table 3 :

Table 3 : Total number of psychological/psychiatric disorders among the 982 (1002-20 stillborn) DESexposed and post-DES children.