Bertrand Kirsch
email: bertrand.kirsch@ecole-air.fr

Olivier Montagnier

Emmanuel Bénard

Thierry M Faure

Tightly coupled aeroelastic model implementation dedicated to fast aeroelastic tailoring optimisation of high aspect ratio composite wing

Keywords: aeroelasticity, composite materials, reduced order model, aeroelastic tailoring, High Altitude Pseudo-Satellite, drones

come

Introduction

The growing efforts made in the fields of solar cells, energy storage and composite materials pave the way to a new concept of aircraft, namely High Altitude Pseudo Satellite (HAPS). Alongside stratospheric airship, a particular type of solar powered High Altitude Long Endurance (HALE) fixed wing Unmanned Aerial Vehicles (UAV) aims to meet a virtually infinite endurance. Examples of HAPS projects, under progress or discontinued, are sketched on Fig. 1. To achieve this ambitious goal, because of the low on-board power, aerodynamic and structural performances are stretched to their limits. It implies, on the aerodynamic side, a high aspect ratio wing optimising the lift-to-drag ratio and, on the structural side, a lightweight composite airframe. This results in a very flexible airframe, particularly vulnerable to destructive fluid/structure interactions like torsional divergence and flutter. A review on non-linear aeroelasticity of high aspect-ratio wings can be found in [START_REF] Afonso | A review on non-linear aeroelasticity of high aspect-ratio wings[END_REF]. Classical solutions designed to push further aeroelastic critical speed mostly rely on the stiffening of the airframe or the adjustment of mass distribution. Both options are detrimental to mass balance which is a key feature of HAPS. In that context, alternative solutions should be explored, among these are aeroelastic tailoring discussed below.

This technology, mentioned for the first time by Munk concerning propeller with optimised wood fibres orientation [START_REF] Munk | Propeller containing diagonally disposed fibrous material[END_REF], has moved ahead with the forward-swept wing X-29 research aircraft, prone to torsional divergence Figure 1: Sketch of some HAPS projects with a comparison against Airbus A320-200 dimensions [START_REF] Weisshaar | Aeroelastic tailoring of forward swept composite wings[END_REF]. Aeroelastic tailoring consists in exploiting laminate composite anisotropy by setting a proper layup. For most applications, a laminate layup presents a mirror symmetry with respect to its middle plane. If we consider a laminate plate of unidirectional (UD) plies, the layup [-45°, 0°, 45°, 90°, 90°, 45°, 0°, -45°], with angular values representing fibre orientation in each ply, presents a mirror symmetry [START_REF] Hahn | Introduction to Composite Materials[END_REF]. The aim of this type of symmetry is to dissociate membrane behaviour (inplane loadings imply inplane displacements and vice versa) from bending behaviour. Another common rule concerning the UD plies orientation is to set balanced layup, namely to put for instance as many -45°plies as 45°plies [START_REF] Hahn | Introduction to Composite Materials[END_REF]. It is meant to avoid coupling between the different bending behaviour. The true principle of aeroelastic tailoring is to ignore these rules, thus allowing coupling between the different behaviours of the laminate. It mainly consists in creating a link between the bending and the warping of the laminate. On the aerodynamic side, it induces the coupling of the bending due to lift forces and the twisting of the wing which determines the local angle of attack (AoA). Thus, the aeroelastic tailoring is a way of establishing a feedback loop on the aeroelastic behaviour of a very flexible aircraft. This phenomenon has been specifically studied on the Thin Wall Beam (TWB) composite configuration, for instance in the recent work of Farsadi et al. [START_REF] Farsadi | Geometrically nonlinear aeroelastic behavior of pretwisted composite wings modeled as thin walled beams[END_REF]. We could also mention studies conducted by Qin and Librescu [START_REF] Qin | Aeroelastic instability of aircraft wings modelled as anisotropic composite thin-walled beams in incompressible flow[END_REF] or Haddadpour and Zamani [START_REF] Haddadpour | Curvilinear fiber optimization tools for aeroelastic design of composite wings[END_REF] concerning the aeroelastic design of a TWB composite wing and research conducted on tow steered composites by Stodiek et al. [? ]. A more general review of the ongoing efforts concerning the aeroelastic tailoring technology can be found in [START_REF] Jutte | Aeroelastic tailoring of transport aircraft wings: State-of-the-art and potential enabling technologies[END_REF].

The difficulty lays on the modelling of this particular behaviour. In fact, different approximations can be made depending on the motion of the structure relatively to the fluid. In order to quantify this phenomenon, dimensionless numbers are used. In our case, the degree of fluid/structure coupling is characterised by the reduced frequency f r = T f /T s with T f the period of fluid motion and T s the period of structure motion [START_REF] Hemon | Vibrations des structures couplees avec le vent[END_REF]. Because of their low structural frequencies and their low flight speed, solar powered HALE UAV are characterised by a reduced frequency near 1. In fact, it means that both dynamic are tightly linked. Therefore, high fidelity simulation of this kind of fluid structure/interaction with a finite volume method coupled with a finite element method leads to a prohibitive computational cost. Consequently, aeroelastic reduced order model are still widely used, especially for Very Flexible Aircraft (VFA). The large displacement and rotation encountered by a very flexible wing associated with a low flight speed, make the aeroelastic behaviour all the more difficult to simulate, because of the induced geometrical non linearities and the fluid motion unsteadiness.

A relevant illustration of the need for an accurate modelling of theses phenomena is the accident of the Helios UAV which occurred on the 26th June 2003 due to a large deflection and rotation of the wing leading to an unstable pitch oscillation [START_REF] Noll | Investigation of the Helios prototype aircraft mishap[END_REF]. To meet this need, several reduced order aeroelastic models have been developed. For computational efficiency, most of them are based on inviscid, incompressible potential flow theory coupled with beam or plate theories. We could mention computation code NANSI (Nonlinear-Aerodynamics/ Nonlinear-Structure Interaction) [START_REF] Wang | Nonlinear-aerodynamics/nonlinear-structure interaction methodology for a high-altitude long-endurance wing[END_REF] which combines an Unsteady Vortex Lattice Method (UVLM) and a nonlinear beam theory. The UVLM is particularly useful in case of low aspect ratio wing or delta wing, because the method is able to predict 3D effects. Another solution is proposed by Murua in SHARP program (Simulation of High Aspect Ratio Planes) [START_REF] Murua | Assessment of wake-tail interference effects on the dynamics of flexible aircraft[END_REF] using UVLM with a displacement based geometrically exact beam theory. Some models are dedicated to high aspect ratio wings like Drela's program ASWING [START_REF] Drela | Integrated simulation model for preliminary aerodynamic, structural, and control-law design of aircraft[END_REF]. This VFA conception tool combines a nonlinear isotropic beam formulation with an unsteady lifting line theory. More recently, Shearer and Cesnik have developed a Matlab toolbox called UM/NAST (University of Michigan/ Nonlinear Aeroelastic Simulation Toolbox) [START_REF] Shearer | Nonlinear flight dynamics of very flexible aircraft[END_REF] made up of a strain-based geometrically nonlinear beam formulation linked with a finite state two-dimensional incompressible flow aerodynamic theory proposed by Peters et al. [START_REF] Peters | Finite state induced flow models. I-Two-dimensional thin airfoil[END_REF]. A similar formulation is used by Ribeiro in the Matlab toolbox Aeroflex [START_REF] Ribeiro | Aeroflex: a toolbox for studying the flight dynamics of highly flexible airplanes[END_REF]. One last example is the Matlab toolbox proposed by Patil and Hodges called NATASHA (Nonlinear Aeroelastic Trim and Stability of HALE Aircraft) [START_REF] Chang | Flight dynamics of highly flexible aircraft[END_REF] coupling an intrinsic beam formulation with Peters' theory.

The present paper deals with the development of a code well suited for VFA aeroelastic tailoring and the evaluation of aeroelastic tailoring effect on critical speeds. The number of configurations to be studied to optimise the laminate stacking sequence of an anisotropic wing box requires an extremely low computation cost and therefore a very efficient code. For this purpose, the choice for structural computation fell on an open source tool in Fortran 90/95 named GEBT (Geometrically Exact Beam Theory) developed by Yu and Blair [START_REF] Yu | GEBT: A general-purpose nonlinear analysis tool for composite beams[END_REF] and Wang and Yu [START_REF] Wang | Geometric Nonlinear Analysis of Composite Beams Using Wiener-Milenkovic Parameters[END_REF]. Since this code does not integrate the aerodynamic forces, the first part deals with the construction of the tight coupling between the GEBT and the unsteady model of Peter's. In case of complex cross section shape or anisotropic composite materials, flexibility matrix is determined using an homogenisation tool consisting in a three-dimensional finite element calculation using periodic boundary conditions (unique element spanwise) in CalculiX. The second part focuses on the actual calculation of the critical speeds in modal and temporal domain. A particular focus is put on the fast critical speed computation strategy using a non-iterative modal approach about the geometrically non linear deformed shape of the wing with the computation of only a few aeroelastic modes. A Python pre/postprocessor is added to take into account, among others, beam homogenisation step and critical speed computation, and to permit a smooth integration into a multidisciplinary optimisation framework like openMDAO [START_REF] Gray | OpenMDAO: an open-source framework for multidisciplinary design, analysis, and optimization[END_REF]. The resulting open source computation code is called "GEBTAero". In the third part, computation speed and accuracy of this implementation is assessed using widely used aeroelastic test cases and compared to other aeroelastic codes. Finally, since no common experimental or numerical test cases have been found in the literature, the code is used to study aeroelastic tailoring. For this goal, a representative high aspect ratio composite wing test case derived from Patil wing is proposed first. The impact of aeroelastic tailoring is then investigate on a simple two-ply composite laminates by sweeping both ply orientations (lower ply θ 1 and upper ply θ 2 ).

Medium-fidelity nonlinear aeroelastic model

The high aspect ratio assumption gives us the opportunity of neglecting tridimensional effects and thus using a strip theory which can be easily linked to a beam formulation. This implementation can be done in two different ways:

A loose coupling consisting in defining the aerodynamic loads to apply on the beam in accordance with position and speed parameters extracted from the last structural calculation. This method is easy to implement and is well suited for modular architecture like, for instance, the wind turbines conception software FAST (Fatigue, Aerodynamics, Structures, and Turbulence) from the National Renewable Energy Laboratory (NREL) [START_REF] Jonkman | FAST User s Guide[END_REF]. However, it has also drawbacks like the inability of performing coupled eigenvalue analysis or some convergence issues for nonlinear Newton-Raphson algorithms (aerodynamic loads are not taken into account into the Jacobian matrix). It also implies in our ranges of reduced frequencies a very small time step.

A tight coupling is realised by integrating aerodynamic loads directly into the weak formulation of the beam theory. This approach is much more complex and less modular but more powerful. Indeed, the most interesting application is the possibility to determine, for a particular flow velocity, the aeroelastic modes of the wing about a steady sate, namely frequencies, modal shapes and damping factors. The latter is a key parameter for our study because it defines the limit between stable and unstable speed and thereby provides the flutter boundary. For all these reasons, tight coupling has been chosen for our toolbox.

Geometrically exact beam theory

Due to aeroelastic tailoring use and airframe flexibility, it is essential to ensure a proper modelling of the laminate anisotropy and geometrical nonlinearities. For this purpose, the choice fell on an open source tool named GEBT (Geometrically Exact Beam Theory) [START_REF] Yu | GEBT: A general-purpose nonlinear analysis tool for composite beams[END_REF][START_REF] Wang | Geometric Nonlinear Analysis of Composite Beams Using Wiener-Milenkovic Parameters[END_REF] designed for composite slender structures under large deflections and rotations, assuming the strains to be small. This tool coded in Fortran 90/95 implements a mixed variational formulation based on exact intrinsic equations for dynamics of moving beams developed by Hodges [START_REF] Hodges | A mixed variational formulation based on exact intrinsic equations for dynamics of moving beams[END_REF].

The exact intrinsic equations for dynamics are derived from Hamilton's weak principle asymptotically developed along the beam axes:

t 2 t 1 L 0 δ (K -U ) + δW d x 1 d t = δA (1)
where t 1 and t 2 are arbitrary fixed times, and U are the kinetic and strain energy, respectively, δ is the usual Lagrangian variation for a fixed time, δW is the virtual work of applied loads and δA the virtual action on the same period. The resulting formulation is detailed in Hodges [START_REF] Hodges | A mixed variational formulation based on exact intrinsic equations for dynamics of moving beams[END_REF].

The main strength of this method compared to classical displacement based formulation is to avoid the dependency from a coordinate system (intrinsic nature) for the position and rotation parameters. Then, kinematical and constitutive relations are added to the weak formulation with Lagrange multipliers (mixed nature). The resulting formulation allow a finite element implementation with very simple shape functions (constant or linear). According to Hodges [START_REF] Hodges | Finite Element Solution of Nonlinear Intrinsic Equations for Curved Composite Beams[END_REF], we defined three coordinates systems (Fig. 2):

• a unique global body attached frame a x a , y a , z a moving with a given linear and angular velocity v a and ω a in an inertial frame and consistent with flight mechanics conventions (x a pointing towards, y a pointing the right wing and z a pointing downwards).

• at least one undeformed beam frame b x b , y b , z b fixed in frame a : x b is normal to the airfoil and points towards the right. x b is within the plane y a , z a even in case of a swept wing, y b is tangent to the airfoil reference line pointing forward and z b completes the triad. In our case, a frame b is defined for each section of the wing with a different dihedral or/and wing twist.

• a deformed beam frame B x B , y B , z B for each beam element.

Direction cosine matrix describing the rotation between frames are defined using Rodrigues parameters θ i = 2e i tan(α/2) with α the magnitude of the rotation about a unit vector e :

C = 1 -(1/4)θ T θ ∆ -θ + (1/2)θ T θ 1 + (1/4)θ T θ (2)
with ∆ the identity matrix and the tilde notation defining a matrix using the following identity for any vector w :

θ w = θ ∧ w (3)
Fundamentals unknowns of this formulation are the displacement u a and the Rodrigues parameters θ a in frame a , the internal forces and moments F B and M B in frame B and the linear and angular momenta P B and H B in frame B . These unknowns are linked to the weak formulation with the following expressions: The anisotropic nature of the beam is concentrated in Eq. ( 5). I and K matrices are the cross sectional properties of the beam and could be determined by various means. Proprietary program Variational Beam Sectional Analysis (VABS) [START_REF] Cesnik | VABS: a new concept for composite rotor blade cross-sectional modeling[END_REF] is used in [START_REF] Yu | GEBT: A general-purpose nonlinear analysis tool for composite beams[END_REF]. The present work uses a 3D finite element reduction done with the open source solver CalculiX [START_REF] Dhondt | Calculix: a free software three-dimensional structural finite Element Program[END_REF] detailed in section 2.2. Eqs (4)-( 5) are expanded below in frame B :

P B H B =      ∂ K ∂ V B ∂ K ∂ Ω B      = I V B Ω B ( 4
)
F B M B =      ∂ U ∂ V B ∂ U ∂ Ω B      = K γ κ ( 5
)
             γ 11 2γ 12 2γ 13 κ 1 κ 2 κ 3              =        S 11 S
                    F 1 F 2 F 3 M 1 M 2 M 3              (6)              P 1 P 2 P 3 H 1 H 2 H 3              =        µ 0 0 0 µx m3 -µx m2 0 µ 0 -µx m3 0 0 0 0 µ µx m2 0 0 0 -µx m 3 µx m2 i 22 + i 33 0 0 µx m3 0 0 0 i 22 -i 23 -µx m 2 0 0 0 -i 23 i 33                     V 1 V 2 V 3 Ω 1 Ω 2 Ω 3              (7)
with S i j the coefficients of the flexibility matrix (S = K -1 ), µ the mass per unit length, x m 2 , x m3 the coordinates of the mass center about respectively y B and z B , i 22 the mass moment of inertia about y B , i 33 the mass moment of inertia about z B and i 23 the product of inertia in the plane y B , z B . Dividing a beam into N elements with the starting node of the i th element numbered as i and the ending node numbered as i + 1 and using linear or constant shape function for the unknowns, the initial weak formu-lation leads to:

N i =1 δu T i f - u i + δu T i +1 f + u i + δψ T i f - ψ i + δψ T i +1 f + ψ i + δF T i f - F i _ + δF T i +1 f + F i + δM T i f - M i + δM T i +1 f + M i + δP T i f P i + δH T i f H i = δu T N +1 F N +1 + δψ T N +1 M N +1 -δF T N +1 û N +1 -δF T N +1 θ N +1 -δF T N +1 û N +1 -δu T 1 F 1 -δψ T 1 M 1 + δF T 1 û 1 + δM T 1 θ 1 (8)
with δu the virtual displacement, δψ the virtual rotation, () the boundary conditions and:

f ± u i = ±C T C a b F i -f ± i + ∆L i 2 ωa C T C a b P i + Ċ T C a b P i (9) f ± ψ i = ±C T C a b M i -m ± i + ∆L i 2 ωa C T C a b H i + Ċ T C a b H i + C T C a b Ṽ i ⃗ P i -ẽ 1 + γi F i ( 10
)
f ± F i = ± u i - ∆L i 2 C T C a b e 1 + γ i -C a b e 1 ( 11
)
f ± M i = ± θ i - ∆L i 2 ∆ + θ 2 + θ i θ T i 4 C a b κ i ( 12
)
f P i =C T C a b V i -v i -ωa u i -u i ( 13
)
f H i =Ω i -C b a C ω a -C b a ∆ -θ i /2 1 + θ i θ T i /4 θ i ( 14
)
with e 1 = 1 0 0 T , C the matrix direction cosines between frame b and B , C a b the matrix direction cosines between frame b and a , C b a = C a b T , ∆L i the length of the i th element and f ± i , m ± i the discretised distributed forces and moments defined by:

f - i = 1 0 (1 -ξ) f a ∆L i d ξ; f + i = 1 0 ξ f a ∆L i d ξ ( 15
)
m - i = 1 0 (1 -ξ) m a ∆L i d ξ; m + i = 1 0 ξm a ∆L i d ξ (16)
with f a and m a the distributed forces and moments. From Eqs ( 9)-( 16), subscripts a , b or B referred to the development frame. Derivations () are made in the inertial frame. The resulting nonlinear system of 18N + 12 equations consists in:

• 12 (N -1) equations associated to intermediate nodes:

f + u i + f - u i +1 = 0; f + ψ i + f - ψ i +1 = 0; f + F i + f - F i +1 = 0; f + M i + f - M i +1 = 0
• 12 equations associated to the starting node:

f - u 1 -F 1 = 0; f - ψ 1 -M 1 = 0; f - F 1 -û 1 = 0; f - M 1 -θ 1 = 0
• 12 equations associated to the ending node:

f + u N + F N +1 = 0; f - ψ 1 + M 1 = 0; f - F 1 + û 1 = 0; f - M 1 + θ 1 = 0
• 6N equations associated to the elements:

f P i = 0; f H i = 0
The corresponding 18N + 12 unknowns are u i , θ i , F i , M i for each node and P i and H i for each element. Boundary conditions are prescribed to boundary nodes, either in displacement/rotation or in force/moment. This formulation is so far purely structural and corresponds to the implementation of the structural solver GEBT [START_REF] Yu | GEBT: A general-purpose nonlinear analysis tool for composite beams[END_REF]. Based on this work, our toolbox will implement unsteady aerodynamic loads throughout distributed forces and moment f a and m a .

Composite beam homogenisation

In our toolbox, the flexibility matrix can be determined using different methods. In the case of isotropic Euler-Bernoulli beam (shear strains γ 12 and γ 13 are neglected, thus line and column 2 and 3 are set to 0), the flexibility matrix (Eq. ( 6)) is greatly simplified:

S =        1/ (E A) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1/ (G J ) 0 0 0 0 0 0 1/ (E I G 2 ) 0 0 0 0 0 0 1/ (E I G 3 )        (17)
with E A the axial stiffness, G J the torsional stiffness, E I G 2 the bending stiffness spanwise and E I G 3 the bending stiffness chordwise. These parameters could be:

• directly derived from test cases data;

• calculated for simple cross section shape (E , G are the Young and Coulomb modulus of the material, A is the cross section area, J is the torsion constant and I G 2,3 are the second moment of area).

In case of complex cross section shape or anisotropic composite materials, flexibility matrix is determined using an homogenisation tool following a method detailled in [START_REF] Cartraud | Computational homogenization of periodic beam-like structures[END_REF] and applied for instance in the Abaqus preprocessor HOMTOOLS [START_REF] Lejeunes | Une Toolbox Abaqus pour le calcul de propriétés effectives de milieux hétérogenes, 10e colloque national en calcul des structures[END_REF]. As explained in [START_REF] Kirsch | Assessment of aeroelastic tailoring effect on high-aspect-ratio composite wing flutter speed using an open source reduced order model solver[END_REF], It consists in a three-dimensional finite element calculation realised on a Representative Volume Element (RVE) of the beam using periodic boundary conditions along beam axis direction x B . The RVE is a 3D mesh written in Abaqus input format. Different cases are considered:

• simple shape constant cross section (plate or box): the mesh is automatically generated using the pre/postprocessor CalculiX GraphiX (Fig. 3a);

• constant cross section: a 2D mesh is extruded with a unique element spanwise (Fig. 3b);

• periodic cross section: a wing section representing a beam period is meshed (between two ribs for instance).

The frame used for the mesh is the frame B (Fig. 2). Two dummy nodes are created representing the beam strains and curvatures degrees of freedom (one with γ 11 and another with κ 1 , κ 2 and κ 3 ). The beam kinematic choice fell on an Euler-Bernoulli hypothesis: shear strain is neglected (line/column 2 and 3 are set to 0), the middle plane of the cross section is normal to the beam axis x B but warping is nevertheless allowed (Fig. 4). The flexibility matrix S becomes: The different materials (isotropic or orthotropic) of the cross section are then affected to mesh cells.

S =       
The idea behind this method is to consider the real beam kinematics as the superposition of an Euler-Bernoulli kinematics and a periodic motion of the RVE faces. This kinematic is imposed using a set of periodic conditions linking strains and curvatures (assimilated to dummy nodes displacements) to the displacements of nodes from opposite faces of the RVE in the beam axis direction x B ( -refers to the left face and + to the right face in the sense of increasing beam axis coordinate):

u + x -u - x = L x γ 11 + z + κ 2 -y + κ 3 (19) u + y -u - y = L x -z + κ 1 -x κ 3 ( 20
)
u + z -u - z = L x y + κ 1 -x κ 2 (21)
with u + x ,y ,z and u - x ,y ,z the nodes displacement of the two periodic faces, L x the RVE length along beam axis, x the coordinate of the middle of the RVE along the beam axis x B , y + and z + are the nodes coordinates within the right face.

Four elementary loading cases are imposed through dummy node forces/moments and a linear static calculation is performed for each case using the open source finite element solver CalculiX CrunchiX. For instance, the elementary torsional load case corresponds in Eq. 5 to M 1 = 1 and F 1 = M 2 = M 3 = 0. Then, flexibility matrix forth line or column is identified using the strains and curvatures vector (S 41 = γ 11 , S 44 = κ 1 , S 45 = κ 2 and S 46 = κ 3 ). By duality, these strains and curvatures correspond to the forces recorded on the two dummy nodes.

Furthermore, this reduction toolbox allow us to :

• determine the position of the elastic center EC. It corresponds to the reference point where no traction/bending coupling is generated. Thus :

y E C = S 61 S 66 (22)
z E C = S 51 S 55 (23)
• approximate the mass center position (x m2 , x m3 ), the mass moment of inertia i 22 , i 33 and the product of inertia i 23 of the mass matrix (Eq. 4) using the volume and the material density of each mesh cell.

Finite state induced flow unsteady aerodynamic

As required by the typical reduced frequencies of VFA, a two-dimensional unsteady aerodynamic model must be used [START_REF] Haddadpour | Evaluation of quasi-steady aerodynamic modeling for flutter prediction of aircraft wings in incompressible flow[END_REF]. Because of the low Mach number and the relatively high Reynolds number (order of magnitude of 10 6 ) the choice fell on inviscid potential flow theories. Many of them rely on Wagner formulation in time domain [START_REF] Wagner | Über die Entstehung des dynamischen Auftriebes von Tragflügeln[END_REF] and Theodorsen theory [START_REF] Theodorsen | General theory of aerodynamic instability and the mechanism of flutter[END_REF] in frequency domain. The latter is suitable for computing flutter speed but implies the use of iterative loop like the so called p-k method. Indeed, aerodynamic loads depends on reduced frequency through Theodorsen function which, in is turn, is a function of velocity. To compel with computation needs, a more suitable type of method has been developed derived from the latter namely finite state approximations [START_REF] Peters | Two-dimensional incompressible unsteady airfoil theory-an overview[END_REF]. Instead of modelling shed wake effect with Theodorsen function, finite state approximation computes it with a number N S of states resolved from a set of N S Ordinary Differential Equations (ODE) which could be solved at the same time as the non linear system of the beam formulation. The resultant equations gives convergence to both Theodorsen and Wagner functions relatively to N S . However, this formulation does not take into account stall phenomena and is thus not adapted to Limit Cycle Oscillation (LCO) simulation. Moreover, because of the incompressible hypothesis, Mach number should be reasonably limited to 0.3. In addition, even if full scale HAPS presents high Reynolds number, a particular attention should be taken when dealing with scale model. It is notably the case for wind tunnel specimens.

A widely used finite state approximation for VFA or rotor blade aeroelasticity is the method developed by Peters et al. [START_REF] Peters | Finite state induced flow models. I-Two-dimensional thin airfoil[END_REF]. This formulation is implemented in our toolbox with the following aerodynamic loads:

L = πρb 2 ḧ + U α -b a α + 2πρU b ḣ + U α + b 1 2 -a α -λ 0 ( 24
)
M = b 1 2 + a L -πρb 3 1 2 ḧ + u α + b 1 8 - a 2 α ( 25
)
with L the linear lift, M the linear moment around a reference point F , ρ the air density and U the flow velocity.

The semi-chord b , the height h , AoA α and the distance a between the point F and the semi-chord are detailed in Fig. 5. The induced-flow velocity λ 0 is approximated using N S induced-flow states λ 1 , λ 2 , . . . , λ N S by:

λ 0 ≈ 1 2 N S n=1 b n λ n
where the b n are found in [START_REF] Hodges | Introduction to structural dynamics and aeroelasticity[END_REF] by the least-square method. Then, the induced-flow dynamics are derived from the assumption that the shed vorticity stay in the plane of the airfoil and travel downstream with the same velocity as the flow. ⃗ λ is a column matrix containing the values of λ n determined using a set of N S first-order ODEs [START_REF] Hodges | Introduction to structural dynamics and aeroelasticity[END_REF]:

A λ + U b λ = ḧ + U θ + b 1 2 -a θ c ( 26
)
Matrix ⃗ ⃗ A is defined by:

A = D + d b T + c d T + 1 2 c b T ( 27
)
where

D n m =      1 2n n = m + 1 -1 2n n = m -1 0 n ̸ = m ± 1 ( 28
) b n =    (-1) n -1 (N S + n -1)! (N S -n -1)! 1 (n !) 2 n ̸ = N S (-1) n-1 n = N S ( 29
)
d n = 1 2 n = 1 0 n ̸ = 1 ( 30
)
c n = 2 n ( 31
)
W =    -cos α AC cos β AC -sin β AC -sin α AC cos β AC   
According to that point, we defined a forth type of frame, namely a flow frame x F , y F , z F consistent with twodimensional aerodynamic conventions. This orthonormal triad of unit vector is defined as follow (Fig. 6):

x F = -W
The lift forces directed by y F is normal to the wind and in the plane of the airfoil directed by x B which is also the direction of the lift moment :

y F = x B ∧ x F ∥x B ∧ x F ∥
The lift moment directed by z F complements the orthonormal triad of unit vectors:

z F = x F ∧ y F
Using these definitions, the tight coupling requires to link aerodynamic unknowns (section 2.3) with structural one (section 2.1). We defined, under the assumption that wind direction and velocity are constant in time or quasi-steady:

• α the local AoA, α ∈] -90°, 90°[: α = arcsin y B .y F ( 32
)
The time derivatives α, α, ḣ , ḧ :

α = Ω a .z F = C T C a b Ω B .z F (33) α = Ωa .z F + Ωa . ż F (34) ḣ = -V a .y F = -C T C a b V B .y F (35) ḧ = -V a .y F -V a . ẏ F (36)
Then, unsteady aerodynamic loads, considered constant within each element, are injected in distributed beam loads along with gravity forces as:

f a = L y F -µg z a ( 37
)
m a = M z F -µg x m2 z a ( 38
)
using Eqs ( 32)-( 36) to substitute aerodynamic unknowns with structural ones. The aerodynamic model adds N S equations for each member, the coupled aeroelastic system contains (18 + N S ) N + 12 equations and the same number of unknowns, providing that structural unknowns are completed with N ×N S induced-flow states λ n i . Hereinafter the system of equations will be written as:

(X , Ẋ ) = 0 (39)
This system permits a non-iterative computation of the unsteady aerodynamics loads, thanks to the Peters theory and the tight coupling between both models. This is in contrast with classical iterative "k method" and "p-k method" used to determined aeroelastic critical speeds.

In case of low reduced frequency (f r ≪ 1), a simplified aerodynamic model with λ 0 = 0, called "quasi-steady with added mass", can be used. In that case, the system has the same size as the initial purely structural formulation.

Transient motion and critical speed computation methodology

The formulation described above permits different applications:

• study nonlinear transient dynamic behaviour using a time marching scheme;

• study nonlinear steady-state behaviour by neglecting all the time derivatives;

• realise an eigenvalue analysis of small motion about a nonlinear steady-state by linearising about it.

A python pre/postprocessor is developed, including the homogenisation algorithms detailed in section 2.2, in order to exploit the aeroelastic solver.

The resulting open source computation code, thereafter called "GEBTAero"1 , has different capacities summarized on Fig. 7.

Quasi-steady algorithms

The simplest application is the determination of the static shape of the wing at a specific wind velocity using a steady state solution of Eq. ( 39):

(X , 0) = 0 [START_REF] Tang | Experimental aeroelastic models design and wind tunnel testing for correlation with new theory[END_REF] It is also possible with a steady-state solution to determine the torsional divergence speed by monitoring the divergence of the wing root internal forces and moments F B and M B for a very small AoA.

Eigenvalue analysis 3.2.1. Solving the eigenvalue problem

In the frequency domain, one of the main advantage of the tight coupling is the ability to directly compute aeroelastic modes. In that case, the solution vector is written:

X (t ) = X + X (t ) (41)
with X the steady-state solution and X (t ) ≪ X . As a results in Eq. ( 39) and performing a Taylor expansion, the system becomes, keeping first order terms:

X + X (t ), X = (X , 0) + ∂ F ∂ X X + ∂ F ∂ X X = 0 using Eq (40): ∂ F ∂ X X + ∂ F ∂ X X = 0 (42)
For an eigenvalue analysis, we assume that:

X (t ) = X 0 e νt ( 43
)
It gives us the following Generalised Eigenvalue Problem (GEP):

∂ F ∂ X + ν ∂ F ∂ X X 0 = 0 (44)
By analogy with classical dynamic resolutions, ∂ F ∂ X refers to a stiffness matrix and ∂ F ∂ X to a mass matrix . They must not be confused with stiffness matrix K and mass matrix I defined in section 2.1.

Because of the use of Lagrange multiplier inherent to the geometrically exact beam theory and the addition of induced-flow ODE, this GEP is real and non-symmetric. Neither stiffness nor mass matrix is symmetric positive semi-definite, thus our GEP requires a direct transformation into a Standard Eigenvalue (SEP) :

-1

X 0 = - 1 ν X 0 ( 45
)
Direct inversion of is avoided using the resolution of a linear system:

V = Z ( 46
)
with V = -1 X 0 and Z = X 0 .

Computing aeroelastic critical speeds

Once the SEP is properly defined, more work is needed to produce aeroelastic critical speeds. For the sake of efficiency, in particular for optimisation applications, the number of computed modes should be as small as possible. In this regard, the algorithm focuses on the mode showing the smallest damping for a given flow speed, either zero frequency (for the divergence speed) or non zero frequency (for the flutter speed).

To allow the proper tracking of these few modes, the SEP is resolved using the Implicitly Restarted Arnoldi Method (IRAM) implemented into open source ARPACK program (ARnoldi PACKage) [START_REF] Lehoucq | ARPACK Users' Guide: Solution of Large Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods[END_REF]. Matrices are stored in coordinate list sparse format (COO). The sparse direct linear solver required by Eq. ( 46) is MUMPS (MUltifrontal Massively Parallel sparse direct Solver) [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF]. A significant proportion of computation time comes from the resolution of Eq. ( 46), then special attention is devoted to optimise the configuration of MUMPS. Before the resolution of the SEP, rows and columns are reordered using software package Scotch [START_REF] Pellegrini | Scotch: A software package for static mapping by dual recursive bipartitioning of process and architecture graphs[END_REF] and are scaled.

The aeroelastic behaviour of the airframe is described by the lowest frequency structural modes modified by the fluid forces. However, Induced-flow ODEs [START_REF] Cartraud | Computational homogenization of periodic beam-like structures[END_REF] add N × N S eigenvalues to the initial structural problem. Most of them are close to those associated with the uncoupled ODE system derived from Eq. ( 26):

A λ + U b λ = 0 ( 47
)
If ν i are the eigenvalues of matrix b U A, it corresponds to the -1/ν i according to Eq. [START_REF] Brown | Integrated strain actuation in aircraft with highly flexible composite wings[END_REF]. For example, if N s = 6, A has 4 real eigenvalues : ν 1 ≃ 0.3865, ν 2 ≃ 0.5011 , ν 3 ≃ 2.8461, ν 4 ≃ 16.5383 and 2 complex conjugate eigenvalues ν 5,6 ≃ 6.3987.10 -2 ± 0.1546i . It results in negative real part eigenvalues -1/ν i . All these modes are stable but they are proportional to the flow velocity U , then every pair of complex conjugate mode of A generates up to 2N eigenvalues able to take any frequency depending on U . Because ARPACK is designed to quickly compute only a few eigenmodes and in order to avoid hiding unstable modes with induced flow generated stable modes, the parameter "which" of the modal solver is set to "SR" corresponding the the smallest algebraic real part eigenvalues ν. After inverting the system with Eq. ( 45), it corresponds to eigenvalues -1/ν with positive real part ordered by ascending real part followed by eigenvalues with negative real part ordered by decreasing real part. This setting allow us to quickly find the smallest flow velocity associated with a real unstable mode, namely the divergence speed, and the smallest flow velocity associated with a complex unstable mode, namely the flutter speed. Both speed could be computed on undeformed wing (X = 0) or on deformed wing.

Displaying aeroelastic modes

For visualisation purposes, eigenvalue analysis is also used to plot a few number of modes versus flow velocity. Because each flow velocity value corresponds to a specific eigenvalue problem, it is therefore necessary to establish links between eigenvalues and eigenvectors corresponding to the evolution of the same aeroelastic mode. For that purpose, modes of interest are identified at zero flow velocity to the first n i n t e r e s t structural modes. It thus permits to filter modes produced by induced flow speed computation. Then, for each flow velocity increment ∆U , a sufficient number of eigenmodes (n c a l c ul a t e d > n i n t e r e s t ) are calculated with ARPACK parameter "LM" (large magnitude) corresponding after the problem inversion to the smaller magnitude eigenvalues. Then, they are matched to previous modes using eigenvectors correlation. ∆U and n c a l c ul a t e d are dynamically adjusted in the algorithm to allow a proper correlation, in particular to avoid mode crossing phenomenon As illustrated in section 4, eigenvalue analysis permits to obtain:

• the n first aeroelastic modes (frequency and reduced damping) between U mi n = 0 m/s and U ma x ;

• the undeformed wing flutter and divergence speed (X = 0);

• the deformed wing flutter and divergence speed(X ̸ = 0).

Transient dynamic algorithms

Using the Newmark time scheme as detailled in [START_REF] Yu | GEBT: A general-purpose nonlinear analysis tool for composite beams[END_REF] with γ = 1/2 and β = 1/4, a transient dynamic simulation of the airframe can be done. A file with output data for each time step is generated allowing a visualisation of the wing motion in ParaView. This kind of simulation is illustrated in section 4.2 but temporal simulations are not well suited for structural optimisation which is the main objective of this work.

It should be noted that the temporal study of the flexible wing motion is a key issue of the non linear aeroelasticity. Indeed, it allow a proper comparison between experimental measurements and numerical models and permits for example to simulate accident like the one of Helios [START_REF] Su | Dynamic response of highly flexible flying wings[END_REF] or to study very flexible aircraft control laws like it is done in [START_REF] Virgilio Pereira | Model predictive control architectures for maneuver load alleviation in very flexible aircraft[END_REF]. The numerical model described above still presents some limitations for this kind of applications because the flight mechanics is not yet taken into account, neither the stall phenomenon. Hence it cannot, for instance, simulate the limit cycle oscillations observed experimentally in [START_REF] Kirsch | Numerical and experimental study of aeroelastic tailoring effect using flexible composite laminates for haps application[END_REF][START_REF] Tang | Experimental aeroelastic models design and wind tunnel testing for correlation with new theory[END_REF]. However, once these requirements will be met, the efficiency of the model makes it relevant for such temporal simulations.

Validation test cases

Homogenisation algorithm

First of all, the homogenisation routine is tested against analytical results for simple shape constant cross section with isotropic material. Three shapes are tested : a thin plate (length L, height h), a rectangular box (length L, height h, thickness e) and a cylinder (radius R). The first two shapes are automatically generated with CalculiX GraphiX whereas the cylinder mesh is done manually. The flexibility matrix generated is compared to analytical values [START_REF] Young | Roark's formulas for stress and strain[END_REF] in Eq. 17. Results are presented in Tab. 1 and show very good agreement even for the torsional flexibility S 44 estimated within a precision of 0.5%.

Critical speed computation

The critical speed computation is validated against common aeroelastic test cases, namely the Goland wing [START_REF] Goland | The flutter of a uniform cantilever wing[END_REF] and the Patil Wing [START_REF] Patil | Nonlinear aeroelastic analysis, flight dynamics, and control of a complete aircraft[END_REF]. The first is universally used among literature and the second is more appropriate to our program because of its high aspect ratio. Unfortunately, both of them concern isotropic wing since there is no common anisotropic test case among the literature. Characteristics of both wing are detailed in Tab. 2. First of all, a convergence of the undeformed Goland wing flutter speed relatively to N and N s parameters are realised (Fig. 8). It appears that a reasonable value for N start from 6, parameter N S shows good results from 4. Because of the factorial entering into definition of A (Eq. ( 27)) and b (Eq. ( 29)), N s is limited to 8 in order to avoid round-off errors and ill-conditioned matrices.

Then, the flutter speed and frequency are computed at sea level (ρ = 1.225 k g .m -3 ) and at 20000 f t (ρ = 0.6526 k g .m -3 ) using N = 10 elements and N S = 6 induced-flow states and compared to literature. The corresponding Mach number at 20 000 ft is 0.56, which is quite large for an incompressible model, but it is still relevant for the sake of comparison, provided that other models use this hypothesis. Results are presented in Tab. 3.

Results show a good agreement with other strip theory reduced order models whereas UVLM based model (SHARP) predict a higher critical speed because of three-dimensional effects. A more appropriate test case is the Patil high aspect ratio wing whose characteristics are detailed in Tab. 2. Besides undeformed wing flutter speed assessment, flutter speed is also calculated about the wing deformed by its own weight (static wing tip deflection = 2.93 m). Both simulations are made using N = 10 and N S = 6 and demonstrate also a good agreement (Tab. 4). The Mach number is of the order of 0.1, in accordance with incompressible hypothesis. In order to illustrate the impact of the large displacement and rotation of the Patil wing on the aeroelastic behaviour, the first four modes are plotted for the undeformed configuration (Fig. 9a) and the deformed one (Fig. 9b). In our case, the unstable 137.2 70.7 --NATASHA [START_REF] Patil | Nonlinear aeroelastic analysis, flight dynamics, and control of a complete aircraft[END_REF] 135.6 70.2 --UM/NAST [START_REF] Brown | Integrated strain actuation in aircraft with highly flexible composite wings[END_REF] 136.2 70.2 174.9 68.1 SHARP [START_REF] Murua | Assessment of wake-tail interference effects on the dynamics of flexible aircraft[END_REF] 165 69 --Aeroflex [START_REF] Ribeiro | Aeroflex: a toolbox for studying the flight dynamics of highly flexible airplanes[END_REF] 137.0 70.8 177.0 69.2 The longitudinal distance between two snapshots corresponds to the distance covered during the time step and the colorscale represents the wing rotation magnitude (in radian). This display intends to illustrate the large displacement and rotation of the airfoils. It should be recalled that the flight mechanics is not taken into account. For further study of Patil wing LCO, the reader is referred to [START_REF] Ahmad | Limit cycle oscillation prediction based on finite element-modal approach[END_REF].

Furthermore, computation time is assessed on the Patil wing flutter speed with a precision of 10 -1 m/s (Fig. 11). The results are promising with a computation time in the order of O (N ), paving the way for complex optimised composite airframe simulations. Deformed wing computation seems to be faster, which can be explained by the lower value of critical speed to reach by the algorithm. The impact of parameter N S is complex The divergence speed computation time is very similar with for instance 4.93 s (undeformed) and 6.25 s (deformed) for N s = 6 and N = 100.

Aeroelastic tailoring effect assessment

Thanks to the computation speed and automation routines of the program, since no common experimental or numerical test cases have been found in the literature, the aeroelastic tailoring effect is numerically assessed on simple composite laminates. In addition to a sandwich configuration evaluated in [START_REF] Kirsch | Assessment of aeroelastic tailoring effect on high-aspect-ratio composite wing flutter speed using an open source reduced order model solver[END_REF] with the variation of the upper ply orientation, the present paper evaluate an anisotropic flexible wing test case derived from Patil wing and a two-ply composite laminates (carbon/epoxy T700GC/M21) by sweeping both ply orientations (lower ply θ 1 and upper ply θ 2 ).

Representative high aspect ratio anisotropic composite wing

Both wing test cases presented above are isotropic. To illustrate the anisotropic capability of the model, coupled with geometrical non linearities, the Patil wing is modified by introducing a bending/twisting coupling. In order to produce a realistic coefficient, a composite box with a specific layup is defined with the same bending (S 55 ) and twisting (S 44 ) flexibility as Patil wing compatible with the wing section dimensions. The composite box characteristics are detailed in Tab. 5 using carbon/epoxy T700GC/M21. The ply order starts from the center of the box. The composite box is centred on the original wing elastic axis. Using beam homogenisation tool of GEBTAero, the resulting flexibility matrix is:

S =       
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.01.10 -4 5.88.10 -6 0 0 0 0 5.88.10 -6 4.97.10 -5 0 0 0 0 0 0 0

      
This configuration permits a very good agreement with bending and twisting flexibility of the original Patil wing test case (Tab. 2). Thus, bending/twisting coefficient S 45 is added to the Patil wing test cases, all else being equal. The deformed (by its own weight) and undeformed critical speeds are presented in Tab. 6 alongside isotropic reference. The first four aeroelastic modes of both configurations are plotted on Fig. 12.

The divergence speed of this anisotropic test case is about twice the speed of the isotropic wing for the undeformed configuration and almost three time the speed for the deformed configuration. Concerning the flutter, the behaviour of the undeformed configuration is quite similar as the isotropic test case. However, for the deformed configuration, the impact on the unstable mode (green) is obvious. In fact, the bending twisting coupling tends to cancel the static deflection of the wing, as it can be seen in Fig. 13. Thus, this example confirms that the aeroelastic tailoring is a satisfactory way to compensate the aeroelastic performance decrease due to the static deflection.

Two-ply flexible laminate

The characteristics of the two-ply flexible laminate are detailed in Tab. 7. To ensure a good critical speed convergence, the number of beam elements N is set to 20 and the number of induced-flow states N s is set to 6, the gravity effect is taken into account and the AoA is set to 0°. A critical speed computation is performed every 3°for both ply orientations (3600 simulation points computed in 2h 45min on a single core) . At each point, the flexibility matrix is determined using the homogenisation algorithm, particularly the torsional flexibility S 44 , the bending flexibility S 55 and the bending/twisting coupling S 45 from Eq. 18. Then, using the method detailed in section 3.2, the divergence speed, the flutter speed and the flutter frequency is computed. Simulations are made at sea level (ρ = 1.225 kg.m -3 ), knowing that the Reynolds number is rather small for the use of an inviscid potential flow theory (from 10 4 to 10 5 for a speed from 5 to 50 m/s).

Results are plotted on Fig. 14 and Fig. 15. Concerning the bending flexibility (Fig 14a) and the torsional flexibility (Fig. 14b), results are consistent with classical laminate theory (CLT). In fact, the bending flexibility is minimal (stiffer configuration in bending) if all the fibres are oriented along the beam direction (0°/0°) and maximal if all the fibres are oriented transversely (90°/90°). The torsional flexibility is minimal (stiffer configuration in torsion) for a balanced layup oriented at ±45°(-45°/45°or 45°/-45°) and relatively low if fibres are all oriented at -45°or 45°(unbalanced layups). The bending/twisting coupling shape (Fig. 15a) is more complex. This parameter can be positive or negative and, according to the CLT, balanced layups are devoid of coupling. The maximum in magnitude is obtained for unbalanced layups 30°/30°and -30°/-30°. All these flexibility coefficients are symmetrical about the mirrored layup axis (θ 1 = θ 2 ) because the homogenisation problem is symmetrical by switching both ply orientations. The balanced layup axis (θ 1 = -θ 2 ), corresponding to the change of sign of both ply orientations, is a symmetry axis for bending and twisting but an anti-symmetry axis for the coupling coefficient. Moreover, all coefficients show a smooth and regular aspect consistent with the continuous nature of the homogenisation problem.

Because the cross section parameters are unchanged by switching both ply orientations, the mirrored layup axis is also a symmetry axis for the critical speed (minimal speed between divergence and flutter; Fig. 15b) and the associated critical frequency (Fig. 15c). A strong correlation appears between the critical speed and the bending/twisting structural coupling, namely the aeroelastic tailoring effect. Indeed, according to the frame of Fig. 2, a positive lift will bend the wing negatively. In case of a positive bending/twisting coupling it will twist the wing negatively corresponding to a decrease of the AoA which has the effect of stabilising the wing deflection amplitude. High critical speed are mainly found for layups stiff in torsion but the correlation is less obvious and are always found for layups flexible in flexion, which seems counter-intuitive. Finally, optimal configurations are found for balanced layup (without coupling), but it is important to note the strong impact of the coupling coefficient change of sign with a halving of the critical speed within a tiny range of 5°ply orientation. Then, to overcome this difficulty, the design range should preferably be chosen in the positive coupling area. The critical frequency (Fig. 15c) permits to know whether the instability is a divergence (0 H z ) or a flutter. Layups with negative structural coupling are prone to torsional divergence as in the case of the geometrical coupling generated by a forward swept wing [START_REF] Weisshaar | Aeroelastic tailoring of forward swept composite wings[END_REF].

Conclusion

Design challenges induced by HAPS in terms of aeroelastic performances show the need for an accurate reduced order model able to simulate nonlinear behaviour of an anisotropic high aspect ratio wing. Moreover, innovative solutions like aeroelastic tailoring require a fast computation of critical speeds compatible with an optimisation process. The present work presents the computation code GEBTAero, a solution based on the geometrically exact beam theory coupled with a two-dimensional unsteady finite state aerodynamic model implemented into an open source solver. The homogenisation step required by the beam model is either analytical or performed by a 3D FEM calculation with periodic boundaries using CalculiX solver. Accuracy and swiftness of critical speed computation on both undeformed and deformed wing has been evaluated using common aeroelastic test cases with different problem sizes, proving the ability to quickly compute complex design, and thus to be used within a multidisciplinary optimisation platform. Configurations using aeroelastic tailoring, which are the core target of this solver, are then evaluated numerically on a representative high aspect ratio anisotropic composite wing and a simple 2-ply composite laminates with both variable ply orientations. It illustrates the strong correlation between the structural bending/twisting coupling of an unbalanced composite laminates and its critical aeroelastic speed. It also shows the high sensitivity of ply orientation on the aeroelastic behaviour. Furthermore, the lack of experimental data for this type of anisotropic specimen show the need for a wind tunnel campaign dedicated to very flexible composite laminates using aeroelastic tailoring.

Figure 2 :

 2 Figure 2: Structural frames definition

S

  11 0 0 S 14 S 15 S 16 0 14 0 0 S 44 S 45 S 46 S 15 0 0 S 45 S 55 S 56 S 16 0 0 S 46 S 56 S 66 Composite box section with two plies vertically and three plies horizontally (b) Full cross section NACA0015 airfoil mesh

Figure 3 :

 3 Figure 3: Examples of meshes used for the homogenisation

Figure 4 :

 4 Figure 4: Warping of a NACA0015 airfoil section (amplified); CalculiX GraphiX

Figure 5 :

 5 Figure 5: Airfoil parameters

Figure 6 :

 6 Figure 6: Flow frame definition

Figure 7 :

 7 Figure 7: GEBTAero computation features

Figure 8 :

 8 Figure 8: Convergence tests for N and N s parameters for the Goland test case

  speed (m/s ) frequency (r a d /s ) speed (m/s ) frequency (r a d /s ) present (N = 10; N S = 6

Figure 9 :

 9 Figure 9: Patil wing aeroelastic modes

Figure 10 :Figure 11 :

 1011 Figure 10: Real time flutter instability of the Patil wing; altitude = 20 k m ; velocity = 38 m /s

Figure 12 :Figure 13 :

 1213 Figure 12: Anisotropic Patil wing aeroelastic modes plot

Figure 14 :

 14 Figure 14: Composite laminates numerical simulation depending on lower and upper ply orientation : bending and torsional flexibility (symmetry axis in white dashed line)

  (a) bending/twisting coupling (b) critical speed (c) Critical frequency

Figure 15 :

 15 Figure 15: Composite laminates numerical simulation at sea level (ρ = 1.225 kg.m -3 ), depending on lower and upper ply orientation : bending twisting coupling, critical speed and frequency (symmetry axis in white dashed line; anti-symmetry axis in white dotted line)

  12 S 13 S 14 S 15 S 16 S 12 S 22 S 23 S 24 S 25 S 26 S 13 S 23 S 33 S 34 S 35 S 36 S 14 S 24 S 34 S 44 S 45 S 46 S 15 S 25 S 35 S 45 S 55 S 56 S 16 S 26 S 36 S 46 S 56 S 66

Table 1 :

 1 Validation of the 3D FEM with periodic boundaries homogenisation algorithm against analytical values[START_REF] Young | Roark's formulas for stress and strain[END_REF]

			Thin plate	Rectangular box	Cylinder
	Dimensions	m L=5 × 10 -2 ; h=10 -3 L=0.1; h=5 × 10 -2 ; e=5 × 10 -4	R=1
	Area	m 2	5 × 10 -5	× 10 -4	3.1416
	Area moment of inertia I G 2 m 4 Area moment of inertia I G 3 m 4 Torsion constant J m 4	4.1666 × 10 -12 1.0416 × 10 -8 1.6456 × 10 -11	7.1062 × 10 -8 2.0461 × 10 -7 1.6281 × 10 -7	0.7854 0.7854 1.5708
	Number of elements		100	200	2792
	Computation time	s	0.54	1.33	7.37
	S 11 relative error	%	3.37 × 10 -6	8.35 × 10 -6	4.14 × 10 -2
	S 44 relative error	%	-5.03 × 10 -3	-0.44	0.13
	S 55 relative error	%	3.38 × 10 -6	-4.08 × 10 -6	8.25 × 10 -2
	S 66 relative error	%	2.82 × 10 -5	-5.22 × 10 -6	8.37 × 10 -2

Table 2 :

 2 Characteristics of Goland and Patil wingsGoland[START_REF] Goland | The flutter of a uniform cantilever wing[END_REF] Patil[START_REF] Patil | Nonlinear aeroelastic analysis, flight dynamics, and control of a complete aircraft[END_REF]

	Semi-span	L	m	6.096	16
	Chord	2b	m	1.8288	1
	Mass per unit length	µ	kg/m	35.71	0.75
	Elastic axis EA (from leading edge)		% chord	33	50
	Center of gravity CG leading edge)		% chord	43	50
	Distance between CG and EA	v	m	0.18288	0
	Bending stiffness (spanwise)	E I G 2	N.m 2	9.77 × 10 6	2.10 4
	Bending stiffness (chordwise)	E I G 3	N.m 2	/	4.10 6
	Torsional stiffness	G J	N.m 2	0.99 × 10 6	1.10 4
	Mass moment of inertia around e.a.	i 11	kg.m	8.64	0.10

Table 3 :

 3 Goland wing flutter speed and frequency

Table 4 :

 4 Patil wing flutter speed and frequency

	program	Undeformed wing speed (m/s ) frequency (r a d /s ) speed (m/s ) frequency (r a d /s ) Deformed wing
	present (N = 10; N S = 6)	32.2	22.6	23.3	11.9
	NATASHA[43]	32.2	22.6	-	-
	UM/NAST[46]	32.2	22.6	23.2	10.3
	Aeroflex[16]	32.6	22.3	23.4	12.2

Table 5 :

 5 Anisotropic Patil wing composite box characteristics

	Parameter	Unit	Value
	width	mm	200
	height	mm	67
	ply thickness	mm	0.125
	upper layup	-	[0, 45, 45, 0]
	lower layup	-	[0, -45, 45, 0]
	left layup	-	[0, 45, -45, 0]
	right layup	-	[0, -45, 45, 0]
	E l	G p a	148
	E t	G p a	10
	G l t	G p a	4.6
	ν l t	-	0.3
	ρ	k g .m -3	1.6

Table 6 :

 6 Anisotropic Patil wing critical speeds and frequencies computed by GEBTAero

			flutter	divergence
	configuration	speed frequency	speed
		(m/s)	(rad/s)	(m/s)
	isotropic undeformed (without weight)	32.2	22.6	38.0
	isotropic deformed (with weight)	23.4	11.9	38.0
	anisotropic undeformed (without weight)	32.3	22.6	78.3
	anisotropic deformed (with weight)	31.6	22.1	97.0

Table 7 :

 7 2-ply composite laminates characteristics

	Parameter	Unit	Value
	Chord	mm	30
	Length	mm	300
	total thickness	mm	0.5
	ply thickness	mm	0.25
	E l	G p a	148
	E t	G p a	10
	G l t	G p a	4.6
	ν l t	-	0.3
	ρ	k g .m -3	1.6

https://framagit.org/BertrandK/GEBTAero


THIN FILMS OF LITHIUM GLASSES

A. LEVASSEUR Laboratoire de Chimie du Solide du CNRS Université de Bordeaux I 351, cours de la Libération -33405 TALENCE CEDEX (France).

Before approaching the lithium glasses in form of thin films the two main methods to obtain these materials will be des The pressure of €he atmosphere in which deposition is car ried out optimally depends on the applied technology.

Evaporation

The number of atoms or molecules evaporating from a liquid or solid surface depends strongly on the temperature. As is generally known, the equilibrium vapour pressure would be obtained in a thermodynamically closed system. However, in practical evaporation, no equilibrium is obtained because the environment of the vapour source acts as a vapeur sink. The evaporant atoms condense on all parts that are at lower temperature than the vapour source.

Systematic investigations of evaporation rates in vacuum have been performed mainly by Hertz, Knudsen and Langmuir. In these experiments it was found that a liquid has a specific ability to eva porate and cannot exceed a certain maximum evaporation rate at , a given temperature even if the heat supply is unlimited. The theoretical maximum evaporation rates are obtained only if the number of vapeur molecules leaving the surface corresponds to that required to exert the equilibrium pressure p e on the same surface, with none re turning to it.

From these experiments and considerations, the following equation for the molecular evaporation rate could be formulated

• • • • • • • • • • • • • • • • • • • •• • . . .. . . . .. . . . . •• • • •• • • • Nuclei grow

Coalescence

Stages of film growth.

Shape changes dunng coalescence.

FIG. 4

Coalescence Channels

Holes

Continuous film

Schematic representation of the stages of metal film growth and approximate relationship between dislocati. on density and film thickness according to Pashley.

Most films deposited even at room temperature are in a non equilibrium state and highly imperfect containing vacancies, disloca tions, stacking faul ts a• na grain boundaries, unless there is some me chanism for achieving equilibrium. The method of approaching equili brium is by movement of atoms in and on the surface layers. The most important parameter controlling the mobility of atoms in a solid film is diffusion. Therefore if the condensation process occurs closer to the melting point of the film material, a better o�dered solid film is formed. This can be achieved, for example, by increasing the substrate temperature.

ordering procesies, grain size.

In �ddition to influencing surface mobility and the substrate temperature will also affect the c) Reactive_evaporation

For the deposition of a stoichiometric oxide film by reactive evaporation, a relatively high o 2 partial pressure and a slow metal-atom condensation rate are required, so that completely oxidized metal-oxide films can be formed. The partial pressure of the reactive gas cornponent is usually few 10-4 mbar. The significant technology of reactive evaporation is applied in all cases where direct evaporation of a chemical compound is not possible because of thermal dissociation or too low a vapour pressure. In practice, oxide films are usually produced using sub-oxides or metallic starting ma terials. However, basically it is also possible to produce sulphides and nitrides or other compounds in this manner.

General considerations

Without going into <letail, we will attempt to discuss thr� parameLers ùnd system components of significance for cathode sput tering, and thus give a goneral idea of the fundamental correlations of cathode sputtering and also of the opportunities offered by this proccss.

Fig. 5 shows a schematic drawing of a greatly simplified sct-up for cathode sputtering. The process takcs place in a vacuum chamber, which has been evacuated as well as possible before coating.

In order to prevent contamlnation of the films ta be produced by incorporation of less dcfined should be 10-6 mbar or lower. rcsidual gas, the starting pressuru

The working pressure is then achievu<l with the working gas. The sputtering process itself takes place in a gas discharge which is ignited at a pressure between 10-J and 10-2 mbar, dcpending on the special variant of the method. A requirem0nt of the vacuum pumps is therefore an ultimate vacuum as high as pos sible and also a high and constant pumping speed in the mbar range. All of these processes take pl�ce, both at the target and at the substrate and determine the properties of the growing films.

About 95 % of the energy of the bombarding ions on the target is lost as heat in the solid material and only 5 % is passed onto the secon dary particles.

With this technique the deposition rate is generally low A similar phenomenon is observed for the activation ener gies (Fig. 10) .

The ionic conductivity varies exponentially with Li 2 o content and the activation energy varies linearly with Li2o content (Fig. 9, 10). Similar behavior had already been observed for glasses The starting matarials and the obtained thin films are very hygroscopie.

The deposition of such glasses requires a small evaporation apparatus, which has been designed and built up in order to be put easily into a glove box. This glove box is equipped with water vapor and oxygen purification devices.

The deposition rate controlled by a piezoelectrical balance is about 140 A s -1 , five times higher than that determined with oxide glasses.

The obtained films are checked by SEM.

Electrical measurements

The high ionic conductivity of the glasses does not allow using a sandwich electrode structure for electrical measurements.

The arrangement of the electrodes used for ionic and electronic conductivity measurements of -�i ---- -during the first temperature increase from 25 ta 90 °C the measured conductivity is similar ta that of the starting glasses (Fig. 14).

This conductivity is due essentially to the motion of the ions in the bulk of the thin film. Thus despite the uncertainty of the exact iodine content the films seem to have the same composition and the same behavior as the starting glasses.

-the ten times higher conductivi ty measured when the films are heated for several hours at 90 °C seems to be independent of the composition of the thin films (Fig. 14). Such a type of annea ling never produces any conduc tivity enhancement on similar composition bulky glasses.

  EVAPORATION AND CONDENSATION IN HIGH VACUUM Thin films of materials of such different properties as metals, halides, oxides and sulphides can be obtained in the crystal line or amorphous state by condensation of the vapour. The mechanism involved in forming the film may be a pure physical e.g. simple con densation or may also involve chemical reactions. Evaporated films were probably first made by Faraday in 1857 when he exploded metal wires in an inert atmosphere. The deposi tion of metal films in a vacuum by resistance heating of platinum wires was performed by Nahrwold in 1887 and only a year later this technique was used by Kundt to produce films for measuring the re fractive indices of metals. In the following period, evaporated thin films had only academic interest, although the vacuum evaporation of metal films by Pohl and Pringsheim in 1912 was performed under better developed technological conditions. The techniques of forming thin films by condensation at very low pressures have been developed in parallel with industrial development of techniques of producing high vacua in large volumes, using pumps of very high speed and suitable cieanliness, Films obtained by condensation alone are usually made at pressures between 10-6 and 10-B mbar and down to the lowest values that modern vacuum technology can attain. This is combined with fast deposition rates to produce deposits in which only few foreign gas atoms are incorporated. Of great technical importance was the discovery and deve lopment of the reactive evaporation process by Auwarter 1952 and Brinsmaid 1953. Such processes which in volve a chemical reaction between evaporated constituents and the gas atmosphere are carried out mainly in the 10-4 mbar region with slow deposition rate. The reaction is often activated either thermally or by uv radiation or alternatively by ionic oi electronic bombardment. The industrial development of these evaporation and conden sation techniques is due largely to the relative ease of obtaining either pure metal films or exact stoichiometric compound deposits of uniform thickness, because the laws which determine the phenomena at low pressure are better defined than close to or at atmospheric pres sure. At very low pressures, the mean free path of the vapour atoms or molecules exceeds the usual distance between evaporation source and the substrate. Thus the substrate receives a flux of vapour which travels in straight lines with no or only very few gas/vapour colli sions and the spatial distribution of the condensate obeys purely geometrical laws. If however, the residual pressure has such a value that the source to substrate distance•corresponds to the mean free path of the vapeur atorns, then only 37 % of the atoms travel without collision with the residual gas molecules. At residual pressures above• l0-3 mbar, the propagation of the evaporated atéms between source and the substrates is subject to an isotropie diffusion effect because of increasing frequency of intermolecular collisions.

  FIG.1

  FIG. 2Profile area A, of the thickness profile

  2. Film deposition by cathode sputtering Cathode sputtering or more exactly the sputtering of (usually solid) materials by bombardment with positive noble gas ions is the oldest vacuum process for producing thin films. Sputtering, the cause of the erosion of the cathode in glow discharge, often an undesirable effect, was discovered more than 130 years ago by Grove 1852 in England and Plücker in Germany during gas discharge experi ments. Saon afterwards in 1877, metal sputtering was applied in the production of rnirrors and later it was used for decorating va rious articles with noble metal films. Around 1930 it was used for applying electrically conducting films of gold onto the wax masters of the Edison phonographs. It then became of less importance for the next 30 years compared wi th the rapidly developing deposi tian of films by evaporation and condensation in high vacuum. 11<,'v/ùVCr, sincc about 1955 sputtcrinq has un<lergonc il rr,llê.t.i.fJ/3illîCC, Intensive studies of the phenomena occurring during the sputtcring process and hence better control of the process together with the tcchnical requirement for high quality films having very good adhesion and possessing specific properties have certainly con tributed ta the wide_ application.

  FIG. 5Simple diode type sputtering system.

  FIG. 7Glass thin film sample between two gold electrodes.Ionie conductivity was measured by a classical complex im pedance method using a SOLARTRON 1170 frequency pulse analyser in the frequency range of 0.1 mHz up to 10 4 Hz.

  FIG.11

  FIG. 12Arrangement of the electrodes for electrical measurements.

  FIG. 14 3 -l Variation of log a vs 10 .T for B 2 s 3 -Li 2 S-LiI glasses