Christophe Blondeau
email: christophe.blondeau@onera.fr

Timothée Achard

Philippe Girodroux-Lavigne

Roger Ohayon
email: roger.ohayon@cnam.fr

Recent achievements towards aero-structure gradient computation using high-fidelity CFD-CSM in the Onera ELSA software

Keywords: optimization, aeroelastic gradient, adjoint technique

INTRODUCTION

Design of aeroelastic and aerodynamic performance has benefited not only from recent achievements in terms of numerical tools for computation of flexible gradients through tangent or adjoint approach but also from the re-organization of the optimization process itself to deal with high-fidelity aerodynamic and structural models [START_REF] Blondeau | A Bi-Level High Fidelity Aero-Structural Integrated Design Methodology -A Focus on the Structural Sizing Process[END_REF][START_REF] Salah El Din | Transonic Wing-body Civil Transport Aircraft Aero-Structural Design Optimization using a Bi-Level High Fidelity Approach -A Focus on the Aerodynamic Process[END_REF]. It is now recognized that rigid aerodynamic shape optimization for drag reduction has to consider the impact of structural flexibility on the performance especially in the context of multi-point optimization. On the other hand, aeroelastic design (i.e. structural weight optimization by tailoring static or dynamic aeroelastic loads) has evolved towards CFD modelling 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]. This paper presents the latest achievements about the development of a high-fidelity based aero-structure gradient computation capability in the framework of the elsA software [START_REF] Cambier | Status of the elsA Software for Flow Simulation and Multi-Disciplinary Applications[END_REF][START_REF] Girodroux-Lavigne | Recent Navier-Stokes aeroelastic simulations using the elsA code for aircraft applications[END_REF].

By aero-structure we mean gradient of any function of interest with respect to aerodynamic shape or structural design parameters. Up to now, structural modelling for flexible gradient computation was limited to a beam-like equivalent model for the wing box [START_REF] Marcelet | Sensitivity Analysis of a Strongly Coupled System Using the Discrete Direct and Adjoint Approach[END_REF]. More recently this approach was extended to get sensitivities of an aerodynamic function of interest (drag component, lift coefficient, L/D ratio …) with respect to structural design parameters such as wall thicknesses of skins and spars for an ideal rectangular wing section, but still limited to an equivalent beam model kinematics [START_REF] Dumont | Overview of Recent Development of Aeroelastic Adjoint Method for Civil Aircraft Wing Optimization[END_REF][START_REF] Ghazlane | Aerostructural Adjoint Method for Flexible Wing Optimization[END_REF]. Although useful for preliminary design studies and design space exploration, this structural modelling 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 the present work was to relief all these limitations by re-developing an aerostructure gradient capability from scratch by systematic differentiation of discretized aeroelastic equations, and all associated operators, in the elsA aeroelasticity module (elsA/Ael). This way, the new capability inherits the structural paradigm embedded into the elsA/Ael module (structural flexibility matrix), the whole catalogue of fluid-structure transfer methods, as well as the underlying parallel architecture. High-fidelity aero-structure gradient computations and optimization studies using the adjoint approach have been reported in the literature for structured and unstructured grids [START_REF] Mavriplis | A Discrete Adjoint-Based Approach for Optimization Problems on Three-Dimensional Unstructured Meshes[END_REF][START_REF] Nielsen | Recent Improvements in Aerodynamic Design Optimization on Unstructured Meshes[END_REF][START_REF] Fazzolari | Efficient aerodynamic shape optimization in MDO context[END_REF][START_REF] Kenway | Scalable Parallel Approach for High-Fidelity Steady-State Aeroelastic Analysis and Adjoint Derivative Computations[END_REF][START_REF] Gauger | Efficient Deterministic Approaches for Aerodynamic Shape Optimization[END_REF]. However, none of them seems to have a core memory efficient on-the-fly implementation for deriving algebra manipulations such as matrix-vector product involving the geometric or state residual Jacobian matrix. In the following, equations and solving procedure for the aeroelastic gradient are discussed as well as the numerical control parameters and their impact on the solution convergence. Preliminary results for static aeroelastic gradient computations with respect to aerodynamic shape parameters, in linear and adjoint mode, for the Onera M6 wing will be presented. Euler and RANS fluid modelling will be considered. On going extension work of the aeroelastic gradient tool in order to consider structural variables will also be presented. Two approaches can be considered. The first one, called intrusive approach, needs a substantial modification of the elsA source code. The second one, called non-intrusive, considers a projection of the static displacement field on the structural mode shapes in order to benefit from the classical harmonic forced response mode of linearized Euler or RANS simulation tools. Preliminary results for this second approach will be presented for the Onera M6 wing with respect to stiffness parameters of the primary wing box structure.

SENSITIVITY OF A STRONGLY COUPLED AEROELASTIC SYSTEM

In order to introduce notations and remind typical solving techniques, a brief review of sensitivity analysis methods for aerodynamic and structural systems of discretized equations is given first. The interested reader should refer to [START_REF] Peter | Numerical sensitivity analysis for aerodynamic optimization: A survey of approaches[END_REF] for a complete review of numerical sensitivity analysis.

As the subsequent derivations will always be related to discretized equations, it will not be reminded elsewhere in the following discussion. The generalization to a strongly coupled fluid-structure system will then be presented. The subsequent derivations are valid for any set of design parameters (i.e. impacting local aerodynamic shape only, structural stiffness only or both aerodynamic shape and stiffness as global design variables).

Gradient evaluation for a mono-disciplinary system

Any discretized set of steady-state nonlinear equilibrium equations can be written in the following residual form [START_REF] Blondeau | A Bi-Level High Fidelity Aero-Structural Integrated Design Methodology -A Focus on the Structural Sizing Process[END_REF] Where and represent the vector of design variables and the corresponding dependent state vector. represents either a second-order accurate cell-centered finite volume spatial discretization of the fluid residual or a finite element spatial discretization of the structural residual. Note that Eq. ( 1) has an implicit dependency on the computational aerodynamic grid or structural grid . To simplify notations, this dependency is only reproduced when necessary in subsequent expressions.

is supposed to be once continuously differentiable with respect to the state variables and computational mesh in the vicinity of and .

We are interested in the sensitivity derivatives of any set of function with respect to a design variable. Applying the chain rule differentiation to yields:

(2)

The total sensitivity derivatives are computed from the differentiation of the discrete residual Eq. ( 1)

( ) [ ] (3)
where the Jacobian matrix is assumed invertible. Note that the gradient vanishes because the residual equation has to be satisfied for any value of . Usually terms like , or are obtained analytically (full differentiation by hand) or numerically (automatic differentiation, finite differences or complex step).

Adjoint equations for a scalar function of interest are derived from Eq. ( 2) by substituting for from Eq. (3) to obtain [START_REF] Girodroux-Lavigne | Recent Navier-Stokes aeroelastic simulations using the elsA code for aircraft applications[END_REF] where is a vector of adjoint variables that satisfies

[ ] (5)
This expression is independent of the design variables , which makes it particularly attractive for aerodynamic optimization problems with a large number of design parameters and relatively few responses of interest.

Considering an aerodynamic problem Eq. ( 1) is a system of nonlinear equations with unknowns. As an example, let's consider the functional dependency of a drag coefficient with respect to a set of shape design parameters: [START_REF] Dumont | Overview of Recent Development of Aeroelastic Adjoint Method for Civil Aircraft Wing Optimization[END_REF] where represents the fluid conservative variables. Functional is the surface mesh parameterization that links design parameters to wetted surface grid positions.

is the volumic mesh deformation operator that propagates the surface displacements to the interior fluid domain:

. The fluid states are obtained by solving the discretized fluid equations under prescribed boundary conditions. Eq. ( 6) exhibits an additional direct dependency to the grid and is thus a generalization of Eq. ( 1). Direct differentiation of Eq. ( 6) wrt. gives [START_REF] Ghazlane | Aerostructural Adjoint Method for Flexible Wing Optimization[END_REF] Then differentiating the discrete residual wrt. leads to the following expression for the total derivatives of the fluid states:

[ ] (8)
Terms and are computed analytically around the equilibrium solution that satisfies .

The adjoint counterpart of Eq. ( 7) reads ) with the vector of adjoint variables that satisfies

[ ] (9
[ ] (10)
A standard Backward Euler iterative scheme is used to solve for the static equilibrium. The discrete residual is linearized around the flow solution at current iteration and a pseudo time step is introduced along with a suitable approximation of the Jacobian matrix [START_REF] Peter | Large Stencil Viscous Flux Linearization for the Simulation of 3D Compressible Turbulent Flows with Backward-Euler Schemes[END_REF] to form the implicit term in the following preconditioned fixed point iterative scheme:

( ) (11)
Applying this algorithm to the linearized form of the residual (around the steady state solution) in Eq. ( 8) gives the following iterative scheme for the total derivative of the flow field :

* + ( ) (12)
In a similar manner the iterative formula for the adjoint vector is derived,

* + [ ] (13)
Clearly expressions in Eq. ( 12) and ( 13) are very similar and it is straightforward to show that the following property holds [START_REF] Peter | Large Stencil Viscous Flux Linearization for the Simulation of 3D Compressible Turbulent Flows with Backward-Euler Schemes[END_REF] An exact dual algorithm has been proposed in [START_REF] Nielsen | An implicit, exact dual adjoint solution method for turbulent flows on unstructured grids[END_REF] to enforce this property at every iteration.

Let's now consider the classical structural mechanics example of finding the total gradient of a constraint response in a finite element. The system is discretized by the finite element method to give the residual discrete system of equations [START_REF] Nielsen | An implicit, exact dual adjoint solution method for turbulent flows on unstructured grids[END_REF] where and represent the nodal displacements and the applied loads. The stiffness matrix is positive semi-definite and self-adjoint.

Assuming a general linear elastic behavior, the relationship between stresses and strains will be linear of the form , where denotes the elasticity matrix and and the initial strain and stress fields. The deformation field is obtained from the nodal displacements through the usual relation , being a linear differential operator and the matrix of shape functions. Finally the constraint function reads [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF] Direct differentiation of Eq. ( 16) wrt. a scalar parameter yields

( ) (17)
In the above expressions the geometric dependency of the structural mesh to the design parameter has been emphasized. Once again, for sake of clarity and without any loss of generality it will be omitted in the following when necessary.

From Eq. ( 15) direct differentiation gives

( ) (18)
leading to the adjoint form of the sensitivity equation:

( ) ( ) (19)
where, using the symmetric property of , the adjoint matrix satisfies [START_REF] Gallard | Aerodynamic aircraft design for mission performance by multipoint optimization[END_REF] Unlike the aerodynamic systems of equations, Eq. ( 15) and ( 20) are usually solved efficiently using direct parallel sparse multifrontal algorithms, see [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF].

Extension to an aeroelastic problem

For an aeroelastic system the state vector is defined as and the objective function is of the form . Design parameters are linked to structural properties, control the aerodynamic shape and global variables affect simultaneously stiffness and shape. The design parameter vector is denoted as . The vector of state equations reads

( ) (21)
Direct differentiation of Eq. ( 21) wrt.

gives

[ ] ( ) ( ) (22)
and solving for allows the computation of the objective function gradient [START_REF] Mortchéléwicz | Prediction of aircraft transonic aeroelasticity by the linearised Euler equations[END_REF] The two blocks in Eq. ( 22) are coupled through aerodynamic loads applied on the structural skin, which induce an elastic deformation which in turn alters fluid grid positions.

A suitable transferring technique is applied to to compute structural loads of the form [START_REF] Wang | Improved approximate methods for computing eigenvector derivatives in structural dynamics[END_REF] The reference aerodynamic mesh aligned with the unloaded configuration (i.e. the jigmesh) has been introduced. This distinction is useless for the structural mesh which always matches the jig-shape.

denotes the linear load transfer operator. Knowing updated structural grid positions from displacements , it is possible to transfer them back as on the aeroelastic interface through the linear splining step [START_REF] Hamm | Grassmann Discriminant Analysis: A Unifying View on Subspace-Based Learning[END_REF] Once again, the linear displacement transfer operator only depends on reference meshes. Finally the fluid interior domain deformation is attained in a second step by propagating the surface deformation with (26) where the volumic operator merely depends on the reference aerodynamic mesh.

Thanks to Eq. ( 24) through (26) it is possible to expand Eq. ( 22) and ( 23) in terms of unknown vectors and . Firstly Eq. ( 26) is linearized wrt. to give

(27)
From Eq. ( 25) the total derivative of the boundary displacement field is written as

(28)
Substituting for and from Eq. ( 27) and (28) into Eq. ( 23) leads to the following expanded form of the objective function:

*( ) ( )+ (29)
Terms and only depend on the equilibrium steady state and are calculated analytically in a pre-processing step. Grid sensitivity terms and represent the geometric parameterization which is also part of an external dedicated procedure. Other terms like and and their associated partial derivatives are related to the transfer techniques used during the aeroelastic steady state computation.

By performing similar algebra manipulations, the discrete block system in Eq. ( 22) can be cast into the compact form

{ ([ ] [ ] [ ] ) [ ] [ ] [ ] [ ] [ ] (30)
Constant matrices [A], …, [H] are defined analytically with the following formulas

[ ] [ ] (31) [ ] [ ] [ ][ ] [ ] [ ][ ] [ ] [ ] [ ][ ] [ ]
where the utility matrix has been introduced.

The system in Eq. ( 30) is presented in a form suitable for an iterative block scheme resolution similar to the lagged-block strategy formerly proposed in [START_REF] Martins | A coupled adjoint method for high-fidelity aero-structural optimization[END_REF][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 Eq. ( 29) becomes

([ ] [ ] [ ] ) (32)
Yet again, for each residual state equation in Eq. ( 21) an adjoint vector is introduced which multiplies the linearization of that equation, and used to formulate the augmented objective function derivative such that for any vector and

([ ] [ ] [ ] ) (33)
Expanding and and factorizing out the difficult terms

⁄ et ⁄ yields ( [ ]) ( [ ] [ ] [ ] ) ([ ] [ ] ) ([ ] [ ] ) ([ ] [ ] [ ] ) (34)
where adjoint vectors and are chosen such that the following system is satisfied:

[ [ ] [ ] [ ] [ ] [ ] ] ( ) ( [ ] [ ] [ ] ) (35)
Finally, once adjoint unknowns are determined, the total gradient is obtained with

([ ] [ ] ) ([ ] [ ] ) ([ ] [ ] [ ] ) (36)

GRADIENT COMPUTATION FOR AERODYNAMIC SHAPE OPTIMIZATION

The different test cases for inviscid and viscous flow are based on the ONERA M6 wing which has been extensively used for three-dimensional transonic CFD validation purposes, see [START_REF] Schmitt | Pressure Distributions on the ONERA M6 Wing at Transonic Mach Numbers[END_REF].

Aeroelastic analysis and design models

The computational structured meshes depicted in Figure 1 The finite element model has a classical wing box layout with spars, ribs and stiffeners. Member thicknesses and sections have been designed in a pre-processing optimization step. This model can be easily tuned in order to control flexibility and consequently aeroelastic effects. The Mach number for this case is M=0.734 and the incidence is 2.08°. The wing tip displacement is 6.27cm (i.e. 5.14% of span) and the associated lift coefficient is C z =0.196 compared to a rigid coefficient of C z =0.221. An upwind Roe scheme with a MUSCL interpolation associated to a Van Albada limiting function is applied. The Spalart-Allmaras turbulence model has been selected.

Figure 3 below shows the pressure distribution at equilibrium and the residual norm convergence for density and turbulent variables. A first step of 600 fluid iterations is performed followed by a coupling every 200 iterations. The aeroelastic simulation exhibits an excellent convergence up to machine precision. The idea is to ensure as much as possible the cancellation of the residual that is the basic assumption for accurate discrete gradient computations. Three design parameters have been defined for specific assessment of gradient accuracy. First parameter x 01 controls the section camber at 50% span and affects linearly the region from root to 90% of wing span. Second parameter x 02 is the quarter chord twist angle at the wing tip section, and the last parameter x 03 is the angle of attack.

Inviscid flow -Euler computations

Two gradient computations have been carried out in tangent and adjoint mode for the Euler test case. The objective is to compute total gradients for lift and pressure drag coefficients wrt. the three design parameters previously defined. The convergence of the discrete residual norm of the fluid block equations is reported in Figure 4. In direct mode, three right hand sides are considered, i.e. one for each design parameter whereas only two right hand sides are considered in adjoint mode, i.e. one for each objective function. In both direct and adjoint mode an excellent convergence is reached. Table 1 summarizes the gradient values and allows a comparison with rigid results and finite differences. It is worth noting that flexibility effects dramatically change rigid gradient values and even reverse the sign of the drag sensitivity wrt. the camber parameter. One can observe a residual duality gap between direct and adjoint results. This is because some geometrical terms are still approximately linearized in the adjoint formulation. These approximations are currently being removed.

Viscous flow -RANS computations

The same transonic flow conditions M=0.734 and =2.08° have been retained for the viscous test case.

Figure 5 presents the convergence histories for the density residual in tangent mode. All computations have been performed using a constant eddy and turbulent viscosity, i.e. using the standard "frozen  t approximation". Gradient values are summarized in Table 2 and can be directly compared with Euler values in Table 1. Obviously viscous effects significantly alter inviscid gradient values even changing the sign of the lift coefficient derivative with respect to the camber design variable. Finite difference values compare quite well with tangent results. Currently, the viscous adjoint mode is still in a validation phase and we were not able to provide reference results for the M6 wing in this paper.

GRADIENT COMPUTATION FOR STRUCTURAL SIZING

The adjoint approach developed in section three is efficient whenever a reduced number of responses is considered compared to a large set of design parameters. Typically for an aerodynamic design dedicated to performance improvement, drag components and lift coefficient gradients for several flight points might be needed [START_REF] Gallard | Aerodynamic aircraft design for mission performance by multipoint optimization[END_REF][START_REF] Reuther | Constrained Multipoint Aerodynamic Shape Optimization Using an Adjoint Formulation and Parallel Computers[END_REF]. The set of design parameters for controlling a wing shape will not exceed a few hundreds. Unlikely, structural sizing for certification has to consider failure criteria at many locations under 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 as an attempt to reduce the size of the set of admissible constraints, and thus the number of gradients to compute. However, this technique has several drawbacks and may lead to sub-optimal designs. Instead we propose a gradient reconstruction based on a modal approach. This technique allows a complete decoupling of the structure and aerodynamic blocks. The idea stems from the (frequency) linearized Euler or RANS solvers used in the construction of state-space models for flutter stability analysis in aeroelasticity. These solvers compute complex harmonic responses (pressure coefficients, generalized harmonic forces) to prescribed boundary displacement fields (i.e. structural mode shapes in our case).

Modal reconstruction of loads sensitivity around the steady state solution

Using previous notations, the steady equilibrium corresponding to a fluid state ̅ , a mesh position ̅ and a structural displacement field is written as

̅ ̅ (37)
If the full set of structural eigenvectors were available, the modal projection holds and then substituting in Eq. (37) gives the corresponding generalized coordinates as (38) where denotes the diagonal generalized stiffness matrix.

In practice only a restricted set of first eigensolutions is computed and the modal approximation to becomes

∑ with ∑ (39)
However, remind that is computed exactly from Eq. ( 37) such that the residual term in Eq. ( 39) is known from the simple difference . Now Eq. ( 37) can be reformulated as (40)

Pre-multiplying this equation by , then substituting for and subsequently differentiating wrt. a single design parameter p gives (41)

In the above expression it is assumed that the sensitivity of the residual load vector is zero. This makes sense if the residual vector is thought of being an assumed deformation under a prescribed static load thus leading to .

From linearized Euler or Navier-Stokes equations theory, see for instance [START_REF] Jameson | In Time dependent calculation using multigrid, with application to unsteady flows past airfoil and wings[END_REF][START_REF] Mortchéléwicz | Prediction of aircraft transonic aeroelasticity by the linearised Euler equations[END_REF], it is possible to introduce a suitable first order approximation to the aerodynamic loads ̅ ̅ in terms of a small perturbation of the steady state equilibrium.

Recalling that only design parameters linked to physical properties are considered here, the system in Eq. (30) simplifies to

{ [ ] [ ] [ ] (42)
Solving for from the first block of equations and then substituting in the structural block gives the following system for the displacement derivatives

{ [ ] [ ] [ ] [ ] (43)
Now direct differentiation of Eq. ( 37) yields (44)

and from direct comparison with the first block in Eq. ( 43) the gradient of structural loads takes the simple form (45)

Inserting the modal decomposition of in this expression leads to

( ) (46)
In order to exploit this relation, the gradient of the generalized coordinates has to be determined first. This is done by inserting Eq. (46) into Eq. (41) to obtain

[ ] ( ) (47)
where the generalized aerodynamic forces matrix is defined as .

This equation can be simplified further if we introduce the additional assumption that design parameters only affect the stiffness matrix, i.e. , being the mass matrix. Then assuming that eigenmodes are mass-normalized, the following property is easily demonstrated (48)

Recognizing that

, then the first term in brackets in the right hand side of Eq. (47) vanishes. Finally the expression for reads

[ ] ( ) (49)
The residual vector gradient is approximated from the usual static sensitivity equation

( ) (50)
To summarize the proposed strategy for the gradient reconstruction of aerodynamic forces, each vector requires at least linearized rigid aerodynamic responses: responses to mode shapes plus one response to the generalized vector . It is possible to reduce this cost to linearized simulations by taking advantage of the work by Wang in [START_REF] Wang | Improved approximate methods for computing eigenvector derivatives in structural dynamics[END_REF], who proposed an accurate approximation of an eigenvector derivative that is assumed to be spanned by a truncated set of normal modes augmented by a residual static mode such that ∑

where and

and the pseudo-load matrix . Note that in this case under previous assumption . If a sufficient number of retained modes is considered, then the residual term becomes negligible and the linearized aerodynamic response is readily obtained from a linear combination of already available responses to mode shapes.

Aeroelastic analysis and design models

The structural and Euler fluid meshes have already been presented in Figure 2, but the Mach number for this case is M=0.83 and the incidence is 3.0°. Only an inviscid Euler flow is considered here. The wing tip displacement is 9.13 cm (i.e. 7.48% of span) and the associated lift coefficient is C z =0.232. The steady aeroelastic computation uses an upwind Roe scheme with a Van Albada limiter. The time scheme is a Backward-Euler, with a classical pseudo time step. The density residual convergence history is presented in the figure below along with the pressure coefficient distribution at equilibrium.

Gradient computation and results

Figure 10 and Figure 11 represent the results obtained respectively for the reconstruction of the pressure coefficient gradient w.r.t first and second design parameter, compared with the gradients obtained by finite difference method. For both parameters the reconstruction is very close to the finite difference results, which validates the linearized approach proposed in this study. In order to reduce the number of computations, the accuracy of the gradient reconstruction is assessed without the contribution of eigenvector derivatives in Eq. ( 46) and (49). Corresponding results are presented in Figure 12 for parameter . Table 3 below summarizes the reconstructed gradients for lift coefficient. Obviously, this reconstruction is of less quality which points out the relevant contribution of eigenvector derivatives. In order to understand the relevancy of adding eigenvector derivatives to the equations, the distance between subspaces spanned by the eigenvectors and their derivatives has been calculated. As these subspaces are of the same dimension, the Grassmann distance using the principal angles can be computed [START_REF] Hamm | Grassmann Discriminant Analysis: A Unifying View on Subspace-Based Learning[END_REF] to give: 0.479 0.493 Table 4 : Relative Grassmann distance between subspaces spanned by eigenvectors and eigenvector derivatives wrt. skin thickness and spar web thickness .

The relative Grassmann distance (i.e. normalized by maximum value) for is less than the one for . This means that the information provided by the eigenvector derivatives is more useful in the latter case than in the former.

CONCLUSION

The elsA block-structured highly scalable software has been enriched with a new capability for aeroelastic gradient computation for aerodynamic shape optimization. This new module is no more restricted to simplified structural models and benefits from the highly parallel architecture of elsA. Theoretical details of the proposed direct and adjoint formulation are provided as well as numerical preliminary results for the M6 wing test case for inviscid and viscous fluid models. Thanks to advanced relaxation and acceleration techniques, excellent convergence rates have been observed. In addition, a non-intrusive modal approach is proposed to obtain an accurate reconstruction of gradients with respect to structural stiffness parameters. This technique only requires harmonic responses to prescribed boundary displacements on the wetted area and benefits from the existing linearized Euler or RANS solvers used in the construction of state-space models for flutter stability analysis in aeroelasticity. Results for an inviscid transonic test case are reported for the M6 wing.

  contain 1.11 million cells divided into five blocks for the Euler grid and 3.8 million cells divided into 42 blocks for the Navierstokes grid.

Figure 1 :

 1 Figure 1: M6 wing Euler and RANS mesh.

Figure 2 :

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

Figure 3 :

 3 Figure 3: Cp distribution (left), convergence of density and turbulence residual norm for RANS aeroelastic simulation (right).

Figure 4 :

 4 Figure 4: Convergence of linear density residual norm for three parameters (left) and adjoint density residual norm for two objective functions (right).

Figure 5 :

 5 Figure 5: Convergence of linear gradient density residual for three parameters.

Figure 6 :

 6 Figure 6 : Aeroelastic equilibrium shape colored by pressure coefficient (left). History of density residual norm during the fixed-point aeroelastic computation (right)

Figure 10 :

 10 Figure 10 : Reconstruction of the pressure coefficient gradient wrt. compared to finite differences (dashes).

Figure 11 :

 11 Figure 11 : Reconstruction of the pressure coefficient gradient wrt. compared to finite differences (dashes).

FD

Figure 12 :

 12 Figure 12 : Reconstruction of the pressure gradient wrt. second design parameter compared to finite difference computation, without eigenvector derivatives contribution

Table 1 :

 1

		Rigid Adjoint	FD	Flex. Tangent	Flex. Adjoint
	x 01	+1.56e+00	+8.13e-01	+8.141e-01	+8.057e-01
	x 02	+5.95e-03	+8.51e-04	+8.539e-04	+9.688e-04
	x 03	+1.30e-01	+7.13e-02	+7.142e-02	+6.943e-02
		Rigid Adjoint	FD	Flex. Tangent	Flex. Adjoint
	x 01	+1.92e-01	-1.30e-02	-1.415e-02	-1.900e-02
	x 02	+2.85e-04	+5.62e-05	+5.648e-05	+5.367e-05
	x 03	+1.14e-02	+3.37e-03	+3.386e-03	+3.000e-03

Euler gradients for lift and pressure drag coefficients wrt. camber (x 01 ), tip twist (x 02 ) and angle of attack (x 03 ). Comparison between rigid, finite difference and flexible gradients in tangent and adjoint mode.

Table 2 :

 2

Navier-Stokes gradients for lift and pressure drag coefficients wrt. camber (x 01 ), tip twist (x 02 ) and angle of attack (x 03 ). Comparison between rigid, finite difference and flexible gradients in tangent mode.

Table 3 :

 3 Gradients for lift coefficient wrt. wing skin and spar thickness. Values with and without eigenvector derivatives contribution. Reference finite difference values are provided.

The first aeroelastic coupling is done at iteration n. 600, then the coupling operation repeats every 200 iterations, for a total number of 3000 iterations. Two design parameters have been defined for the structural model. First parameter controls the thickness of the linked lower and upper skin, and second parameter affects the thickness of the spar webs. Finally, the smoothed mode shapes used to approximate the structural displacements are presented in Figure 9.

Skins Spars

Additional complex 3D geometry aircraft applications such as wing-body and complete aircraft are currently being considered to assess the full potentiality of the new module and will be reported elsewhere. Further developments will extend the embedded aeroelastic gradient module to support structural design parameters as well.

ACKNOWLEDGEMENTS

The authors wish to thank Airbus Group for partially funded this work. We also address many thanks to Antoine Dumont from the Applied Aerodynamics Dept. for providing the aerodynamic meshes, shape parameterization inputs and finite difference results. We also greatly appreciated fruitful discussions with Jacques Peter from the Computational Fluid Dynamics Dept. for his knowledge of the elsA software and advices for debugging strategies.

COPYRIGHT STATEMENT

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


Walid Larbi
email: walid.larbi@cnam.fr

Jean-Franc ¸ois Deü

Roger Ohayon

Vibro-Acoustic Analysis of Laminated Double-Wall: Finite Element Formulation and Reduced-Order Model

Keywords: double-wall, viscoelastic, vibroacoustic, finite element, modal reduction

This paper presents a finite element model for sound transmission analysis through a double sandwich panels with viscoelastic core inserted in an infinite baffle. The proposed model is derived from a multi-field variational principle involving structural displacement of the panels and acoustic pressure inside the fluid cavity. To solve the vibro-acoustic problem, the plate displacements are expanded as a modal summation of the plate's real eigenfunctions in vacuo. Similarly, the cavity pressure is expanded as a summation over the modes of the cavity with rigid boundaries. Then, an appropriate reduced-order model with mode acceleration method by adding quasi-static corrections is introduced. The structure is excited by a plane wave. The radiated sound power is calculated by means of a discrete solution of the Rayleigh Integral. Fluid loading is neglected. Various results are presented in order to validate and illustrate the efficiency of the proposed reduced finite element formulation.

Introduction

Double-wall structures are widely used in noise control due to their superiority over single-leaf structures in providing better acoustic insulation. Typical examples include double glazed windows, fuselage of airplanes and vehicles, etc. Various theoretical, experimental and numerical approaches have been investigated to predict the sound transmission through double walls [START_REF] Fahy | Sound and structural vibration, 1st edn[END_REF][START_REF] Panneton | Modélisation numérique tridimensionnelle par éléments finis des milieux poroélastiques: application au problème couplélasto-poro-acoustique[END_REF][START_REF] Quirt | Sound transmission through windows I. Single and double glazing[END_REF].

By introducing a thin viscoelastic interlayer within the panels, a better acoustic insulation is obtained. In fact, sandwich structures with viscoelastic layer are commonly used in many systems for vibration damping and noise control. In such structures, the main energy loss mechanism is due to the transverse shear of the viscoelastic core [START_REF] Akrout | Vibro-acoustic behaviour of laminated double glazing enclosing a viscothermal fluid cavity[END_REF][START_REF] Basten | Noise reduction by viscothermal acousto-elastic interaction in double wall panels[END_REF].

1
In the first part of this paper, a non-symmetric finite element formulation of double-wall sandwich panels with viscoelastic core is derived from a variational principle involving structural displacement and acoustic pressure in the fluid cavity. Since the elasticity modulus of the viscoelastic core is complex and frequency dependent, this formulation is complex and nonlinear in frequency. Therefore, the direct solution of this problem can be considered only for small size models. This has severe limitations in attaining adequate accuracy and wider frequency ranges of interest. An original reduced order-model is then proposed to solve the problem at a lower cost. The proposed methodology, based on a normal mode expansion, requires the computation of the uncoupled structural and acoustic modes. The uncoupled structural modes are the real and undamped modes of the sandwich panels without fluid pressure loading at fluid-structure interface, whereas the uncoupled acoustic modes are the cavity modes with rigid wall boundary conditions at the fluid-structure interface. It is shown that the projection of the full-order coupled finite element model on the uncoupled bases, leads to a reduced order model in which the main parameters are the classical fluid-structure and residual stiffness complex coupling factors. Moreover, the effects of the higher modes of each subsystem is taken into account through an appropriate so-called static correction. Despite its reduced size, this model is proved to be very efficient for simulations of steady-state and frequency analyses of the coupled structural-acoustic system with viscoelastic damping and the computational effort is significantly reduced.

In the last part, numerical examples are presented in order to validate and analyse results computed from the proposed formulation.

Finite Element Formulation of the Coupled Problem

Local Equations

Consider a double-wall structure coupled to an acoustic enclosure shown in Fig. 1. Each wall occupies a domain Ω Si , i ∈{ 1, 2} such that Ω S =( Ω S1 ,Ω S2 ) is a partition of the whole structure domain. A prescribed surface force density F d is applied to the external boundary Γ t of Ω S and a prescribed displacement u d is applied on a part Γ u of Ω S . The acoustic enclosure is filled with a compressible and The harmonic local equations of this structural-acoustic coupled problem can be written in terms of structure displacement u andfluidpressurefieldp [START_REF] Morand | Fluid-structure interaction[END_REF] divσ(u)+ρ S ω 2 u = 0 in Ω S (1)

Σ Ω s1 Ω s2 Ω F F d Γ t u d u d
σ(u)n S = F d on Γ t (2) σ(u)n S = pn on Σ (3) u = u d on Γ u (4) Δp + ω 2 c 2 F p = 0 in Ω F (5) ∇p • n = ρ F ω 2 u • n on Σ ( 6
)
where ω is the angular frequency, n S and n are the external unit normal to Ω S and Ω F ; ρ S and ρ F are the structure and fluid mass densities; c F is the speed of sound in the fluid; and σ is the structure stress tensor.

Constitutive Relation for Viscoelastic Core

In order to provide better acoustic insulation, damped sandwich panels with a thin layer of viscoelastic core are used in this study (Fig. 1). When subjected to mechanical vibrations, the viscoelastic layer absorbs part of the vibratory energy in the form of heat. Another part of this energy is dissipated in the constrained core due to the shear motion.

The constitutive relation for a viscoelastic material subjected to a sinusoidal strain is written in the following form:

σ = C * (ω)ε (7)
where ε denotes the strain tensor and C * (ω) is the complex elasticity tensor. It is generally complex and frequency dependent ( * denotes complex quantities). It can be written as:

C * (ω)=C ′ (ω)+iC ′′ (ω) (8)
where i = √ -1. Furthermore, for simplicity, a linear, homogeneous, viscoelastic constitutive equation will be used in this work. In the isotropic case, the viscoelastic material is defined by a complex and frequency dependence shear modulus in the form:

G * (ω)=G ′ (ω)+iG ′′ (ω) (9)
where G ′ (ω) is know as shear storage modulus, as it is related to storing energy and G ′′ (ω) is the shear loss modulus, which represents the energy dissipation effects.

With these assumptions, the stress tensor of the sandwich structure is complex and frequency dependent.

Finite Element Equation

After variational formulation using the test-function method [START_REF] Morand | Fluid-structure interaction[END_REF][START_REF] Larbi | Finite element formulation of smart piezoelectric composite plates coupled with acoustic fluid[END_REF] and discretization by the finite element method, we obtain the following matrix system of the fluid/elastic structure with viscoelastic damping coupled problem:

K * u (ω) -C up 0K p -ω 2 M u 0 C T up M p U P = F 0 (10)
where U and P are the vectors of nodal values of u and p respectively; M u and K * u (ω) are the mass and stiffness matrices of the structure; M p and K p are the mass and stiffness matrices of the fluid; C up is the fluid-structure coupled matrix; F is the applied mechanical force vector. Note that since the elasticity modulus of the viscoelastic core of the sandwich panels is complex and frequency dependent, the stiffness matrix K * u (ω) is also complex and frequency dependent.

Reduced Order Model

In this section, we introduce a reduced-order formulation based on a normal mode expansion with an appropriate static correction.

Eigenmodes of the Structure in Vacuo

In a first phase, the first N s eigenmodes of the structure in vacuo are obtained from

K * u (ω) -ω 2 M u U = 0 (11)
Due to the frequency dependent of the stiffness matrix, this eigenvalue problem is complex and nonlinear. It is assumed that vibrations of the damped structure can be represented in terms of the real modes of the associated undamped system if appropriate damping terms are inserted into the uncoupled modal equations of motion. Thus, the complex stiffness matrix is decomposed in the sum of two matrices:

K * u (ω)=K u0 + δK * u (ω) (12)
where K u0 = Re[K * u (0)] is the real and frequency-independent stiffness matrix calculated with a constant Young module's of the viscoelastic core and δK * u (ω) is the residual stiffness matrix.

The ith real eigenmode is obtained from the following equation

K u0 -ω 2 si M u Φ si = 0 for i ∈{1, ••• ,N s } (13)
where (ω si , Φ si ) are the natural frequency and eigenvector for the ith structural mode.

Eigenmodes of the Internal Acoustic Cavity with Rigid Walls

In this second phase, the first N f eigenmodes of the acoustic cavity with rigid boundary conditions are obtained from the following equation

K p -ω 2 fi M p Φ fi = 0 for i ∈{1, ••• ,N f } (14)
where (ω fi , Φ fi ) are the natural frequency and eigenvector for the ith acoustic mode.

It is important to note that the physical acoustic modes in a rigid fixed cavity are such that ΩF p dv =0which excludes the ω fi =0 , p = constant solution which is not physical but had to be introduced in the formulation of the coupled problem. Thus, the variational formulation of this kind of system, in order to be regularized for zero frequency situation, i.e. valid for a static problem, has to be modified by adding the following constraint ρ F c 2 F Σ u.n ds + ΩF p dv =0(see [START_REF] Morand | Fluid-structure interaction[END_REF] for details). When doing this, on one hand the static pressure is defined precisely by

p s = - ρ F c 2 F |Ω F | Σ u.n ds ( 15
)
and on the other hand, the reduced order formulation will be carried only by projection on the physical acoustic modes.

Modal Expansion of the General Problem

By introducing the matrices

Φ s =[ Φ s1 ••• Φ sNs ] of size (M s × N s )a n dΦ f = Φ f 1 ••• Φ fN f of size (M f × N f
) corresponding to the uncoupled bases (M s and M f are the total number of degrees of freedom in the finite elements model associated to the structure and the acoustic domains respectively), the displacement and pressure are sought as

U = Φ s q s (t) and P = Φ f q f (t) (16)
where the vectors q s =[ q s1 ••• q sNs ] T and q f = q f 1 ••• q fN f T are the modal amplitudes of the structure displacement and the fluid pressure respectively.

Substituting these relations into Eq. ( 10) and pre-multiplying the first row by Φ T s and the second one by Φ T f , we obtain

• N s mechanical equations -ω 2 q si + Ns k=1 γ * ik (ω)q si + ω 2 si q si - N f j=1 β ij q fj = F i (17)
• N f acoustic equations

-ω 2 q fi + ω 2 fi q fi -ω 2 Ns j=1 β ij q sj =0 (18)
where F i (t)=Φ T si F is the mechanical excitation of the ith mode; β ij = Φ T si C up Φ fj is the fluid structure coupling coefficient and γ * ik (ω)=Φ T si δK * u (ω)Φ sk the reduced residual stiffness complex coefficient.

At each frequency step, the reduced system (Eqs. ( 17) and ( 18)) is solved by updating γ * ik (ω). After determining the complex amplitude vectors q si and q fi , the displacement and pressure fields are reconstructed using the modal expansion (Eqs. ( 16)).

Static Corrections

The process of mode truncation introduces some errors in the response that can be controlled or minimized by a modal truncation augmentation method. In this method, the effects of the truncated modes are considered by their static effect only. First the applied loading vector F is composed as:

F = L i=1 α i (t)F 0i (19)
where F 0i is the invariant spatial portion and α j (t) is the time varying portion. For each invariant spatial load, the static modal eigenvector Ψ si is given by:

Ψ si = K -1 u0 F 0i (20)
The truncated basis containing the real and undamped structure modes is enriched by the static modal eigenvectors such that

Φs =[Φ s1 ••• Φ sNs , Ψ si , ••• Ψ sl ]=[Φ s Ψ s ] (21)
The truncated fluid basis is enriched with the static pressure P s computed from Eq. ( 15):

Φf =[Φ f P s ] (22)
Thus, the displacement and pressure are sought as

U = Φ s q s (t)+Ψ s q 0 s (t) and P = Φ f q f (t)+P s q 0 f (t) (23)
where the vectors q 0 s and q 0 f are the quasi-static modal amplitudes of the structure displacement and the fluid pressure respectively. Similar coupled differential equations than Eqs. ( 17) and ( 18) can of course be obtained after modal projection on the enriched bases Φs and Φf .

Numerical Examples

Sound Transmission through an Elastic Double-Panel System

In this section, the validation of the proposed coupled finite elements formulation for sound radiation is presented. The problem under consideration is shown in Fig. 2. A normal incidence plane wave excites a double-plate system filled with air (density ρ F = 1.21 kg/m 3 and speed of sound c F = 340 m/s). The plane wave has a pressure amplitude of 1 N/m 2 and is applied to plate 1 as the only external force to the system. The plates are identical and simply supported with thicknesses of 1 mm. The density of the plates is 2814 kg/m 3 , the Youngs modulus is 71 GPa, the Loss factor is 0.01 and Poisson ratio 0.33. The surrounding fluid is the air. This example was originally proposed by Panneton in [START_REF] Panneton | Modélisation numérique tridimensionnelle par éléments finis des milieux poroélastiques: application au problème couplélasto-poro-acoustique[END_REF]. When the excitation is applied to the first plate, the second one vibrates and radiates sound caused by the coupling of air and plate 1. The normal incidence sound transmission loss is then computed using the Rayleigh's integral method [START_REF] Fahy | Sound and structural vibration, 1st edn[END_REF] which needs the finite element solution of surface velocities of plate 2. For this purpose, the resolution of the coupled system is done with a modal reduction approach using the first 10 in vacuo structural modes and the first 10 acoustic modes of the fluid in rigid cavity. The truncated bases are enriched by the static modal eigenvectors.

Fig. 3 shows the normal incidence transmission loss through a simply supported plate (dashed line). Due to the modal behavior of the plate, dips in the transmission loss curve are observed at its resonance frequencies (modes (1, 1), (3, 1) and (1, 3)). When a second plate is used to form an airtight cavity (continuous line), an increase in the transmission loss is achieved except in the region of the so-called plate-cavity-plate resonance (mode (1, 1)*). At this frequency, the two plates move

Fig. 3 Comparison of the normal incidence sound transmission (nSTL) through an air-filled double panel and a simple panel out of phase with each other and the effect of the cavity on the plates is mostly one of added stiffness. This frequency is similar to the mass-air-mass resonance of unbounded double panels.

In addition, the variation of the nSTL of an air-filled panels and a simple panel is in very good agreement with the published date from [START_REF] Panneton | Modélisation numérique tridimensionnelle par éléments finis des milieux poroélastiques: application au problème couplélasto-poro-acoustique[END_REF].

Sound Transmission through a Double Laminated Glazing Window

The proposed reduced order finite-elements formulation is applied now to calculate the transmission loss factor of a double laminated glazing window. The system consists of two identical clamped laminated panels of glass separated by an air cavity of 12 mm thickness. Each laminated glass is composed of two glass plates bonded together by a Polyvinyl Butyral (PVB) interlayer. The thickness of outer and inner glass ply is h 1 = h 3 =3mm and those of the PVB interlayer is h 2 =1 .14 mm. The glass ply is modeled as linear elastic material (density 2500 kg/m 3 , Youngs modulus 72 GPa, and Poisson ratio 0.22). The material properties of the PVB are both thermal and frequency dependent. From dynamic and thermal tests, Havrillak and Negami have found an empirical law describing this dependence. The resulting complex frequency dependent shear modulus of the PVB is given at 20 • Cas [START_REF] Havriliak | A complex plane analysis of?-dispersions in some polymer systems[END_REF]:

G * (ω)=G ∞ +(G 0 -G ∞ ) 1+(iωτ 0 ) 1-α0 -β0 (24)
where G ∞ =0 .235 GPa, G 0 =0 .479 Mpa, α 0 =0 .46, β 0 =0 .1946, τ 0 = 0.3979. The Poisson ratio of the PVB is 0.4 and density is 999 kg/m 3 . Concerning the excitation and the finite element discretization, we used the same ones as in the previous example. A comparison between a simple glass and a laminated glass with PVB interlayer with an equivalent surface mass is shown in Fig. 4. Calculation was limited to 2000 Hz maximum. This comparison shows that laminated glass has a much lower acoustic radiation compared to conventional glass at resonance frequencies du to the effect of the viscoelastic layer. The reduction of sound radiation power is around 10 dB in lower frequencies and around 20 dB in higher frequencies. In fact, at low frequencies, the viscoelastic material is soft and the damping is small. At higher frequencies, the stiffness decreases rapidly and the damping is highest. Moreover, flexural vibrations causes shear strain in the viscoelastic core which dissipates energy and reduces vibration and noise radiation. Note that the thickness of the viscoelastic layer has a significant influence in terms of attenuation. Fig. 5 nSTL through an air-filled double panel: comparison between the modal reduction approach and the direct nodal method Fig. 5 shows a comparison between the nSTL of the coupled problem, obtained with the proposed accelerated modal reduction approach with a truncation on the first twenty structural modes (N s =20) and first twenty acoustic modes (N f =20) and the direct nodal method (Eq. ( 10)) where the displacement and pressure vectors are calculated for each frequency step. The structural modes are calculated from Eq. (13) using the constant shear storage modulus G ∞ . As can be seen, a very good agreement between the two methods is proved. In this respect, it should be noted that the resulting reduction of the model size and the computational effort using the reduced order method are very significant compared to those of the direct approach.

Conclusions

In this paper, a finite element formulation for sound transmission through double wall sandwich panels with viscoelastic core is presented. A reduced-order model, based on a normal mode expansion, is then developed. The proposed methodology requires the computation of the eigenmodes of the undamped structure, and the rigid acoustic cavity. Quasi-static corrections are introduced in order to accelerate the convergence. Despite its reduced size, this model is proved to be very efficient for simulations of steady-state analyses of structural-acoustic coupled systems with viscoelastic interlayers when appropriate damping terms are inserted into the modal equations of motion.

Fig. 1

 1 Fig. 1 Double sandwich wall structure coupled to an acoustic enclosure

Fig. 2

 2 Fig. 2 Double-plate system filled with air: geometric data

Fig. 4

 4 Fig.[START_REF] Akrout | Vibro-acoustic behaviour of laminated double glazing enclosing a viscothermal fluid cavity[END_REF] Comparison of radiated sound power from a simple glass pane and a laminated glass with the same mass


domain: spi chim

text:
Frédérick De Meyer

Eric Cloarec

Xavier Rozanska

Karen Gonzalez-Tovar

Alain Valtz

Chris- Tophe Coquelet

Inhibition of aqueous chemical reactions by strong liquid structure makers: experimental demonstration with amides and acid-base reactions

Keywords: acid-base reactions, amides, hydrogen bond, acid gas, aqueous amine, liquid structure maker, reaction inhibitor

The chemical absorption of CO 2 and H 2 S in aqueous tertiary amines are well-known acid-base reactions. Kinetic and vapor-liquid equilibrium experiments show that the addition of an amide such as HMPA, that is known to be a strong liquid structure maker, significantly inhibits the acid-base reactions. The impact is more pronounced for CO 2 than for H 2 S absorption. Despite the presence of water in the solvent, the absorption becomes almost physical. Due to hydrogen bonding and the hydrophobic effect, each amide molecule is involved in a cluster containing several water molecules thus rendering the water molecules less available to participate in the reaction and to solvate HS -and HCO 3 -ions. This effect is absent when ethylene glycol, a weak structure maker, is added, even in large quantities. This study clearly demonstrates the importance of solvent structure in the study of chemical reactions. State-of-the-art molecular dynamics simulations of the water-HMPA system could not reproduce the strongly negative excess volume of the mixture. This illustrates the need for more accurate force fields to simulate the structuring effect and their impact on chemical reactions.

Introduction

The effect of the solvent on reaction rates and chemical equilibrium has been studied for a long time [START_REF] Buncel | The Role of the Solvent in Chemical Reactions[END_REF]. The intrinsic properties of the solvent, namely, dipole moment, polarizability, dielectric constant, nucleophilic and electrophilic properties, among others, are not always sufficient to describe the solvent effects, rather, one needs to look at the complex physicochemical interactions of solvents, reactants, transition states, and products. Therefore, there are many ways solvents can influence the kinetic and thermodynamic properties of a reaction [START_REF] Varghese | Origins of complex solvent effects on chemical reactivity and computational tools to investigate them: a review[END_REF]. In aqueous solutions, hydrogen bonding as well as the hydrophobic effect, play a crucial role [START_REF] Breslow | Hydrophobic effects on simple organic reactions in water[END_REF], [START_REF] Sijbren | Hydrophobic interactions and chemical reactivity[END_REF]. They can result in the formation of clusters or complexes between water and co-solvent molecules [START_REF] Godhino | Use of Catalytic Thermometric Titrimetry To Investigate the Interaction of Dipolar Aprotic Solvents with Water[END_REF]. Their impact on chemical reactions has not been studied extensively and is not sufficiently understood. In this paper we report a case for which the addition of a strong liquid structure maker to an aqueous solvent significantly inhibits acid-base chemical reactions, even in the presence of excess water.

System and Experimental Methods:

In this study, the chemical absorption of CO 2 , a Lewis acid, and H 2 S, a Brønsted and Lewis acid, in an aqueous tertiary amine (Am) were chosen because they are well-known acid-base reactions in which, according to the literature, the role of water is different between CO 2 and H 2 S absorption [START_REF] Versteeg | On the kinetics between CO2 and alkanolamines both in aqueous and non-aqueous solutions-II. Tertiary amines[END_REF]. In the chemical absorption reaction of CO 2 , water participates via the reactions OH -+ CO 2 ↔ HCO 3 -, Am + CO 2 + H 2 O ↔ AmH + + HCO 3 -, and, as a solvation medium for the ions. In the chemical absorption reaction of H 2 S, water mainly solvates the ions, and does not directly participate in the direct and very fast proton transfer via Am + H 2 S ↔ AmH + + HS -, but via OH -+ H 2 S ↔ H 2 O + HS -. Following solvents have been studied (see SI for suppliers and purity): aqueous MDEA (MethylDiEthnolAmine), aqueous MDEA with physical cosolvent EG (Ethylene Glycol) and, aqueous MDEA with physical co-solvent HMPA (HexaMethylPhosphorTriAmide). Aqueous DEAE-EO (2-(2-Diethylaminoethoxy)ethanol), pure, or with EG or HMPA, has also been investigated (results presented in SI). Pure HMPA is a polar aprotic solvent with a very strong dipole moment µ = 5.54D, a dielectric constant ε=30 and a relatively low acceptor number (AN=10.6), making it a Lewis base, selectively solvating cations and Lewis acids [START_REF] Mayer | A semiempirical model for the description of solvent effects on chemical reactions[END_REF]. Pure EG is a polar protic solvent with a low dipole moment µ = 2.28D, a high dielectric constant ε=41.4 and a relatively high acceptor number (AN=44.9).

Vapor-Liquid Equilibria (VLE) were measured using the static-analytical method as described in ref. [START_REF] Orlov | Chemoinformatics-Driven Design of New Physical Solvents for Selective CO2 Absorption[END_REF]. CO 2 and H 2 S absorption rate were measured in a Lewis cell, according to the method described in ref. [START_REF] Orlov | Computational screening methodology identifies effective solvents for CO2 capture[END_REF]. For all tests, CO 2 and H 2 S absorption were measured simultaneously. In all experiments the amine concentration is 13 mol%. Viscosity measurements were carried out using an Anton Paar Lovis 2000 ME visco/densitometer based on the falling sphere principle. According to the manufacturer, the accuracy is 0.5% of the measured value. Figure 1 shows the VLE data for the absorption of CO 2 in aqueous MDEA. In the absence of HMPA there is a strong absorption of CO 2 because the CO 2 partial pressure in the gas phase is relatively low. When 9 mol% HMPA is added to the solvent (HMPA/water (1:9) molar ratio), there is a significant shift of the absorption isotherm to the left, in comparison to the aqueous MDEA isotherm, indicating a much weaker absorption of CO 2 , even if for every HMPA molecule there are still 9 H 2 O molecules in the solvent. Replacing half of the water molecules by HMPA (HMPA/water (1:1) molar ratio) further shifts the isotherm to the left, but the shift is smaller. Replacing all the water by HMPA (MDEA/HMPA isotherm), results in an isotherm that is almost identical to the isotherm of the physical CO 2 absorption in pure HMPA. The MDEA/HMPA and pure HMPA isotherms are close to and in between the isotherms for 9 mol% HMPA (MDEA/HMPA/water (1:9)) and 43.5 mol% HMPA (MDEA/HMPA/water (1:1)). This indicates that even with a large excess of water (9 water molecules per HMPA molecule), the CO 2 absorption is close to physical absorption in non-aqueous MDEA/HMPA. The complete absence of water even favors the physical absorption of CO 2 (water is a bad physical solvent for CO 2 ). This explains why the 43.5 mol% HMPA isotherm (MDEA/HMPA/water (1:1)) is left of the two nonaqueous solvents. Interestingly, when half of the water molecules is replaced by ethylene glycol (EG), the shift of the isotherm to the left is less pronounced than when only 1 out of 9 molecules is replaced by a HMPA molecule.

In Figure 2 the results are given for H 2 S absorption. Likewise for CO 2 , the impact of replacing one out of nine water molecules by HMPA results in a stronger shift to the left of the isotherm than when half of the water molecules are replaced by ethylene glycol. The isotherms for 43.5 mol% HMPA, and the nonaqueous isotherms are all linear, indicating a physical absorption of the H 2 S. HMPA is reported be an excellent solvent for physical H 2 S absorption [START_REF] Orlov | Computer-aided design of new physical solvents for hydrogen sulfide absorption[END_REF]. This is confirmed here by the low partial pressures of H 2 S at higher H 2 S mole fraction. At lower H 2 S mole fractions the chemical solvents absorb more H 2 S (see Figures at very low partial pressures in SI). Note that experimental measurements at very low H 2 S partial pressures are challenging.

In the SI the same figures are given, with the loading of the solvent expressed in mol of acid gas per mol of amine. The observations are unchanged.

In Figure 3 the absorption of CO 2 is shown as a function of time. Replacing half of the water molecules by ethylene glycol has barely an impact on the overall absorption rate of CO 2 . However, the CO 2 absorption rate slows down significantly when only 1 out of 9 water molecules is replaced by a HMPA molecule. In Figure 4 shows the results for H 2 S, for which the observation is very similar to CO 2 . Note that, as expected, H 2 S absorbs much faster than CO 2 , because it is a rapid proton transfer. The kinematic viscosities of both the EG and HMPA containing solvents are very similar and thus don't explain their different impact on the CO 2 absorption rate. The results for DEAE-EO (pKa=10,15) are given in the SI and are similar to MDEA (pKa= 8,65), indicating that they don't depend on the reactivity of the tertiary amine. The pH of the solvents has also been measured at ~323 K and is 11,7 for aqueous MDEA (13 mol%), 10,0 for aqueous MDEA (13 mol%) with 9 mol% HMPA and 10,6 for aqueous MDEA (13 mol%) with ethylene glycol (42,5 mol%). The addition of HMPA thus lowers the basicity of the solvent more than the addition of ethylene glycol. The change in basicity impacts the acid gas absorption properties.

Discussion

To explain the significant difference in impact between HMPA and ethylene glycol one needs to include their interactions with water in the analysis. In Figure 6 the reduced excess molar volumes [START_REF] Desnoyers | Treatment of Excess Thermodynamic Quantities for Liquid Mixtures[END_REF] for the binary mixtures HMPA-water [START_REF] Ivanov | Thermodynamics of (water + hexamethylphosphoramide) mixtures: Heat capacity properties in the temperature range between 283[END_REF] and EG-water [START_REF] Tsierkezos | Thermodynamic Properties of Water + Ethylene Glycol at 283[END_REF] are compared. The addition of HMPA molecules to water quickly results in a strongly non-ideal behavior (strongly negative reduced excess volume), indicating a strong structuring effect of HMPA due to the formation of HMPA-water clusters [START_REF] Desnoyers | Treatment of Excess Thermodynamic Quantities for Liquid Mixtures[END_REF], [START_REF] Ivanov | Thermodynamics of (water + hexamethylphosphoramide) mixtures: Heat capacity properties in the temperature range between 283[END_REF], [START_REF] Redlich | Algebraic representation of the thermodynamic properties and the classification of solutions[END_REF]. Note that this corre-sponds to a destruction of the water structure. Similar observations have been made for other aprotic polar solvents like NMP, DMI, and so on [START_REF] Uosaki | Excess molar volumes of (cyclic amide + water) at 298.15 and 308.15K[END_REF], [START_REF] Lazarev | Structure of hydration shell of 1,3-dimethyl-2-imidazolidinone according to FT -Raman spectroscopy[END_REF]. It has been reported that such molecules form stronger hydrogen bond clusters with several water molecules and are therefore structure makers. The water-amide hydrogen bonds are stronger than those between water molecules and their strength increases with increasing donor number (lower acceptor number) of amides [START_REF] Ivanov | Thermodynamics of (water + hexamethylphosphoramide) mixtures: Heat capacity properties in the temperature range between 283[END_REF], [START_REF] Zaichikov | Determination of the intermolecular interaction parameters in the water-amide systems based on the data of the excess thermodynamic functions[END_REF], [START_REF] Zaichikov | yaninov, Structural thermodynamic characteristics and intermolecular interactions in aqueous solutions of hexamethylphosphortriamide[END_REF]. Despite steric effects, the water-HMPA hydrogen bonds should thus be much stronger than the waterwater hydrogen bonds. A Prigogine-Flory-Patterson analysis of the HMPA-water mixture (see SI) confirms that the overall negative excess volume originates mainly from the intermolecular interactions, and to some extend from the free volume. The internal pressure contribution is positive. Moreover, HMPA shows a pronounced hydrophobicity and the hydrophobic effect will also significantly contribute to non-ideal behavior. Interestingly, the acceptor number of the HMPAwater mixture decreases almost linearly from 54.8 (pure water) to 10.6 (pure HMPA) as a function of HMPA content [START_REF] Mayer | NMR-Spectroscopic Studies on Solvent Electrophilic Properties, Part II: Binary Aqueous-Non Aqueous Solvent Systems[END_REF]. The non-ideality of the mixture in terms of excess volume, excess heat capacity and acid gas absorption, thus doesn't seem to be captured by the acceptor number (Lewis basicity). The non-ideality of the EG-water mixture is much less pronounced: EG is only a weak structure maker [START_REF] Tsierkezos | Thermodynamic Properties of Water + Ethylene Glycol at 283[END_REF]. Similarly, for the MDEA-water system the excess volumes are only slightly negative [START_REF] Bernal-Garcia | Densities and Excess Molar Volumes of Aqueous Solutions of n-Methyldiethanolamine[END_REF]. Figure 6. Reduced excess molar volume of EG-water [START_REF] Tsierkezos | Thermodynamic Properties of Water + Ethylene Glycol at 283[END_REF] and HMPA-water [START_REF] Ivanov | Thermodynamics of (water + hexamethylphosphoramide) mixtures: Heat capacity properties in the temperature range between 283[END_REF] mixtures as function of the EG or HMPA molar fraction x 1 , at 293,15 K and at 303,15 K. Circles are experimental datapoints, the full lines are the Redlich-Kister equations [START_REF] Redlich | Algebraic representation of the thermodynamic properties and the classification of solutions[END_REF].

The acceptor number of the HMPA-water mixture decreases from 54.8 (pure water) to 10.6 (pure HMPA) as a function of HMPA content [START_REF] Mayer | NMR-Spectroscopic Studies on Solvent Electrophilic Properties, Part II: Binary Aqueous-Non Aqueous Solvent Systems[END_REF] indicating the mixture becomes gradually more basic with the addition of HMPA and thus that the activation energy for acid-base reactions with amines should change [START_REF] Mayer | A semiempirical model for the description of solvent effects on chemical reactions[END_REF]. However, in the presence of a tertiary alkanolamine like MDEA, a Lewis base, measurements of the pH show a decrease by ~1,7 units after the addition of 9 mol% HMPA (Lewis base), indicating a lower overall basicity. Water, which plays the role of the acid, thus seems to have a lower activity, probably due to its strong hydrogen bonds with HMPA. Indeed, because of the preferential H-bonding between water and HMPA, a significant part of the water is only present passively, thus not participating in any acid-base reaction. In thermodynamics this non-ideality is modelled by an activity that represents the active concentration of each solvent component available to participate in the chemical reaction, as a reactant or as a product, or as a solvation medium [START_REF] Poling | The properties of gases and liquids[END_REF]. The activity coefficient of water in the aqueous HMPA-amine solvent is thus lower than one (the activity of water is lower than expected from its concentration). The equilibrium of the dissociation reaction of water into H + and OH -is thus shifted towards less dissociation, resulting in a lower basicity of the solvent. In a previous analysis using Hansen solubility parameters Wanderly et al. didn't include the hydrophobic effect in his analysis and suggested that the solvent hydrogen-bonding might not play an important role in CO 2 absorption [START_REF] Wanderley | CO2 solubility and mass transfer in water-lean solvents[END_REF] This effect is almost absent in the presence of ethylene glycol. The gradual replacement of water by ethylene glycol results in a linear reduction of CO 2 absorption. Moreover, due to its hydroxyl groups, EG can even stabilize the reaction products, although weaker than water.

To better understand and visualize the structuring effect of HMPA, state-of the art molecular dynamics simulations were performed the aqueous HMPA mixture (see SI). The simulated temperature dependency of pure water and pure HMPA agree well with experimental results. The non-ideal experimental behavior of the density of the mixture of water and HMPA could not be reproduced by the molecular dynamics simulations. The strong experimental increase in density (negative excess volume) when a small amount of HMPA is added to water (structuring effect), for example, is not captured by the atomistic model, which can thus not be used to gain more insight in the molecular mechanisms. This demonstrates the need for more accurate force fields to simulate the structuring effect of solvent molecules and their impact on chemical reactions.

Conclusion:

The lower thermodynamic absorption capacity as well as the slower absorption rate of CO 2 , and, to a lesser extent H 2 S, in an aqueous tertiary amine solution in the presence of excess water and limited amount of a strong structure maker such as HMPA can be explained by the strong nonideality between water and HMPA. The formation of clusters of HMPA molecules and several water molecules due to strong hydrogen bonding results in water, although present in the solvent, being unavailable to participate in the chemical reaction (specifically for CO 2 absorption) and to stabilize the products of the absorption reactions, mainly HCO 3 -and HS -. This effect dominates the other intrinsic properties of the system. The impact of EG is much weaker.

ASSOCIATED CONTENT

Supporting Information. The supporting information section includes information on the chemicals used, experimental VLE data for aqueous MDEA and aqueous DEAE-EO with and without co-solvents, kinetic data for DEAE-EO, viscosity of DEAE-EO and the molecular simulations. This material is available free of charge via the Internet at http://pubs.acs.org.

Figure 1 .

 1 Figure 1. Experimental vapor-liquid equilibrium of CO 2 absorption at 323 K. (a:b) is the rounded molar ratio a/b of the last two components. The amine concentration is always 13 mol%. The solid lines are a guide to the eye.Figure1shows the VLE data for the absorption of CO 2 in aqueous MDEA. In the absence of HMPA there is a strong absorption of CO 2 because the CO 2 partial pressure in the gas phase is relatively low. When 9 mol% HMPA is added to the solvent (HMPA/water (1:9) molar ratio), there is a significant shift of the absorption isotherm to the left, in comparison to the aqueous MDEA isotherm, indicating a much weaker absorption of CO 2 , even if for every HMPA molecule there are still 9 H 2 O molecules in the solvent. Replacing half of the water molecules by HMPA (HMPA/water (1:1) molar ratio) further shifts the isotherm to the left, but the shift is smaller. Replacing all the water by HMPA (MDEA/HMPA isotherm), results in an isotherm that is almost identical to the isotherm of the physical CO 2 absorption in pure HMPA. The MDEA/HMPA and pure HMPA isotherms are close to and in between the isotherms for 9 mol% HMPA (MDEA/HMPA/water (1:9)) and 43.5 mol% HMPA (MDEA/HMPA/water (1:1)). This indicates that even with a large excess of water (9 water molecules per HMPA molecule), the CO 2 absorption is close to physical absorption in non-aqueous MDEA/HMPA. The complete absence of water even favors the physical absorption of CO 2 (water is a bad physical solvent for CO 2 ). This explains why the 43.5 mol% HMPA isotherm (MDEA/HMPA/water (1:1)) is left of the two nonaqueous solvents. Interestingly, when half of the water molecules is replaced by ethylene glycol (EG), the shift of the isotherm to the left is less pronounced than when only 1 out of 9 molecules is replaced by a HMPA molecule.In Figure2the results are given for H 2 S absorption.

Figure 2 .

 2 Figure 2. Experimental vapor-liquid equilibrium data points for H 2 S absorption at 323 K. (a:b) is the molar ratio a/b of the last two components. The amine concentration is always 13 mol%. The solid lines are a guide to the eye.

Figure 3 .

 3 Figure 3. Time-evolution of the fraction of CO 2 in the solvent relative to its fraction at equilibrium. (a:b) is the molar ratio a/b of the last two components. The amine concentration is always 13 mol%. T=323,15 K. The solid lines are a guide to the eye.

Figure 4 .

 4 Figure 4. Time-evolution of the fraction of H 2 S in the solvent relative to its fraction at equilibrium. (a:b) is the molar ratio a/b of the last two components. The amine concentration is always 13 mol%. T=323,15 K. The solid lines are a guide to the eye.The kinematic viscosities of the three solvents are shown in Figure5.

Figure 5 .

 5 Figure 5. Kinematic viscosity as a function of temperature. (a:b) is the molar ratio a/b of the last two components. The amine concentration is always 13 mol%.

ABBREVIATIONS

MDEA MethylDiEthanolAmine; DEAE-EO 2-[2-'Diethylamino) ethoxy] -ethanol; HMPA HexaMethylPhosphorAmide; EG Ethylene Glycol; NMP N-MethylPyrrolidone; DMI 1,3-Dimethyl-2imidazolidinone; AG Acid Gas; AN Acceptor Number; DN Donor Number; Am Amine.

Funding Sources

ACKNOWLEDGMENT

This work was supported by the CO 2 and Sustainability R&D program of TotalEnergies S.E.

AUTHOR INFORMATION

Corresponding Author * frederick.de-meyer@totalenergies.com Present Addresses †Université de Toulouse, IMT Mines Albi, CNRS UMR 5302, Centre Rapsodee Campus Jarlard, 81013 Albi CT Cedex 9, France

Author Contributions

The manuscript was written through contributions of all authors.

All authors have given approval to the final version of the manuscript.