Timothée Achard
email: timothee.achard@onera.fr

Christophe Blondeau
email: christophe.blondeau@onera.fr

Roger Ohayon

COMPARISON OF HIGH-FIDELITY AERO-STRUCTURE GRADIENT COMPUTATION TECHNIQUES. APPLICATION ON THE CRM WING DESIGN

Keywords: Sensitivity analysis, Multidisciplinary Optimization, Adjoint method

Aero-structural optimization is a keystone to concurrently improve aerodynamic performance and reduce the structural mass of an aircraft. Gradient-based multi-disciplinary design optimization is actually efficient if gradients computations are fast and accurate enough. This paper presents two high-fidelity aero-structure gradient computation techniques for design variables impacting the stiffness of the structure. A first part details the new module developed and implemented in the Onera elsA software capable of computing aero-structure gradients using intrusive direct and adjoint methods. In a second part an alternative improved uncoupled non-intrusive approach is proposed. Both approaches are evaluated and compared on the Common Research Model (CRM) test-case considering criteria such as accuracy, efficiency, and applicability on practical industrial problem. Finally, preliminary results are presented on an inverse design problem. The final objective will be to match a target twist distribution on the Common Research Model wing using this new coupled gradient capability.

INTRODUCTION

The raise of composite materials in modern aircraft structures increases the structural flexibility of wings [START_REF] Weisshaarr | Divergence of forward swept composite wings[END_REF]. Consequently, fluid-structure interaction has to be considered when designing such flexible wings. With increasing computational resources, static aeroelastic simulations based on high-fidelity models are now commonly accessible and capable of tackling complex aeroelastic couplings. To improve the design of flexible wings, gradient based multi-disciplinary design optimization techniques are now effective and widely used [START_REF] Kenway | Aerostructural optimization of the common research model configuration[END_REF][START_REF] Jameson | Aerodynamic-structural design studies of low-sweep transonic wings[END_REF][START_REF] Brezillon | Development and application of multi-disciplinary optimization capabilities based on high-fidelity methods[END_REF]. However, gradients calculation is not trivial and can be costly when high-fidelity models are considered. Indeed, in an integrated optimization process, this difficulty increases for required cross-sensitivity derivatives such as the derivatives of aerodynamic performance with respect to structural sizes and derivatives of structural response with respect to changes in aerodynamic shape [START_REF] Grossman | Integrated aerodynamic/structural design of a sailplane wing[END_REF].

Several developments have been made over the last decade to compute aero-structure gradients using adjoint method [START_REF] Ghazlane | Aerostructural adjoint method for flexible wing optimization[END_REF][START_REF] Kenway | Scalable parallel approach for high-fidelity steady-state aeroelastic analysis and adjoint derivative computations[END_REF][START_REF] Zhang | High-fidelity aerostructural optimization with integrated geometry parameterization and mesh movement[END_REF]. At ONERA, active research work dedicated to aero-structure optimization was achieved during the past decade. In 2008, Marcelet [START_REF] Marcelet | Sensitivity analysis of a strongly coupled system using the discrete direct and adjoint approach[END_REF] developed a module capable of computing static aeroelastic equilibrium and sensitivity analysis using adjoint method. However, the structural model behavior was limited to an equivalent beam model kinematics, and the aeroelastic gradients were only computed with respect to aerodynamic shape parameters. In a follow-on work, Ghazlane extended this work to take into account structural parameters (typically stiffener cross-section, plate thickness...) for aero-structure gradients computation [START_REF] Ghazlane | Aerostructural adjoint method for flexible wing optimization[END_REF]. Lastly, Viti [START_REF] Viti | Aero-structural approach coupled with direct operative cost optimization for new aircraft concept in preliminary design[END_REF] demonstrated a bi-step aero-structure preliminary design of a forward-swept wing applying these gradient computation tools. Nonetheless, the structural model was still limited to a beam model in all these works. Although useful for preliminary design studies and design space exploration, this structural modeling is not able to tackle complex aeroelastic couplings, nor deal with composite structures for efficient aeroelastic tailoring design, and is obviously limited to large aspect ratio wings.

The objective of a previous work [START_REF] Blondeau | Recent achievements towards aero-structure gradient computation using high-fidelity cfd-csm in the onera elsa software[END_REF] by the authors was to relieve all these limitations by re-developing an aero-structure gradient capability from scratch, by systematic hand differentiation of discretized aero-elastic equations, and all associated operators, in the elsA aeroelasticity module (elsA/Ael) [START_REF] Girodroux-Lavigne | Recent navier-stokes aeroelastic simulations using the elsa code for aircraft applications[END_REF]. This way, the new capability inherits the structural paradigm embedded into the elsA/Ael module (structural flexibility matrix), the whole catalog of fluid-structure transfer methods, as well as the underlying parallel architecture. The equations and the solving procedure for the aeroelastic and the aero-structure gradients have been discussed in [START_REF] Blondeau | Recent achievements towards aero-structure gradient computation using high-fidelity cfd-csm in the onera elsa software[END_REF], where only aeroelastic gradients computation capability (i.e. gradients of aeroelastic functions of interest with respect to aerodynamic shape parameters) was available. This module, named elsA/Aoc, has been recently extended in order to take into account structural design parameters in the direct and the adjoint mode.

The efficiency of the adjoint approach holds when a reduced number of responses is considered compared to a large set of design parameters. However, an industrial wing box sizing for preliminary design considers up to several hundreds of design parameters and up to several tens of thousands of structural constraints. In order to still benefit from the elegant adjoint formulation, some authors have proposed a constraint aggregation approach [START_REF] Akgn | Efficient structural optimization for multiple load cases using adjoint sensitivities[END_REF][START_REF] Poon | An adaptive approach to constraint aggregation using adjoint sensitivity analysis[END_REF] as an attempt to reduce the size of the set of admissible constraints, thus the number of gradients to compute. However, this technique has several drawbacks and may lead to sub-optimal designs. Considering these difficulties, the authors have proposed in a previous paper [START_REF] Achard | An uncoupled approach to compute aero-structure gradients using high-fidelity cfd-csm[END_REF] a new strategy to compute high-fidelity aero-structure gradients (i.e. gradients of any function of interest with respect to structural design parameters). This strategy is based on an uncoupled non-intrusive approach benefiting from the linearized aerodynamic formulation. The main advantage of the method is the independency of the computational cost with respect to the number of constraints and potentially with respect to the number of structural design parameters. In addition this approach only requires rigid CFD computations.

The final objective of this work is to compare the efficiency of the uncoupled approach for aerostructure gradients computation, with the classical intrusive tangent or adjoint method, through the Common Research Model (CRM) wing design. The comparison aim at addressing the applicability of the optimization strategy on realistic industrial applications. Both approaches are investigated on the CRM test-case and criteria such as accuracy, efficiency and applicability on practical industrial problem are discussed.

The remainder of this paper is organized as follows: we first detail the intrusive direct and adjoint approaches. Then, an improved version of the non-intrusive method originally proposed in [START_REF] Achard | An uncoupled approach to compute aero-structure gradients using high-fidelity cfd-csm[END_REF] is presented. Next, aero-structure gradients are computed on the CRM using these approaches. Finally, a strategy is proposed for an efficient structural sizing using high-fidelity aerostructure gradients, and preliminary results on the CRM wing are discussed.

The first theoretical parts presenting the two gradient computation methods have already been detailed in [START_REF] Achard | High-fidelity aero-structure gradient computation techniques. application to the onera m6 wing[END_REF]. They are presented again to preserve the global consistency of this paper.

INTRUSIVE APPROACH

Let's assume that a steady aeroelastic equilibrium is available. In the following, we will denote the corresponding mesh as "flight shape", as opposed to the unloaded reference mesh called "jig shape". This is illustrated in Figure 1, where X a is the equilibrium aerodynamic mesh, X s the structural mesh, and X a0 the aerodynamic reference mesh. Our objective is to compute the sensitivity of a function of interest J with respect to a set of structural design variables p affecting the stiffness matrix only. To this end, direct and adjoint approaches derived from the linearization of the aeroelastic equations around the equilibrium position will be used.

R a (X a , W, U) = 0 R s (X s , W, U) = 0 (1)
with R a and R s a nonlinear and a linear system of discrete equations. These two blocks of equations are coupled through aerodynamic forces Q a loading the structure skin and the structural displacements U deforming the fluid mesh. The structural loads Q s are obtained with a suitable load transfer technique applied to Q a :

Q s (Q a (W, X a ), X a0 , X s ) = T Q surf (X a0 , X s )Q a (W, X a ) (2)
where T Q surf (X a0 , X s ) represents the linear load transfer operator. The structural displacements alter the fluid grid positions through the relation:

X a = X a0 + δX a (δX a,surf , X a0 ) = X a0 + T vol (X a0 )δX a,surf (3)
with T vol (X a0 ) the volume operator performing the deformation of the fluid interior domain.

The vector δX a,surf corresponds to the displacements of the fluid nodes at the aeroelastic interface:

δX a,surf = δX a,surf (X a0 , X s , U) = T U surf (X a0 , X s )U (4)
where T U surf (X a0 , X s ) represents the linear displacement transfer operator.

Direct approach

Let's consider the scalar aeroelastic objective function J(W, X a , X s ) and a structural parameter p. We assume that the implicit function theorem is applicable to our problem. In this case, the state variables W and U can be expressed as implicit functions W(p) and U(p) of the optimization parameters and have the same properties of regularity than R a and R s . Direct differentiation of system 1 with respect to p yields the tangent problem for the aero-structure gradient:

     ∂R a ∂W ∂R a ∂U ∂R s ∂W ∂R s ∂U           dW dp dU dp      =      - ∂R a ∂X a dX a dp - ∂R s ∂X s dX s dp      (5)
The resolution of this system makes the computation of the gradient of the objective function J possible:

dJ dp = ∂J ∂W dW dp + ∂J ∂X a dX a dp + ∂J ∂X s dX s dp (6)
Our objective is now to expand equations 5 and 6 only with respect to the unknown vectors dW/dp et dU/dp. First, we remind that p affects only the stiffness of the structure, so that the following relations hold:

∂X a0 ∂p = ∂X s ∂p = 0 (7)
We then differentiate equation 3 with respect to p, and using equation 4, we get:

dX a dp = dδX a dp = T vol T U surf dU dp (8)
Substituting these two equations into 6 yields the expression of the total gradient of the objective function:

dJ dp = ∂J ∂W dW dp + ∂J ∂X a T vol T U surf dU dp (9)
Terms ∂J/∂X a and ∂J/∂W only depend on the equilibrium steady state and are calculated analytically in a pre-processing step.

The structural residual can be written K(p)U(p) -Q s (p) = 0, where K is the stiffness matrix, and Q s the structural loads. Using the functional dependency of the structural loads (see equation 2), the structural block reads:

K dU dp = T Q surf ∂Q a ∂W dW dp + ∂Q a ∂X a dX a dp - ∂K ∂p U (10)
The resolution of the coupled system 5 is done in an iterative way, alternatively on each block of equations. Using this last relation, system 5 can be cast into the compact form:

             ∂R a ∂W dW dp (k+1) = - ∂R a ∂X a [A] dU dp (k) K dU dp (k+1) = [B] dW dp (k) + [C] dU dp (k) -[D]U (11)
Constant matrices [A],[B],[C] and [D] are defined analytically with the following formulas

[A] = T vol T U surf [B] = T Q surf ∂Q a ∂W [C] = [M][A] [D] = ∂K ∂p (12)
where the utility matrix

[M] = T Q surf ∂Q a ∂X a has been introduced.
The system in equation 11 is presented in a form suitable for an iterative block scheme resolution similar to the lagged-block strategy formely proposed in [START_REF] Maute | Coupled analytical sensitivity analysis and optimization of three-dimensional nonlinear aeroelastic systems[END_REF], superscript k being the current iteration number. Following the same notations, equation 9 becomes:

dJ dp = ∂J ∂W dW dp + ∂J ∂X a [A] dU dp (13)

Adjoint approach

In order to obtain the adjoint equations of the aero-structure gradient, the objective function is formulated by adding the total variation of the residuals R a and R s to dJ/dp. These variations with respect to the design variables being null, we then write for any vector Λ a and Λ s :

dJ dp = ∂J ∂W dW dp + ∂J ∂X a [A] dU dp + Λ T a dR a dp + Λ T s dR s dp (14)
Expanding dR a /dp and dR s /dp and factorizing out the difficult terms dW/dp and dU/dp yields:

dJ dp = ∂J ∂W + Λ T a ∂R a ∂W -Λ T s [B] dW dp + ∂J ∂X a [A] + Λ T a ∂R a ∂X a [A] + Λ T s (K -[C]) dU dp + Λ T s [D]U (15)
where adjoint vectors Λ a and Λ s are chosen such that the following system is satisfied:

      ∂R a ∂W T -[B] T [A] T ∂R a ∂X a T K T -[C] T         Λ a Λ s   =       - ∂J ∂W T -[A] T ∂J ∂X a T       (16)
Finally, once adjoint unknowns are determined, the total gradient is obtained with:

dJ dp = Λ T s [D]U = Λ T s ∂K ∂p U (17)
Direct and adjoint equations have been written here for the particular case of a scalar aeroelastic function of interest (lift or drag coefficient), with optimization parameters controlling only the stiffness of the structure (e.g. stiffener cross-section, plate thickness). In a previous work by the authors [START_REF] Blondeau | Recent achievements towards aero-structure gradient computation using high-fidelity cfd-csm in the onera elsa software[END_REF], these equations have been detailed for the more general case of parameters affecting both the aerodynamic and the structure meshes. It is very interesting to note that the adjoint system 16 is exactly the same as the one presented in [START_REF] Blondeau | Recent achievements towards aero-structure gradient computation using high-fidelity cfd-csm in the onera elsa software[END_REF]. This means that this adjoint system is independent of the nature of the design parameter. Consequently, one resolution of this system gives access to the gradient of the aeroelastic function of interest with respect to all the parameters, regardless if they control aerodynamic shape, structural properties or even affect simultaneously stiffness and shape. Indeed, once the system is solved, the computational cost to obtain the gradient of the objective function with respect to these parameters is negligible. This constitutes a major advantage in the simultaneous optimization of both aerodynamic and structure.

NON-INTRUSIVE APPROACH

The non-intrusive technique allows a complete decoupling of the structure and the aerodynamic analysis. The idea consists in building a projection of the static displacement field on the structural modeshapes in order to re-use the rigid linearized kernel of the module elsA/Opt. This approach has been discussed in [START_REF] Achard | An uncoupled approach to compute aero-structure gradients using high-fidelity cfd-csm[END_REF]. With this strategy, aero-structure gradients have been computed on the ONERA M6 wing and then on the CRM configuration with satisfactory results. Both configurations used a RANS aerodynamic fluid model and a tridimensional FEM. In this section, after a brief reminder of the theory, some improvements are proposed to reach a better accuracy for the gradients.

Working process with the linearized solver

The linearized solver determines the perturbed field of conservative variables δW with respect to the steady state flow, when the fluid domain boundary is perturbed by a prescribed displacement induced by a perturbation of a structural design parameter δp. The corresponding perturbed structural loads δQ s are then deduced directly from the solution δW. This process is illustrated in Figure 2. Flexibility is taken into account into the steady solution and the computation of δW is performed by a linearization around the equilibrium mesh position denoted Flight shape.

Let's consider the tangent system 11. Solving for the first block of equations and substituting dW/dp in the structure block yields:

           K -[Z] dU dp = - ∂K ∂p U [Z] = -[B] ∂R a ∂W -1 ∂R a ∂X a [A] + [D] (18)
K -[Z] appearing in the first block of system 18 is the Schur complement of the matrix block ∂R a /∂W. We remind that direct differentiation of the static equation KU = Q s with respect to p gives:

∂K ∂p U + K dU dp = dQ s dp (19)
Comparing this expression with the first block of equation in system 18 yields the relation:

dQ s dp = [Z] dU dp (20)
where [Z] is the kernel matrix of the linearized flow solver.

Modal reconstruction of aero-structure gradients from an aeroelastic configuration

Let's consider the steady equilibrium corresponding to a fluid state W, a mesh position X a and a structural displacement field U:

KU = Q s = [T Q surf ]Q a (W, X a ) (21)
K is the stiffness matrix of the finite element model, Q s the nodal structural loads, Q a the aerodynamic loads computed on the aerodynamic skin. If the full set Φ of structural eigenvectors was available, the modal projection U = Φq would hold exactly. Then substituting in equation 21 gives the corresponding generalized coordinates as:

q = γ -1 Φ T Q s ( 22
)
where γ = Φ T KΦ denotes the diagonal generalized stiffness matrix. In practice only a restricted set of first n Φ eigensolutions is computed and the modal approximation to U becomes:

U ≈ U Φ = n Φ i=1 Φ i q i (23)
However, remind that U is computed exactly from equation 21 such that the residual term in equation 23 is known from the simple difference U res = U -U Φ . Equation 21 can now be reformulated as

K(U Φ + U res ) = Q s (24)
Inserting the modal decomposition of U in the expression of the gradient of structural loads (equation 20) leads to

dQ s dp = Z(Φ dq dp + ∂Φ ∂p q) + Z ∂U res ∂p (25)
The use of partial derivatives in the equation above means that Φ and U res only depend on p. This latter assumption is valid if U res is considered as a static residual mode (i.e. a structural displacement under an assumed prescribed load case). In order to exploit this relation, the gradient of the generalized coordinates has to be determined. With some algebra manipulation, it is possible to demonstrate that:

dq dp = [γ -GAF] -1 - ∂γ ∂p q + Φ T Z ∂Φ ∂p q + ∂U res ∂p (26)
where the generalized aerodynamic forces matrix is defined as GAF = Φ T ZΦ.

The computational cost of this method can be easily determined with a reformulation of equation 25:

dQ s dp = ZΦ dq dp + Z( ∂Φ ∂p q + ∂U res ∂p ) (27)
The first part (ZΦ) of equation 27 costs one linearized computation per mode. The second part (Z( ∂Φ ∂p q + ∂Ures ∂p )) requires one linearized computation per parameter. Once these two parts are determined, the gradient of the generalized coordinates is readily available.

Usually, few mode shapes are necessary to correctly approximate the structural displacements (equation 24), however many design variables can be defined in the optimization problem. In this case, the method loses interest since it is independent of the number of objective functions, but dependent on the number of design variables. In an attempt to reduce the computational cost, gradients can be approximated by dropping the terms which depend on design variable number. In this case equations 26 and 27 become:

dq dp = [γ -GAF] -1 - ∂γ ∂p q (28) dQ s dp = ZΦ dq dp (29)

CRM TEST-CASE

In this section, we apply the intrusive and non-intrusive approaches to the CRM. Our objective is to compute lift and pressure drag coefficient sensitivities with respect to structural parameters, as well as gradients of structural loads. The aerodynamic skin of the selected wing/body/HTP (WBH) configuration is depicted in Figure 3. For all subsequent results, a RANS fluid model is used with an upwind Roe scheme and a MUSCL interpolation associated to a Van Albada limiting function. The Spalart-Allmaras one-equation turbulence model is selected. In the attempt to find a lift coefficient of 0.5 for a Mach number of 0.85 at 35000 ft, an angle of attack of 2.127 is applied.

Several finite element models of the CRM wing are available on the NASA CRM website. All these models correspond to the flight shape of the CRM. The model V14 exhibits an expected static behavior, with a realistic vertical displacement of 2.15 m observed at the wing tip in nominal cruise conditions. It is a simple finite element model of the wing box and centre box with spars, skins, ribs, implicit stiffeners and a constant wing skin thickness of 8.89 mm. This model can be considered as a good starting point for a structural sizing process. The provided CRM configuration corresponds to the 1 g cruise shape. In order to find the corresponding jig shape, an inverse procedure has been set up and presented in [START_REF] Achard | An uncoupled approach to compute aero-structure gradients using high-fidelity cfd-csm[END_REF]. To correctly predict the structural displacement due to gravity, a realistic distribution of concentrated masses representing the fuel in the wing is added, along with concentrated masses for the engine, pylon and landing gear. Figure 3 presents the finite element model colored by the structural optimization groups, along with the aerodynamic skin of the CRM. Ten structural parameters are defined, split into two groups of five parameters driving the thickness of the lower skin and the upper skin respectively. In the following, the results will be presented only for the lower skin (see Figure 4), since for this exercise upper and lower skin have the same thickness distribution.

The first nine eigenmodes computed on the initial FEM of the CRM wingbox and smoothed on the aerodynamic skin are presented in Figure 5. As detailed in [START_REF] Achard | An uncoupled approach to compute aero-structure gradients using high-fidelity cfd-csm[END_REF], twenty modes will be used to compute the gradients with the non-intrusive approach.

For the intrusive approach, the convergence of the discrete residual norm of the fluid block is reported in Figure 6. In direct mode, five right hand sides are considered, i.e. one for each design parameter. In adjoint mode, two right hand sides are considered, that is one for each objective function.

Tables 1 and2 summarize the total lift and drag coefficient derivatives respectively, and allows a comparison between intrusive and non-intrusive approach. Gradients obtained with the tangent approach, taken as reference here, have been validated with finite difference results on an other test-case in [START_REF] Achard | High-fidelity aero-structure gradient computation techniques. application to the onera m6 wing[END_REF]. The very low discrepancies observed between tangent and adjoint approach demonstrate that the duality between the two methods is numerically satisfied. These results validate for this test-case the implementation of the direct and adjoint method into elsA for gradients of aeroelastic function of interest with respect to structural design variables. The non-intrusive approach with and without the eigenvector derivatives contribution (∂Φ/∂p) gives satisfactory results, since the sign and the order of magnitude are always well predicted compared to the intrusive approach. We remind that without this contribution, the method is independent of the number of design variables, reducing dramatically the computational cost of the gradients.

For the particular case of aeroelastic coefficient derivatives with respect to structural parameters, there are usually few objective functions but many design variables. Therefore, the adjoint method is the preferable approach, since it requires to solve the adjoint system 16 only once per objective function.

Now considering the total derivatives of the structural loads, the computational cost of the adjoint approach is indexed to the number of individual force components. Usually, a realistic structural loading consists in up to several hundreds of individual force components. This is comparable to the typical number of design parameters for a realistic optimization process. In this context, none of the tangent or adjoint method really stands out. However, computing the gradients of structural loads using the non-intrusive approach but without the eigenvector derivatives contribution requires significant lower computational cost, compared to the intrusive approaches. Again, in this case, aero-structure gradients only cost one rigid linearized computation per mode shape.

In order to assess the accuracy of the non-intrusive approach, total derivatives of the structural loads are computed with the tangent approach and compared to the non-intrusive approach (with and without eigenvector derivatives). To have a representative example, gradients of horizontal and transverse forces with respect to p 1 and p 4 are considered. For sake of clarity, normalized gradients are plotted on nodes belonging to the front spar only (see . First of all, results for intrusive and non-intrusive approach qualitatively match pretty well. Taking the intrusive approach as reference, it is observed that the non-intrusive approach is always predictive when eigenvector derivatives contribution is added. Gradients computed with respect to p 1 without eigenvector derivatives are still close to the reference values. However, those computed with respect to p 4 are either inaccurate or wrong. This points out the relevancy of considering the eigenvector derivatives for some design variables. In [START_REF] Achard | An uncoupled approach to compute aero-structure gradients using high-fidelity cfd-csm[END_REF], the authors propose a criterion that measures the amount of relevant information contained in the eigenvector derivatives compared to the information contained in the modal basis. It reveals whether a specific parameter requires the addition of eigenvector derivatives in the load sensitivity reconstruction. Nevertheless, it can be observed that these discrepancies occur mainly for low-valued gradients, which do not significantly affects gradient-based optimization since the optimizer's path is mostly driven by high-valued sensitivities.

STRUCTURAL SIZING OF THE CRM WING USING HIGH-FIDELITY AERO-STRUCTURE GRADIENTS

Force approximation method

High-fidelity structural sizing requires derivatives of structural responses (typically constraints) with respect to structural parameters. These derivatives are costly, even when computed with the adjoint method since it requires as many resolution of the adjoint system as the total number of constraints in the whole structure. In an attempt to reduce the size of the set of admissible constraints, thus the number of gradients to compute, some authors have proposed to use a constraint aggregation approach in conjunction with the adjoint method. As already mentioned, this approach has several drawbacks and needs to be tuned appropriately.

Approximating constraint responses is not easy, since they may exhibit irregular distribution in the structure. In contrast, the structural loads are easier to approximate due to their smoother distribution. It is therefore proposed to use the force approximation method [START_REF] Vanderplaats | New approximation method for stress constraints in structural synthesis[END_REF] as an alternative to a constraint aggregation approach. Assuming a general linear elastic behavior, the relationship between stresses and strains would be linear of the form σ = H(ǫǫ 0 ) + σ 0 , where H denotes the elasticity matrix and ǫ 0 and σ 0 the initial strain and stress fields. The deformation field is obtained from the nodal displacements through the usual relation ǫ = SNU(p) = BU(p), S being a linear differential operator and N the matrix of shape functions. Direct differentiation of σ with respect to a scalar structural parameter p yields:

dσ dp = ∂σ ∂U dU dp = HB dU dp (30)
Deriving the static equation KU = Q s with respect to p gives the expression of the structural displacement derivatives:

dU dp = -K -1 ∂K ∂p U + K -1 dQ s dp (31)
Equation 30 becomes:

dσ dp = HBK -1 - ∂K ∂p U + dQ s dp = dσ U dp + dσ Q dp (32)
where dσ U /dp is a constraint response under a standard static pseudo-load F U = -(∂K/∂p)U, and dσ Q /dp is a constraint response under the pseudo-load F Q = dQ s /dp. Therefore, for each static load Q s , a companion load dQ s /dp must be defined. With this strategy, the total gradient of a constraint response can be recombined from two separate standard static sensitivity analyses. However, we adopt an even simpler approach by directly providing a linear approximation of the aeroelastic loads to the structure optimizer. Given a structural parameter p varying around its current value p 0 , and knowing the load vector total derivative, we simply formulate the linear approximation:

Q s (p) ≈ Q s (p 0 ) + dQ s dp δp (33)
Then by linearity assumption the total gradient of the constraints is directly obtained, leaving the computational burden to the structure solver.

Proposed strategy for an efficient structural sizing

Given a Finite Element Model (FEM) of the CRM wing and the associated fluid mesh, our objective is to match a target twist distribution extracted from the original geometry of the CRM wing (see [START_REF] Vassberg | Development of a common research model for applied cfd validation studies[END_REF]). This inverse design problem can be formulated as:

Minimize J(p) such that:

                     J(p) = N i=1 U T i -U i (p) U T i 2 ǫ(p) ≤ ǫ max p min ≤ p ≤ p max (34)
where U T is a set of target displacements leading the initial structural mesh to the desired shape. The structure is assumed to be manufactured from Aluminum 2024. The corresponding limit allowable stress is σ max = 324M P a, from which the corresponding maximum allowable strain ǫ max is derived. The set of structural parameters p is defined as the thickness of 27 structural groups depicted in Figure 11. Eleven parameters control the lower and upper skins, ten parameters the front and rear spars, and five parameters the stringers. The optimization process depicted in Figure 12 relies on the aerodynamic software elsA for the computation of aeroelastic loads and associated gradients, and on the MSC NASTRAN software for the structural analysis and optimization part. The process has two levels of convergence: one at the upper level on the load cases and one at the optimizer level on physical properties and design objective function. The structural optimization process is converged for a given set of external non-linear aeroelastic loads. In this optimization cycle the complete FEM is used and the optimization parameters are the physical properties of various structural members.

Then, an equivalent reduced order model (flexibility matrix) of the sized structural model is automatically derived in order to perform the adequate flexible static aeroelastic computations at the current iteration. Once the new flexible non-linear loads are obtained, they are transferred to the structural grid and a new optimization process is performed. This organization is very flexible as it does not require any modification to the NASTRAN core solver.

Figure 12 presents different optimization processes to solve problem 34, according to the strategy chosen to compute the loads. If they come from Doublet Lattice Method (DLM) only (first case), then the optimization strategy merely relies on standard NASTRAN capabilities. The second case named "DLM+Correction" corresponds to the DLM prediction corrected by data from rigid CFD, thus requiring off-line elsA aerodynamic computations. The third case is based on an external CFD based Aeroelastic Analysis (A.A.) to estimate the loads. Thus, at each external iteration, the flexibility matrix S (which ensures the fluid/structure coupling such that U = SQ s ) and the displacements due to gravity have to be updated. These aerodynamic loads are kept constant during the internal optimization loop. Finally, the last case extends the previous one with the contribution of load sensitivities computed with elsA (Aeroelastic + Sensitivity Analysis). The optimizer benefits from the additional contribution of the linearized loads using the force approximation method.

In the first two cases, a strong coupling exists between the aerodynamic and the structural discipline, but the fidelity of the fluid model is low. For the last two cases, of higher-fidelity, there is a weak coupling between these two disciplines, since the loads and their sensitivities are not updated at each re-analysis of the FEM.

A strategy based on a strong coupling and high-fidelity models would involve the merge between the re-analysis loop and the external loads loop. Consequently, for each re-analysis asked by the structural optimizer, elsA would be invoked to compute the loads and their sensitivities for the current set of structural parameters and to provide them to the NASTRAN optimizer. In scheme 12, this means to replace the 20 re-analysis per external iteration by only one re-analysis per external iteration. However, this strategy is expected to be of higher computational cost. As explained in section 4, the non-intrusive approach provides an efficient way to access gradients of loads. However, for cross-gradients like aerodynamic coefficient derivatives with respect to structural parameters, the adjoint method is preferable. With both capabilities readily available, an aerodynamic performance monitoring can be achieved during the sizing process, particularly by taking into account dC d /dp or dC l /dp into the optimization.

Last but not least, a strong advantage of treating the sizing process with separate aerodynamic and structural parts lies in the possibility to use an adapted fluid mesh according to the purpose of the simulation. When gradients of loads are needed, then a Euler mesh might be accurate enough. In contrast, to estimate drag coefficient derivatives, a RANS mesh is more appropriate (see Figure 12). This flexibility for the computation of the different terms is beneficial to save CPU time.

Preliminary results

In this section, first results on the optimization problem exposed above are presented and discussed. The strategy chosen is based on a weak coupling, and the loads are updated with high-fidelity aeroelastic analysis at each external iteration. For these preliminary results, load sensitivities are not taken into account, which corresponds to the third case in Figure 12. This strategy is similar to that presented in [START_REF] Blondeau | A bi-level high fidelity aero-structural integrated design methodology -a focus on the structural sizing process[END_REF]. Only one load case is considered for this application, referring to the 1 g cruise flight. The optimization method used in the NASTRAN SOL 200 is the Modified Method of Feasible Directions (MMFD, [START_REF] Vanderplaats | An efficient feasible directions algorithm for design synthesis[END_REF]).

Figure 13 plots the objective function value with respect to the re-analysis number. For each external iteration, twenty FE re-analysis are performed, and there are ten external iterations. Every FE re-analysis loop begins with an evaluation of the objective function. In order to do this, the updated loads are applied on the previous optimal configuration. These evaluations generate the "peaks" noticeable every 21 iterations, resulting from the loads dependency to the structural design variables. The value of those peaks decrease in such a way that after the eighth external iteration, no peak is observable. Thus, the loads are no longer varying significantly between two external iterations. This situation indicates that the optimization process was able to size the structure using realistic loads. Indeed, the sizing loads are in accordance to those computed with the final design. Figure 14 plots the twist laws obtained for the first, the second and the last external iteration respectively. The black line corresponds to the initial twist law generated with the updated loads applied on the former optimum point (before the twenty reanalysis). The blue twist law is obtained with the same loads but this time applied on the optimized structure (end of the twenty reanalysis). Finally, the red line is the target twist law. For each iteration, a good agreement is found between the blue (optimized) twist law and the target twist law. Naturally, the black (initial) twist law has the same behavior as the peaks from Figure 13: the discrepancies between this twist law and the target one decrease while the external iteration number increases, until a perfect matching at the end of the external optimization loop. This simple application proves the effectiveness of the strategy proposed to solve the inverse design problem. A more complex optimization will be set up soon, taking into consideration multiple load cases and, obviously, introducing the gradients into the process.

CONCLUSIONS

In this paper two approaches are presented to compute high-fidelity aero-structure gradients. The first one is based on the classical intrusive direct and adjoint approaches, and the second one is an uncoupled non-intrusive method benefiting from the linearized aerodynamic theory. Total derivatives of structural loads, and lift and drag coefficients, with respect to structural design parameters, are computed on the Common Research Model test-case, using both intrusive and non-intrusive approaches.

Gradients obtained with the tangent approach are taken as reference since it has been validated previously on an other test-case [START_REF] Achard | High-fidelity aero-structure gradient computation techniques. application to the onera m6 wing[END_REF]. The very low discrepancies observed between direct and adjoint approach demonstrate that the duality between the two methods is preserved. These results validate the implementation of the direct and adjoint method into elsA for gradients of aeroelastic function of interest with respect to structural design variables.

The non-intrusive reconstruction gives satisfactory results compared to the other approaches. The strength of this approach relies on its potential low computational cost, but cannot overcome the adjoint method for the particular case of aerodynamic coefficient derivatives. However, it represents a very promising alternative to compute gradients of structural loads for structural sizing. Usually, for this purpose, constraint aggregation technique is used in conjunction with the adjoint method. However, limiting the number of stress responses in the sizing process may oversimplify the problem formulation, leading to poor designs. In contrast, considering a limited number of structural loads can still lead to a realistic load path and then to a predictive set of stress responses. This is obviously in favor of the direct approach which in this case exhibits a comparable cost to the adjoint approach. We remind that computing the gradients of structural loads using the non-intrusive approach without the eigenvector derivatives contribution requires only one linearized computation per mode shape. Satisfactory results have been obtained on the CRM test-case by applying this approximation. In addition, this method provides aero-structure gradients with only rigid linearized aerodynamic computations, which is interesting for those who do not have access to strongly coupled aero-structure tangent or adjoint solver. Finally, uncoupled approaches are very interesting for MDO teams autonomy, which is a determining advantage in an industrial optimization process.

Ultimately, an inverse design problem is presented, with the objective to match a target twist law on the CRM wing. A strategy is proposed for an efficient structural sizing using high-fidelity aero-structure gradients. This strategy is based on the force approximation method, considered as an alternative to a constraint aggregation approach. Indeed, due to their smoother distribution, the structural loads are easier to approximate than the stresses. Several optimization processes are envisaged, according to the strategy selected to compute the loads. Preliminary results are presented on a weak coupling strategy, where the loads are updated with high-fidelity aeroelastic analysis at each external iteration. For these first results, load sensitivities are not taken into account. Nevertheless, it is demonstrated that the target twist law is reached after few external iterations, with the aeroelastic loads corresponding to the actual design. It is expected that this number of iterations will decrease with the additional contribution of the load derivatives.

The final objective of this work is to compare the intrusive adjoint method and the non-intrusive uncoupled approach through the inverse design problem described previously. The comparison will address the applicability of the optimization strategy on realistic industrial applications, in terms of efficiency and accuracy.

Figure 1 :

 1 Figure 1: Definition of the different meshes

Figure 2 :

 2 Figure 2: Illustration of the linearized method for computing δW.

Figure 3 :

 3 Figure 3: FEM colored by structural optimization group, along with the CRM aerodynamic skin.

Figure 4 :

 4 Figure 4: First five structural design parameters.

Figure 5 :

 5 Figure 5: First nine structural mode shapes splined on the wetted surface.

Figure 6 :

 6 Figure 6: Convergence of linear gradient density residual for five parameters (left) and adjoint residual for two objective functions (right).

Figure 7 :

 7 Figure 7: Gradients of horizontal force with respect to p 1 , computed with the non-intrusive approach (green and grey bars). Reference tangent approach values (red markers) are provided.

Figure 8 :

 8 Figure 8: Gradients of horizontal force with respect to p 4 , computed with the non-intrusive approach (green and grey bars). Reference tangent approach values (red markers) are provided.

Figure 9 :

 9 Figure 9: Gradients of transverse force with respect to p 1 , computed with the non-intrusive approach (green and grey bars). Reference tangent approach values (red markers) are provided.

Figure 10 :

 10 Figure 10: Gradients of transverse force with respect to p 4 , computed with the non-intrusive approach (green and grey bars). Reference tangent approach values (red markers) are provided.

Figure 11 :

 11 Figure 11: Structural design parameters used to size the structure.

Figure 12 :

 12 Figure 12: Structural optimization strategy.

Figure 13 :

 13 Figure 13: Structural optimization convergence.

Figure 14 :

 14 Figure 14: Initial (black), optimized (blue) and target (red) twist laws for three external iterations.

Table 1 :

 1 +3.993 × 10 -3 +4.028 × 10 -3 +3.186 × 10 -3 +4.624 × 10 -3 p 2 +1.200 × 10 -2 +1.209 × 10 -2 +1.112 × 10 -2 +9.274 × 10 -3 p 3 +7.682 × 10 -3 +7.728 × 10 -3 +8.467 × 10 -3 +8.290 × 10 -3 p 4 +2.019 × 10 -3 +2.018 × 10 -3 +2.580 × 10 -3 +2.449 × 10 -3 p 5 +1.202 × 10 -4 +1.193 × 10 -4 +1.545 × 10 -4 +2.748 × 10 -4 Gradients for lift coefficients with respect to wing skin thicknesses. +2.262 × 10 -4 +2.287 × 10 -4 +1.844 × 10 -4 +2.420 × 10 -4 p 2 +6.578 × 10 -4 +6.640 × 10 -4 +6.105 × 10 -4 +4.763 × 10 -4 p 3 +3.592 × 10 -4 +3.620 × 10 -4 +3.944 × 10 -4 +4.221 × 10 -4

		Intrusive		Non-Intrusive
		Direct	Adjoint	Φ , ∂Φ/∂p	Φ
	dC l /dp	p 1 Intrusive		Non-Intrusive
		Direct	Adjoint	Φ , ∂Φ/∂p	Φ
		p 1
	dC d /dp

p 4 +8.034 × 10 -5 +8.000 × 10 -5 +1.061 × 10 -4 +1.260 × 10 -4 p 5 +4.960 × 10 -6 +4.983 × 10 -6 +8.612 × 10 -6 +1.404 × 10 -5

Table 2 :

 2 Gradients for drag coefficients with respect to wing skin thicknesses.


T Achard

C Blondeau

R Ohayon

An Uncoupled Approach to Compute Aero-Structure Gradients Using High-Fidelity CFD-CSM

Aero-structural optimization is a keystone to concurrently improve aerodynamic performance and reduce the structural mass of an aircraft. Gradient-based multi-disciplinary design optimization is actually efficient if gradients computations are fast and accurate enough. This paper proposes a new approach to calculate high-fidelity aero-structure gradients, via an uncoupled non-intrusive approach benefiting from the linearized aerodynamic theory. The main advantage of the method is the independency of the computational cost with respect to the number of structural design constraints and potentially with respect to the number of structural design parameters. The efficiency of this method is demonstrated through two test-cases: the ONERA M6 wing and the Common Research Model (CRM). Gradients of aerodynamic functions of interest are computed with respect to structural design variables, and compared with finite difference results. A good agreement between these two methods proves the efficiency and accuracy of the proposed approach. Finally, this method has potential for industrial applications where uncoupled approaches are very interesting for MDO teams autonomy.

I. Introduction

Gradient based multi-disciplinary design optimization of aircraft is currently used to improve the trade-off between mass and drag reduction. However, gradients calculation is not trivial when high-fidelity simulations are considered. Moreover, the use of composite materials is increasing, resulting in more flexible wing structures [START_REF] Weisshaarr | Divergence of Forward Swept Composite Wings[END_REF] . Thus, aerodynamic shape optimization for drag reduction has to consider the impact of structural flexibility on the aerodynamic performance. From the structural point of view, aeroelastic design (e.g. weight optimization taking into account static or dynamic aeroelastic loads) has evolved towards CFD modeling for flexible loads prediction [START_REF] Blondeau | A Bi-Level High Fidelity Aero-Structural Integrated Design Methodology -A Focus on the Structural Sizing Process[END_REF] . Therefore, high-fidelity aero-structure gradients are crucial to efficiently perform structural sizing and aero-structural optimization.

Several developments have been made over the last decade to compute aero-structure gradients using the adjoint method [START_REF] Ghazlane | Aerostructural Adjoint Method for Flexible Wing Optimization[END_REF][START_REF] Kenway | Scalable parallel approach for high-fidelity steady-state aeroelastic analysis and adjoint derivative computations[END_REF][START_REF] Brezillon | Development and application of multi-disciplinary optimization capabilities based on high-fidelity methods[END_REF] . The efficiency of this approach holds when a reduced number of responses is considered compared to a large set of design parameters. For instance, when considering aerodynamic performance optimization, drag components and lift coefficient gradients for several flight points might be needed [START_REF] Reuther | Constrained Multipoint Aerodynamic Shape Optimization Using an Adjoint Formulation and Parallel Computers[END_REF] , and the set of design parameters for controlling a wing shape can reach a few hundreds. However, this is not the case in structural sizing where failure criteria in all the elements of the structural model have to be taken into account for a large set of load cases. Typically an industrial wing box sizing for preliminary design considers up to several hundreds of design parameters and up to several tens of thousands of constraints. Nonetheless, in order to still benefit from the elegant adjoint formulation, some authors have proposed a constraint aggregation approach [START_REF] Akgn | Efficient Structural Optimization for Multiple Load Cases Using Adjoint Sensitivities[END_REF][START_REF] Poon | An Adaptive Approach to Constraint Aggregation Using Adjoint Sensitivity Analysis[END_REF] as an attempt to reduce the size of the set of admissible constraints, and consequently the number of gradients to compute. However, this technique has several drawbacks. Typically, the well-known Kreisselmeier-Steinhauser function requires a parameter tuning study. An improper choice of this parameter can lead to a conservative estimate of the design space, or at the opposite to an ill-conditioned optimization problem [START_REF] Lambe | Multidisciplinary Design Optimization of an Aircraft Wing via a Matrix-Free Approach[END_REF] .

An alternative approach to calculate high-fidelity aero-structure gradients is proposed in this paper, via an uncoupled non-intrusive approach benefiting from the frequency linearized aerodynamic theory [START_REF] Mortchelewicz | Prediction of aircraft transonic aeroelasticity by the linearised Euler equations[END_REF][START_REF] Hall | A Linearized Euler Analysis of Unsteady Transonic Flows in Turbomachinery[END_REF] . By aero-structure gradients, the authors refer to gradients of any function of interest (aerodynamic or structural) with respect to structural design parameters. This strategy has been introduced by Blondeau et al. [START_REF] Blondeau | Recent achievements towards aero-structure gradient computation using high-fidelity CFD-CSM in the Onera elsA software[END_REF] , and is developed further in this paper.

Working with the elsA aeroelastic module, a steady coupled aero-structural system for a cruise condition is typically solved using a fixed point strategy. First, the fluid equations are advanced for a prescribed number of iterations. Then the wall aerodynamic forces are computed and transferred to the structural model. The corresponding structural displacements are deduced, using either a flexibility matrix or a modal approach. Next, these structural displacements are smoothed on the aerodynamic skin. Finally, these boundary displacements are propagated into the interior of the fluid mesh, and the fluid equations are solved again. This procedure is repeated until the structural displacements no longer evolve between two iterations.

The present approach linearizes the flow equations around the steady state aero-structural solution. Formally, given a field of flow conservative variables W and a set of structural parameters p, we have a direct access to the total derivative dW/dp. From here, the sensitivities of aerodynamic functions of interest (e.g. drag, lift coefficients dC d /dp and dC l /dp) are easily calculated, as well as the sensitivities of the aeroelastic loads dQ s /dp. The main advantage of the method is the independency of the computational cost with respect to the number of structural design constraints and potentially with respect to the number of structural design parameters. Another advantage is that the direct and adjoint methods are intrusive, meaning that there is a need to implement these approaches directly into the source code, as opposed to the proposed non-intrusive approach which only makes use of standard I/O. This paper is organized as follows. Section 2 proposes a theoretical description of the method. Sections 3 and 4 present the results for the gradients obtained with the uncoupled approach, applied to the ONERA M6 wing, and to the Common Research Model (CRM) [START_REF] Vassberg | Development of a Common Research Model for Applied CFD Validation Studies[END_REF] respectively.

II. Theoretical description of the uncoupled approach

The technique described in this paper allows a complete decoupling of the structure and aerodynamic disciplines. The idea stems from the frequency linearized Euler or RANS solvers used for flutter stability analysis. In this section, our objective is to demonstrate how to compute the sensitivity of the structural aeroelastic loads Q s with respect to a structural design parameter p.

II.A. Working process with the linearized solver

Let's denote p a set of structural design parameters, U the structural displacements (size 6n s , with n s the number of structural nodes), W the field of aerodynamic conservative variables and Q s the nodal structural loads. The linearized solver determines the perturbed field of conservative variables δW with respect to the steady state flow, when the fluid domain boundary is perturbed by a prescribed displacement induced by a perturbation of a structural design parameter δp. The corresponding perturbed structural loads δQ s are then deduced directly from the solution δW.

This process is illustrated in Figure 1. Flexibility is taken into account into the steady solution and the computation of δW is performed by a linearization around the equilibrium mesh position denoted Flight shape.

The elsA CFD software developed by ONERA for complex external and internal flow simulations and for multi-disciplinary applications involving aerodynamics 14 embeds an aeroelastic solver elsA/AEL (Aeroelasticity) [START_REF] Girodroux-Lavigne | Recent Navier-Stokes aeroelastic simulations using the elsA code for aircraft applications[END_REF] capable of computing in particular the requested flight shape of a coupled aero-structural system, starting with a jig shape and a flow condition. A frequency linearized module elsA/LUR (Linearized Unsteady RANS) [START_REF] Dufour | Contrasting the Harmonic Balance and Linearized Methods for Oscillating-Flap Simulations[END_REF][START_REF] Liauzun | Flutter Numerical Computations Using the Linearized Navier-Stokes Equations[END_REF] is also embedded, which computes a periodic complex solution δ * We iωt to a prescribed wall deformation δX a e iωt . The point here is to compute the solution for a null excitation frequency in order to obtain a static response which naturally corresponds to the real part of the perturbed flow solution.

Formally, the simplified equation of the linearized solver for a zero frequency is given by:

A 0 (W s )δW = B 0 (W s )δX a (1)
where A 0 and B 0 are constant real matrices depending on the field of flow conservative variables at steady state W s , δW is the unknown perturbed flow conservative variable field, δX a is a prescribed perturbation of the mesh and δW depends linearly on δX a .

The mesh perturbation is formulated as δX a = [T vol (X a )]δU int with [T vol (X a )] a linear volumic operator propagating a perturbation at the fluid-structure interface δU int into the interior of the fluid mesh X a .

Projecting the structural displacement field U on a modal basis such that U = Φq with Φ the truncated set of structural eigenvectors and q the associated generalized coordinates, δW can be written as:

δW = A 0 (W s ) -1 B 0 (W s )[T vol (X a )][T U surf ]δ(Φq) (2)
with [T U surf ] a linear operator smoothing a structural displacement field onto the aerodynamic skin. Outputs of the linearized solver are:

• δW(Φ) = Y[T U surf ]Φ, with Y = A -1 0 B 0 [T vol ]
the kernel matrix of the flow solver;

• the generalized aerodynamic forces matrix GAF.

Let's define the matrix Z, of dimension [6n s , 6n s ] as:

Z = [T Q surf ][T W surf ]Y[T U surf ] (3)
[T Q surf ] denotes a linear load transfer operator, and [T W surf ] is the operator computing the primitive variables on the aerodynamic skin from the conservative variables W in the volume mesh. The matrix Z is helpful to link the perturbed structural loads δQ s to a perturbation on the structural displacements δU:

δQ s = ZδU (4)
The GAF matrix, of dimension [n Φ , n Φ ], is defined as GAF = Φ T ZΦ.

II.B. Modal reconstruction of aero-structure gradients from an aeroelastic configuration

Let's now consider the steady equilibrium corresponding to a fluid state W, a mesh position X a and a structural displacement field U:

KU = Q s = [T Q surf ]Q a (W, X a ) (5)
K is the stiffness matrix of the finite element model, Q s the nodal structural loads, Q a the aerodynamic loads computed on the aerodynamic skin. If the full set Φ of structural eigenvectors was available, the modal projection U = Φq would hold exactly. Then substituting in Eq. ( 5) gives the corresponding generalized coordinates as:

q = γ -1 Φ T Q s ( 6
)
where γ = Φ T KΦ denotes the diagonal generalized stiffness matrix. In practice only a restricted set of first n Φ eigensolutions is computed and the modal approximation to U becomes:

U ≈ U Φ = nΦ i=1 Φ i q i (7)
However, remind that U is computed exactly from Eq. ( 5) such that the residual term in Eq. ( 7) is known from the simple difference U res = U -U Φ . Eq. ( 5) can now be reformulated as

K(U Φ + U res ) = Q s (8)
Using the equations of the linearized aerodynamic solver to compute δQ s = ZδU, the gradient of structural loads takes the simple form

dQ s dp = Z dU dp (9)
Inserting the modal decomposition of U in this expression leads to

dQ s dp = Z(Φ dq dp + ∂Φ ∂p q) + Z ∂U res ∂p (10)
The use of partial derivatives in the equation above means that Φ and U res only depend on p. This latter assumption is valid if U res is considered as a static residual mode (i.e. a structural displacement under an assumed prescribed load case).

In order to exploit this relation, the gradient of the generalized coordinates has to be determined. First left multiply Eq. ( 8) by Φ T , then differentiate with respect to p, and finally expanding using Eq. (10) to obtain

dq dp = [γ -GAF] -1 ( ∂Φ T ∂p (Q s -KU res ) + Φ T Z ∂Φ ∂p q - ∂γ ∂p q + Φ T Z ∂U res ∂p ) (11)
As we only consider stiffness sizing in this study, we introduce the usual assumption of a constant mass matrix, that is

∂M ∂p = 0 (12)
Furthermore, we choose to normalize the modes with respect to the mass matrix: Φ T MΦ = I. Then recognizing that:

∂Φ T ∂p (Q s -KU res ) = 0 (13)
Eq. ( 11) becomes:

dq dp = [γ -GAF] -1 (Φ T Z ∂Φ ∂p q - ∂γ ∂p q + Φ T Z ∂U res ∂p ) (14)
The residual vector gradient in Eq. ( 10) is approximated from the usual static sensitivity equation

∂U res ∂p = -K -1 ( ∂K ∂p U res ) (15)
However, in practice, the contribution of this term to the total gradient of loads dQ s /dp is negligible. Equation ( 14) then takes the following form:

dq dp = [γ -GAF] -1 (Φ T Z ∂Φ ∂p q - ∂γ ∂p q) (16)
And the sensitivity of the structural aeroelastic loads Q s with respect to a structural design parameter p is finally written:

dQ s dp = Z(Φ dq dp + ∂Φ ∂p q) (17)

II.C. Distance between vector subspaces of same dimension

The costly terms in the expressions ( 16) and ( 17) are those which need some linearized aerodynamic computations: ZΦ (n Φ calculations) and Z(∂Φ/∂p)q (n p calculations), with n Φ the number of modes, and n p the number of structural parameters. n Φ cannot be reduced but is quite small (10 modes are usually enough to correctly approximate the structural displacement), whereas n p may be very high (up to several hundreds).

In this subsection, we address the question of the assessment of the amount of useful information contained in the eigenvector sensitivities δΦ q = (∂Φ/∂p)q with respect to the eigenvector basis. One way to have access to this important information is to compute the Grassmann distance between the vector subspaces spanned by the modal basis and the eigenvector sensitivities. We remind here two notions:

• the Grassmann manifold G(m, n) is the whole linear subspaces of dimension m from R n ;

• the Grassmann distance between two vector subspaces, considered as two points on the Grassmann manifold, is the length of the smallest geodesic connecting these two points.

The Grassmann distance is then computed as indicated by Hamm et al. [START_REF] Hamm | Grassmann Discriminant Analysis: a Unifying View on Subspace-Based Learning[END_REF] . Let A and B be two

orthonormal matrices of size [m, m]. The principal angles 0 ≤ θ 1 ≤ θ 2 ≤ ... ≤ θ m ≤ Π
2 between two vector subspaces span(A) and span(B) are defined recursively by:

cos θ j = max aj ∈span(A) max bj ∈span(B) a T j .b j such that    a T j .a 1 = a T j .a 2 = ... = a T j .a j-1 = 0, a T j .a j = 1 b T j .b 1 = b T j .b 2 = ... = b T j .b j-1 = 0, b T j .b j = 1 (18)
The first principal angle θ 1 is therefore the smallest angle of all pairs of unit vectors in the two vector subspaces. Wong 19 linked the principal angles to the geodesic distance by:

d 2 G (A, B) = m j=1 θ 2 j ( 19
)
d G is the Grassmann distance that we aim to compute. Using singular value decomposition (SVD) helps to construct the two orthonormal bases from each vector subspaces spanned by Φ and δΦ q respectively, as well as to compute the principal angles.

The Grassmann distance aggregates an important amount of data (two matrices Φ and ∂Φ/∂p) to a scalar for each structural parameter. This scalar measures the relevancy of adding the eigenvector sensitivities in the modal reconstruction of the aero-structure gradient. For instance, let's consider two structural parameters p 1 and p 2 . Let's assume that the Grassmann distance associated to p 1 is lower than the distance associated to p 2 . This information means that if dQ s /dp 1 and dQ s /dp 2 are reconstructed without calculating Z(∂Φ/∂p 1 )q and Z(∂Φ/∂p 2 )q respectively, a better precision should be obtained for dQ s /dp 1 .

The applications presented in sections 3 and 4 will aim to demonstrate in particular the utility of the Grassmann distance criterion.

III. M6 wing test-case

In this section, we apply the uncoupled approach to the ONERA M6 wing. Our objective is to compute lift and drag coefficient sensitivities with respect to structural parameters. First an inviscid fluid (Euler) and then a viscous fluid (RANS) are used. The computational structured meshes depicted in Figure 2 contain 1.11 million cells divided into five blocks for the Euler grid and 3.8 million cells divided into 42 blocks for the RANS grid. For both models, a centered scheme is used for the spatial discretization in conjunction with the Jameson-Schmidt-Turkel artificial dissipation term, and a Backward-Euler time integration scheme. Since the M6 wing is only a pure aerodynamic test-case, a structural model has been built and designed to represent a realistic civil transport aircraft wing. The finite element model (Figure 3) has a classical wing box layout with spars, ribs and stiffeners. Member thicknesses and sections have been designed in a preprocessing optimization step. This model can be easily tuned in order to control flexibility and consequently aeroelastic effects. Finally, the structural mode shapes smoothed on the aerodynamic grid and retained to approximate the structural displacements in Eq. ( 7) are presented in Figure 4. The Mach number for the Euler case is M=0.83 and the incidence is 3.0 • . After a static aeroelastic computation, the corresponding wing tip displacement ends up at 9.13 cm (i.e. 7.77% of span). Both density residual convergence history and pressure coefficient distribution at equilibrium are presented in Figure 5. The first aeroelastic coupling is performed at the 600th iteration, and then the coupling operation repeats every 200 iterations, for a total number of 3000 iterations. Local pressure coefficient sectional distributions are presented in Figure 6, at 50% and 82% of the wingspan respectively. A dual shock, typical of the M6 wing, is clearly observed on the upper skin. The structural flexibility modifies the pressure coefficient by reducing the shock amplitude at the leading edge and moving the second shock toward the trailing edge.

For the viscous case, a Mach number of 0.84 and an incidence of 3.06 • are considered, as well as a Spalart-Allmaras one-equation turbulence model. A wing tip displacement of 14.3 cm (i.e. 12.2% of span) is obtained. The density residual convergence history is presented in Figure 7 along with the pressure coefficient distribution at equilibrium.

In this preliminary study two design parameters have been defined for the structural model. First parameter p 1 controls the thickness of the lower and upper skin, and second parameter p 2 affects the thickness of the spar webs, as illustrated in Figure 8.

III.A. Grassmann distance computation

In order to assess the amount of useful information contained in the eigenvector sensitivities δΦ q = (∂Φ/∂p)q with respect to the eigenvector basis, for the two parameters p 1 and p 2 considered, the Grassmann distance is computed. Results are presented in Table 1. Distances reported in Table 1 have been normalized by the highest possible value of the Grassmann distance:

d 2 Gmax = nΦ i=1 (Π/2) 2 .
The distance is lower for the structural parameter p 1 . As explained subsection II.C, it is expected that a better precision should be found for dQ s /dp 1 than for dQ s /dp 2 if the reconstruction of the gradients was achieved without δΦ q . p 1 p 2 Relative distance 0.487 0.503

III.B. Results for the Euler test-case

For this case, Figure 9 and Figure 10 show respectively the lift and drag coefficient sensitivities in two sections at 50% and 82% of the wing span with respect to the structural parameter p 1 , compared to finite differences. As an attempt to assess the computational cost reduction, the results are presented with and without the contribution of δΦ q . These local results look excellent, the modal reconstruction is very close to the finite difference results for the Euler test-case. The contribution of eigenvector derivatives is not significant for this parameter, since the reconstructions with and without eigenvector derivatives look very similar.

Figure 11 and Figure 12 represent respectively the lift and drag coefficient sensitivities in two sections at 50% and 82% of the wing span with respect to the structural parameter p 2 , compared to finite differences. The results are presented with and without the contribution of eigenvector derivatives. For this second parameter, the reconstruction without eigenvector derivatives is of less quality which points out their relevant contribution in this particular case. This behavior was correctly predicted by the Grassmann distance criterion (see subsection III.A).

The integrated results on the wing skin are presented in Table 2. Reference finite difference values are also provided.

The sign and the order of magnitude are always well predicted. Adding eigenvector sensitivities is beneficial for the second parameter, but this contribution does not seem as relevant as for p 1 . Nevertheless, the uncoupled approach for this test-case is validated. Further investigations are needed in order to understand the contribution of δΦ q = (∂Φ/∂p)q.

III.C. RANS test-case

Figure 13 and Figure 14 show respectively the lift and drag coefficient sensitivities in two sections at 50% and 82% of the wing span with respect to the structural parameter p 1 , compared to finite differences. The results are presented with and without the contribution of eigenvector derivatives. All linearized computations have been performed using a constant turbulent viscosity, i.e. using the "frozen viscosity approximation". However turbulence sensitivity is taken into account within finite difference calculations. This approximation can lead to minor observable deviations, considering that finite difference results are used as reference for comparison. Similarly to the Euler case, the viscous local results are fairly good: the modal reconstruction is very close to the finite difference reference. Figure 15 and Figure 16 represent respectively the lift and drag coefficient sensitivities in two sections at 50% and 82% of the wing span with respect to the structural parameter p 2 , compared to finite differences. The results are presented with and without the contribution of eigenvector derivatives.

Again, similar conclusions to the Euler test-case can be made. The integrated results on the wing are presented in Table 3. Reference finite difference values are reported.

The same conclusions for this test-case can be drawn, since the sign and the order of magnitude of the gradient values are always well predicted. From the observation that the added contribution of the eigenvector sensitivities was not so significant for the first parameter in the Euler test-case (see Table 2), it was decided to improve the smoothing technique by carefully choosing the structural nodes used for building the spline kernel matrix. As the structural wing box model does not extend up to the leading and trailing edges of the aerodynamic model, the smoothing technique needs to extrapolate information at these locations. After modification, it can be observed in Table 3 that the contribution of the eigenvector sensitivities is beneficial for both parameters. This points out that the accuracy of the smoothing step is crucial for the proposed approach. Additional work still needs to be performed in order to improve the results when eigenvector sensitivities are taken into account. Indeed the generalized vector δΦ q exhibits a rather irregular shape with many local peaks, and is thus very challenging for the smoothing technique. Nevertheless, the uncoupled approach for the viscous test-case is also validated.

IV. Application to the Common Research Model (CRM)

The CRM configuration is a modern transonic aircraft designed by NASA in 2009 in order to develop new experimental databases to validate specific applications of CFD [START_REF] Vassberg | Development of a Common Research Model for Applied CFD Validation Studies[END_REF] . This demand was specifically expressed by the AIAA CFD Drag Prediction Workshop 4 (DPW4).These databases provide a better understanding of numerical simulations, typically for drag prediction. NASA developed a transonic aircraft configuration with a supercritical wing profile, a fuselage, a horizontal tail plane (HTP), a nacelle and a pylon. The CRM is designed for a Mach number of 0.85, for a lift coefficient of 0.5 at 37000 ft. In these conditions, the Reynolds number is about 40 million. The fuselage is representative of a wide-body long-range commercial transport aircraft. Beyond the drag and moment prediction, some authors chose the CRM configuration as a starting point for aerodynamic optimization purpose [START_REF] Lyu | Aerodynamic Shape Optimization Investigations of the Common Research Model Wing Benchmark[END_REF][START_REF] Carrier | Gradient-Based Aerodynamic Optimization with the elsA Software[END_REF] . Due to its realistic configuration, it is also of interest for specific applications like structural sizing or aerostructural optimization. However, realistic aircraft structural models are rarely open-source, since aircraft manufacturers protect their know-how. Moreover, the initial objective of the development of the CRM was to obtain a test-case for CFD validation, so that only the deformed 1 g cruise flight shape is available. Some authors proposed to design a finite element model of the CRM wing. For instance, starting from the geometrical definition of the CRM, DLR designed a parametric finite element model for their FERMAT configuration. This model was used to achieve a static aeroelastic trimming analysis, a structural sizing and a flutter prediction analysis [START_REF] Klimmek | Parametric Set-Up of a Structural Model for FERMAT Configuration for Aeroelastic and Loads Analysis[END_REF] . A process deriving the jig shape of the wing was embedded during the structural sizing. Kenway et al. in 2014 23 designed a structural jig shape of the CRM, taking advantage of the similarities between the CRM and the Boeing 777 geometries, in order to achieve an aerostructural optimization of the wing. This section is organized as follow. First, the fluid and the structural models are presented. Then, the process to find the CRM jig shape is described, and finally the results obtained with the uncoupled approach are discussed.

IV.A. Fluid model

In this paper, we have chosen to work with the wing/body/HTP (WBH) configuration. Different refinement levels are available, but in order to limit the computation time, we select the "coarse" model of 5.4 million cells, split into 26 blocs. For all the subsequent results, a RANS fluid model is used with a Jameson centered scheme for the spatial discretization, along with a Backward-Euler algorithm for time integration. The Spalart-Allmaras one-equation turbulence model is selected. In the attempt to find a lift coefficient of 0.5 for a Mach number of 0.85 at 37000 ft, an angle of attack of 2.127 • is found. Figure 17 shows the pressure coefficient distribution on the WBH configuration.

IV.B. Finite element model (FEM)

Several finite element models of the CRM wing are available on the NASA CRM website. 24 All these models correspond to the flight shape of the CRM. Our objective is to prove the efficiency and accuracy of the uncoupled approach to compute gradient of flexible loads with respect to structural parameters. In this context, we are only interested in a structural model which exhibits a realistic static behavior. We tested the different versions of the finite element model of the wing, solving equation 5 to find the displacements corresponding to the aerodynamic loads Q a computed on the aerodynamic skin. The model "V14", exhibited an expected static behavior, with a realistic vertical displacement of 2.15 m observed at the wing tip in nominal cruise conditions. It is a simple finite element model of the wing box and centre box with spars, skins, ribs, implicit stiffeners and a constant wing skin thickness of 8.89 mm. This model can be considered as a good starting point for a structural sizing process. To meet the independence assumption of the mass matrix M with respect to the design variables (see Eq. 12), shell density has been replaced by non structural mass (NSM). The structural displacement and force nodes used to build the flexibility matrix for the static aero-structure coupling are depicted in Figure 19. The provided CRM configuration corresponds to the 1 g cruise shape. In order to find the corresponding jig shape, an inverse procedure is set up and presented in subsection IV.C. To correctly predict the structural displacement due to gravity, a realistic distribution of concentrated masses representing the fuel in the wing is added, along with concentrated masses for the engine, pylon and landing gears. The total lift generated by the wing is equal to the overall mass of the aircraft. The cumulative inertial loads (transverse force, bending and twisting moments) computed along the X aerodynamic axis are presented in Figure 20, as well as the integrated loads from aerodynamic and inertial. The localized effects of the landing gears and of the powerplant masses, respectively located at 2.5 m and 9 m on the reference integration axis, are clearly observable on the transverse force and twisting moment distributions.

Figure 21 presents the finite element model colored by the structural optimization groups, along with the aerodynamic skin of the CRM. Ten structural parameters are chosen for the demonstration of the gradient computation with the uncoupled approach, split into two groups of five parameters driving the thickness of the lower skin and the upper skin respectively. Table 4 describes the first nine eigenmodes computed on our modified FEM of the CRM wingbox. The corresponding mode shapes smoothed on the aerodynamic skin are presented in Figure 22. Figure 23 represents the evolution of the norm < U res , U res > K with respect to the number of retained modes, U res being computed such that U res = U -U Φ (see Eq. 7). It can be observed that n Φ = 8 modes are enough to correctly approximate the structural displacement U. This truncation criterion provides an interesting information as first approximation to compute the gradients with the uncoupled approach.

Mode number Frequency

Another useful practicality, based on the Generalized Cross Validation criterion (GCV), can help to determine the relevancy of each mode in the reconstruction of the structural displacement U. The principle is to compute the reconstruction error norm < U res , U res > K on a modal basis leaving the ith eigenvector out, for i ∈ [1, n Φ ]. Results are plotted Figure 24. Obviously, the first mode has a major contribution to the structural displacement. As illustrated by the close-up view, the contribution of the higher modes can be discarded. This information is very interesting for the gradient reconstruction. However, it has to be checked that using only the first mode yields a sufficient accuracy compared with finite difference results. Therefore, the GCV criterion analysis is essential for potential computational cost reduction.

IV.C. Inverse procedure to determine the CRM jig shape

As mentioned above, only the 1 g cruise shape is given for the aerodynamic and the structural model. An iterative inverse procedure is used to find the corresponding jig shape for both the structure and aerodynamic models. Let's denote X 0 a the initial fluid mesh, corresponding to the flight shape. We assume that a steady state solution W 0 around the flight shape X 0 a is available. The corresponding aerodynamic loads Q a (W 0 , X 0 a ) are computed once and will be noted Q a . At iteration i, the fluid mesh is noted X i a . The structural mesh, aligned with the flight shape, is noted X s and its stiffness matrix is noted K. For the sake of describing as carefully as possible the procedure, we reintroduce more precisely the three operators needed to manage the fluid-structure interactions:

• [T Q surf (X i a , X s )
] the linear operator transferring the aerodynamic loads Q a to the structural mesh X s . This operator is constructed between the current fluid mesh X i a and the structural mesh X s .

• [T U surf (X s , X i a )] the surface operator smoothing the structural displacements U onto the aerodynamic skin of the fluid mesh X i a . This operator is constructed between the structural mesh X s and the fluid mesh X i a .

• [T vol (X i a )] the linear volumic operator propagating the smoothed displacement U a i at the fluid-structure interface into the interior of the fluid mesh X i a .

The following procedure enables to determine both the aerodynamic and structural jig shapes:

1. Computation of the structural loads

Q 0 s = [T Q surf (X 0 a , X s )]Q a applied on X s .
2. Resolution of the static equation KU 0 = Q 0 s and construction of U 0 s = U 0 +U g with U g the structural displacement due to gravity.

Smoothing operation to determine

U 0 a = [T U surf (X s , X 0 a )]U 0 s 4.
Propagation of U 0 a into the interior of the volumic mesh to find the new mesh X

1/2 a = [T vol (X 0 a )]U 0 a 5. Obtention of the new computational fluid mesh X 1 a = X 0 a -(X 1/2 a -X 0 a )
6. Set iteration counter i=1. We define ∆U i = U i -U i-1 and ∆U 0 is initialized such that ∆U 0 > 100ǫ. While ||∆U i || 2 > ǫ:

(a) Computation of the structural loads

Q i s = [T Q surf (X i a , X s + U i ini )]Q a with U i ini = -U i-1 s (b) Resolution of the static equation KU i = Q i s and construction of U i s = U i + U g (c) Smoothing operation to determine U i a = [T U surf (X s + U i ini , X i a )]U i s (d) Propagation of U i
a into the interior of the volumic mesh to find the new mesh

X i+1/2 a = [T vol (X i a )]U i a (e) Obtention of the new computational fluid mesh X i+1 a = X 0 a -(X i+1/2 a -X 0 a ) (f) i=i+1
Following this procedure, the structural and aerodynamic jig shapes are found in three iterations with ǫ = 2.10 -5 , and we control that the structural jig shape matches perfectly with the aerodynamic jig shape. Figure 25 summarizes the inverse procedure described above, and Figure 26 shows the superimposition of the aerodynamic flight shape with the converged aerodynamic jig shape.

IV.D. Grassmann distance computation

In order to assess the amount of useful information contained in the eigenvector sensitivities δΦ q = (∂Φ/∂p)q with respect to the eigenvector basis, the Grassmann distance is computed for the ten parameters considered. Results are presented in Figure 27. The values of the distance are very similar between the parameters from the lower skin and those from the upper skin. Parameters p 1 and p 6 seem to be less sensitive to the information brought by δΦ q compared with the other parameters. On the contrary, parameters p 4 and p 9 seem to require the contribution of δΦ q to obtain a good precision for the gradients.

The information contained in the nominal modal vector subspace is common to all the design parameters. Besides, the Grassmann distance provides an additional information about the modal derivative contribution specific to each parameter. It is thus possible to sort these distances to assess the relevancy of this latter contribution which then provides a relative weight associated to any parameter.

IV.E. Results for the CRM test-case

This subsection presents the first results obtained with the CRM test-case. The objective is to compute gradients of aerodynamic functions of interest, like lift and drag coefficients, with respect to the ten structural design parameters. Remind that these gradients are computed with equations 16 and 17. The only change with respect to the previous test-cases concerns the effect of gravity, which is now taken into account in the structural displacements. Equations 16 and 17 are still valid if the generalized coordinates take the form:

q = γ -1 Φ T (Q s + Q gravity ) (20)
with Q s containing the aerodynamic and the elastic forces transferred to the structural nodes, and Q gravity representing the structural load due to gravity.

Table 5 and Table 6 present the derivatives for the integrated coefficients for parameters driving respectively the thicknesses of the lower and upper skin. Finite difference values are provided as a reference.

The results obtained for the lower skin and for the upper skin are fairly similar, which could be expected since they have the same thickness distribution. For most of the structural design parameters, these inte- grated results are in good agreement with finite difference results: the sign and the order of magnitude are well predicted. However, for parameters p 5 and p 10 situated at the tip of the wing, the gradients differ on the order of magnitude or on the sign. These discrepancies may come from the chosen step for the finite differences: the study realized to tune the best step for these two parameters exhibited strong nonlinearities. Therefore, further investigations are necessary to find a better step to compute the aero-structure gradients.

The smoothing technique could also be improved, especially for the smoothing of the eigenvector derivatives on the aerodynamic skin, as underlined in subsection III.C. In the following, the discussion of these results will skip parameters p 5 and p 10 .

For the lift coefficient derivatives, adding eigenvector derivatives is beneficial for all parameters excepted p 4 . In order to have a better insight of the effect of the eigenvector derivatives for this parameter, local results have been plotted in the section at 82% of the wing span (see Figure 28). This graph clearly confirms that the contribution of eigenvector derivatives is needed to correctly reconstruct gradients for parameter p 4 . This trend was effectively predicted by the Grassmann distance criterion (see subsection IV.D). A careful inspection of local results is therefore necessary to validate the integrated ones.

For the drag coefficient, the contribution of the eigenvector derivatives is beneficial for all parameters excepted for p 1 and p 6 . Drag coefficient sensitivities are plotted in Figure 29 for a section at 50% of the wing span for design variable p 6 . Again, the contribution of eigenvector derivatives enhances the reconstruction. However, this improvement is not as significant as for parameter p 4 , which is in accordance with the Grassmann distance criterion (see Figure 27). It is now proposed to verify the relevancy of the two criteria proposed in subsection IV.B. Table 7 presents the convergence of the lift coefficient derivatives with respect to parameters p 1 to p 4 , when the first mode, the first 8 modes and the first 20 modes are taken into account for the reconstruction (without the contribution of eigenvector derivatives). Parameters p 6 to p 9 are not printed since their behavior is really similar to parameters p 1 to p Two conclusions can be drawn from Table 7. No noticeable change is observed with the additional contribution of the last set of modes. This confirms that the modal truncation criterion is suitable for the choice of the modal basis size (see Figure 23). The second conclusion is even more important: using only the first mode to reconstruct the gradients provides a very fair estimate. This was also pointed out by the GCV criterion (see Figure 24). In this case, it indicates that a single linearized computation should give a good approximation of the gradients with respect to all the structural design parameters.

V. Conclusion

In this paper, a new strategy is proposed to compute high-fidelity aero-structure gradients via an uncoupled non-intrusive approach benefiting from the linearized aerodynamic theory. The theoretical description is provided, and results are discussed for two test-cases: the ONERA M6 wing and the CRM configuration.

The M6 wing validation case considers inviscid and viscous flow simulations. In both cases, local and global values of the aero-structure gradients obtained with the present uncoupled approach are validated with finite difference results. More specifically, the proposed approach considers two contributions for the gradient reconstruction. The first one only uses information related to the nominal modal basis, and the second one takes advantage of the additional information contained in the eigenvector derivatives with respect to structural design parameters. Even when only the modal basis part is considered (i.e. the formulation is independent on the number of design variables), a good accuracy can be obtained for the integrated coefficient gradients.

Results obtained for the CRM demonstrate the potential of the method when a more realistic test-case is selected. In a first attempt, both contributions from nominal modal basis and eigenvector derivatives have been taken into account. For almost all structural parameters, the uncoupled approach results are in good agreement with finite differences. However, for some parameters, it happens that inaccuracies can occur in the integrated results which requires an additional local inspection of the associated pressure gradient distribution. Specifically, gradients associated to parameters driving the skin thickness near the wing tip are not accurate enough. We believe that the following explanations can be formulated. First, the tuning process for the finite difference step size was very nonlinear. Secondly, this method relies heavily on the quality of the smoothing process applied to the modeshape derivatives. Further investigations will be performed in the near future to relieve these limitations. Nevertheless, the accuracy of the method has been demonstrated.

In a second attempt, the question of possibly dropping the costly contribution of eigenvector derivatives has been addressed. We remind that this contribution requires one aerodynamic solution for each parameter. To assess this contribution, the Grassmann distance criterion has been proposed. This criterion measures the amount of relevant information contained in the eigenvector derivatives, compared with the information provided by the modal basis. Although theoretically promising, it is difficult to interpret in practice. However, comparisons of the results with and without eigenvector derivative contributions are always consistent with the criterion prediction. But, in its present form, it does not seem to be sufficiently selective to be used without post-validation. However, apart from this discussion, considering the only contribution of the mode shapes still yields a predictive gradient reconstruction.

The size of the modal basis is selected using a truncation error threshold for the static displacement field approximation. In order to validate this criterion, a convergence study has been performed with respect to the size of the projection basis. It is effectively observed that the contribution of modes above this cut-off is negligible. These results have also been consolidated by the Generalized Cross Validation (GCV) method, which is expected to exhibit the most relevant modes in the modal basis.

To conclude, the strength of this method remains in its potentially very low computational cost, while using high-fidelity CFD/CSM simulations. Indeed, due to the independency of the computational cost with respect to the number of design constraints and possibly the number of structural design parameters, it is a very promising alternative to the constraint aggregation techniques used in conjunction with the adjoint method. Furthermore, uncoupled approaches are very interesting for MDO teams autonomy, which is a determining advantage in an industrial optimization process.

Ongoing work is dedicated to the implementation of the intrusive gradient capability counterpart in the elsA software. Once operational, this new module will allow further comparisons with the present approach and give us several alternatives for high-fidelity aero-structure gradient computation.

Figure 1 .

 1 Figure 1. Illustration of the linearized method for computing δW.

Figure 2 .

 2 Figure 2. M6 wing Euler and RANS mesh.

Figure 3 .

 3 Figure 3. M6 wing structural mesh and selected displacement/force nodes for flexibility matrix construction.

Figure 4 .

 4 Figure 4. First six structural mode shapes splined on the wetted surface.

Figure 5 .

 5 Figure 5. Aeroelastic equilibrium shape colored by pressure coefficient (left). History of density residual norm for aeroelastic computation (right). Euler test-case.

Figure 6 .

 6 Figure 6. Comparison between rigid and flexible Cp distribution at section y=50% and y=82% span. Euler test-case.

Figure 7 .

 7 Figure 7. Aeroelastic equilibrium shape colored by pressure coefficient (left). History of density turbulent viscosity residual norm for aeroelastic computation (right). RANS test-case.

Figure 8 .

 8 Figure 8. Structural design parameters.

Figure 9 .

 9 Figure 9. Lift coefficient sensitivities with respect to p1 (in red), compared to finite differences (in blue). Results presented with the contribution of eigenvector derivatives (left) and without (right). Euler test-case.

Figure 10 .

 10 Figure 10. Drag coefficient sensitivities with respect to p1 (in red), compared to finite differences (in blue). Results presented with the contribution of eigenvector derivatives (left) and without (right). Euler test-case.

Figure 11 .

 11 Figure 11. Lift coefficient sensitivities with respect to p2 (in red), compared to finite differences (in blue). Results presented with the contribution of eigenvector derivatives (left) and without (right). Euler test-case.

Figure 12 .

 12 Figure 12. Drag coefficient sensitivities with respect to p2 (in red), compared to finite differences (in blue). Results presented with the contribution of eigenvector derivatives (left) and without (right). Euler test-case.

Figure 13 .

 13 Figure 13. Lift coefficient sensitivities with respect to p1 (in red), compared to finite differences (in blue). Results presented with the contribution of eigenvector derivatives (left) and without (right). Viscous test-case.

Figure 14 .

 14 Figure 14. Drag coefficient sensitivities with respect to p1 (in red), compared to finite differences (in blue). Results presented with the contribution of eigenvector derivatives (left) and without (right). Viscous test-case.

Figure 15 .

 15 Figure 15. Lift coefficient sensitivities with respect to p2 (in red), compared to finite differences (in blue). Results presented with the contribution of eigenvector derivatives (left) and without (right). Viscous test-case.

Figure 16 .

 16 Figure 16. Drag coefficient sensitivities with respect to p2 (in red), compared to finite differences (in blue). Results presented with the contribution of eigenvector derivatives (left) and without (right). Viscous test-case.

Figure 17 .

 17 Figure 17. Wing-body-HTP configuration colored by pressure coefficient at M=0.85

Figure 18 .

 18 Figure 18. Slices of pressure coefficient on the wing at M=0.85

Figure 19 .

 19 Figure 19. Selected displacement nodes (green) and force nodes (red) for flexibility matrix construction.

Figure 20 .

 20 Figure 20. Inertial loads (left) and running loads (right) integrated along the X aerodynamic axis.

Figure 21 .

 21 Figure 21. FEM colored by structural optimization group, along with the CRM aerodynamic skin

Figure 22 .

 22 Figure 22. First nine structural mode shapes splined on the wetted surface.

Figure 23 .

 23 Figure 23. < Ures, Ures > K norm with respect to the number of retained modes. The red dotted line represents the truncation criterion of the modal basis.

Figure 24 .

 24 Figure 24. Leave-one-out based method applied on the truncated modal basis.

Figure 25 .

 25 Figure 25. Schematic representation of the inverse procedure

Figure 27 .

 27 Figure 27. Relative Grassmann distances for the ten structural parameters considered

Figure 28 .

 28 Figure 28. Lift coefficient sensitivities with respect to p4 (in red), compared to finite differences (in blue). Results presented with the contribution of eigenvector derivatives (left) and without (right). CRM test-case.

Figure 29 .

 29 Figure 29. Drag coefficient sensitivities with respect to p6 (in red), compared to finite differences (in blue). Results presented with the contribution of eigenvector derivatives (left) and without (right). CRM test-case.

Table 1 .

 1 Relative Grassmann distances for structural parameters p1 and p2

Table 2 .

 2 Gradients for lift and drag coefficients with respect to wing skin and spar thickness. Euler test-case.

			FD	Φ,∂Φ/∂p	Φ
	dC l /dp	p 1	+2.34E-02	+2.12E-02	+2.59E-02
		p 2	-5.23E-03	-5.24E-03	-7.61E-03
	dC d /dp	p 1	-3.81E-04	-6.38E-04	-4.39E-04
		p 2	-5.08E-04	-4.98E-04	-5.23E-04

Table 3 .

 3 Gradients for lift and drag coefficients with respect to wing skin and spar thickness. Viscous test-case.

			FD	Φ,∂Φ/∂p	Φ
	dC l /dp	p 1	+3.63E-02	+3.37E-02	+3.98E-02
		p 2	-3.53E-03	-4.05E-03	-7.62E-03
	dC d /dp	p 1	-4.47E-04	-5.11E-04	-3.20E-04
		p 2	-5.37E-04	-5.52E-04	-6.21E-04

Table 4 .

 4 First nine eigen modes

Table 5 .

 5 Gradients for lift and drag coefficients with respect to wing thickness -lower skin. CRM test-case.

			FD	Φ,∂Φ/∂p	Φ
		p 1	+6.67e-01	+4.37e-01	+9.08e-01
		p 2	+3.02e+00	+2.26e+00	+1.92e+00
	dC l /dp	p 3	+2.23e+00	+2.00e+00	+1.73e+00
		p 4	+4.63e-01	+5.41e-01	+4.93e-01
		p 5	-4.64e-02	+7.45e-03	+4.61e-02
		p 1	-2.08e-02	-1.15e-02	-2.61e-02
		p 2	-2.90e-02	-3.66e-02	-5.75e-02
	dC d /dp	p 3	-6.25e-02	-6.82e-02	-5.25e-02
		p 4	-3.32e-02	-3.26e-02	-1.42e-02
		p 5	-3.61e-03	-1.35e-03	-1.46e-03
			FD	Φ,∂Φ/∂p	Φ
		p 6	+5.85e-01	+4.28e-01	+8.91e-01
		p 7	+2.98e+00	+2.29e+00	+1.88e+00
	dC l /dp	p 8	+2.32e+00	+2.15e+00	+1.71e+00
		p 9	+5.56e-01	+5.70e-01	+4.93e-01
		p 10	+8.31e-03	+4.88e-03	+4.64e-02
		p 6	-2.15e-02	-1.06e-02	-2.57e-02
		p 7	-3.77e-02	-3.75e-02	-5.66e-02
	dC d /dp	p 8	-6.51e-02	-7.29e-02	-5.18e-02
		p 9	-3.37e-02	-3.38e-02	-1.42e-02
		p 10	-2.31e-02	-1.29e-03	-1.45e-03

Table 6 .

 6 Gradients for lift and drag coefficients with respect to wing thickness -upper skin. CRM test-case.

Table 7 .

 7 Gradients for lift coefficients with respect to wing thickness -lower skin. CRM test-case.

		p 1	p 2	p 3	p 4
	FD	+6.67e-01	+3.02e+00	+2.23e+00	+4.63e-01
	First mode	+7.40e-01	+1.69e+00	+1.54e+00	+4.23e-01
	First 8 modes	+9.40e-01	+1.92e+00	+1.73e+00	+4.95e-01
	First 20 modes	+9.08e-01	+1.92e+00	+1.73e+00	+4.93e-01

4

.