Linearized Aeroelastic Computations in the Frequency Domain Based on Computational Fluid Dynamics
David Amsallem, Daniel Neumann, Youngsoo Choi, Charbel Farhat
LLinearized Aeroelastic Computations in theFrequency Domain Based on Computational FluidDynamics
David Amsallem , Daniel Neumann , Youngsoo Choi , Charbel Farhat Stanford University, Stanford, CA 94305-4035, USA
An iterative, CFD-based approach for aeroelastic computations in the frequencydomain is presented. The method relies on a linearized formulation of the aeroelas-tic problem and a fixed-point iteration approach and enables the computation of theeigenproperties of each of the wet aeroelastic eigenmodes. Numerical experiments onthe aeroelastic analysis and design optimization of two wing configurations illustratethe capability of the method for the fast and accurate aeroelastic analysis of aircraftconfigurations and its advantage over classical time-domain approaches.
I. Introduction
Computational aeroelasticity [1] has traditionally relied on linear methods such as the doubletlattice method [2] and the piston theory [3]. While these two approaches enable accurate aeroelasticpredictions in the subsonic and supersonic regimes, respectively, they fail to produce an accuratedescription of flutter in the transonic domain. Approaches based on computational fluid dynamics(CFD) and Arbitrary-Lagrangian-Eulerian (ALE) formulations, have subsequently been able toproduce very accurate flutter predictions in the transonic domain [4].When used in the time-domain, CFD-based aeroelastic predictions typically rely on a post Department of Aeronautics and Astronautics Department of Aeronautics and Astronautics Department of Aeronautics and Astronautics Department of Aeronautics and Astronautics, Department of Mechanical Engineering, Institute for Computationaland Mathematical Engineering, AIAA Fellow. a r X i v : . [ phy s i c s . f l u - dyn ] J un rocessing step [5] that extracts the aeroelastic characteristics of the wet structure, such as dampingratio and frequencies. This post processing step can however be inaccurate and miss aeroelasticcharacteristics associated with certain modes altogether. Furthermore, depending on the initialconditions of the time simulation, certain modes may not be excited.In the present work, an alternative approach based on linearized CFD and defined in the fre-quency domain is developed. It shares some similarities with the fixed-point iteration approachused in the doublet lattice method that can extract flutter frequencies. However, because it relieson CFD, it is accurate in the transonic regime and can also extract damping ratios for non-flutterpoints. As such, it provides the capability of leading to accurate aeroelastic predictions at anyoperating point. In turn, it will be demonstrated in this paper that such a procedure can be easilyintegrated in a design optimization procedure under aeroelastic constraints.This paper is organized as follows. Section II presents the framework of CFD-based linearizedaeroelasticity, focusing on the frequency domain formulation. The proposed approach for the fastcomputation of the system aeroelastic characteristics is developed in Section III. Application to theaeroelastic analysis and design of two wing systems is presented in Section IV. Finally, conclusionsare given in Section V. II. CFD-Based Linearized Aeroelasticity in the Frequency DomainA. Three-field Formulation
The following three-field ALE formulation is considered in this work to describe a fluid/structureinteraction problem: ∂ ( J W ) ∂t + J ∇ · (cid:18) F ( W ) − ∂x∂t W (cid:19) − J ∇ · R ( W ) = 0 ρ ∂ U ∂t − div ( σ ( (cid:15) ( U ))) − B = 0 (cid:101) ρ ∂ x∂t − div (cid:16) (cid:101) E : (cid:101) (cid:15) ( x )) (cid:17) = 0 . (1)The first equation in (1) is the ALE form of the Navier-Stokes equations, the second one, theelastodynamics equation for the structural subsystem and the third one provides the dynamics ofthe fluid grid, represented as a pseudo-structure. All tilde notations refer to that pseudo-structure.In all equations, t denotes time, x ( t ) denotes the vector of time-dependent fluid grid position, J
2s the determinant of the Jacobian of x with respect to the reference configuration and ∇ is thegradient with respect to x . W denotes the vector of conservative fluid variables and U the structural displacements. F and R are the vectors of convective and diffusive fluxes, respectively. ρ denotes material density, (cid:15) thestrain tensor and σ the stress tensor. div is the divergence operator. B is a vector of volume forcesacting on the structure. Finally E denotes a tensor of elasticities.Appropriate Dirichlet and Neumann boundary conditions as well as interface conditions applyto the set of coupled partial differential equations. The interface conditions between the fluid andstructural subsystems enforce the compatibility before the two velocity fields at the interface as wellas the equilibrium of tractions on the wet surface. Finally continuity equations are enforced at thewet interface between the structural and fluid mesh subsystems. The author is referred to [6] formore details on these interface conditions. B. Semi-discretization
The fluid subsystem is discretized in space by the finite volume method and the structural andfluid mesh subsystems discretized by the finite element method. The fluid is assumed to be inviscidin the remainder of this paper. The semi-discretization, after inclusion of the appropriate boundaryconditions, leads to the following set of coupled ordinary differential equations ˙ (cid:92) ( A ( x ) w ) + φ ( w , x , ˙ x ) = ¨ u + f int ( u , ˙ u ) = f ext ( u , w ) (cid:101) Kx = (cid:101) K c u . (2) x is the vector of fluid mesh nodes positions and A is a diagonal matrix containing the correspondingcell volumes. A dot denotes a time derivative. w is the vector of discretized conservative fluidvariables and u the vector of structural displacements. φ denotes the numerical flux function. M isthe mass matrix resulting from the finite element discretization and f int is the corresponding vectorof internal forces. f ext is the vector of external forces applied to the structure. Finally, the fluidmesh motion is modeled as being piece-wise static. (cid:101) K is the corresponding fictitious stiffness matrixand (cid:101) K c is a transfer matrix representing the effect of the structural motion u at the wet interface.3 . Linearization The set of nonlinear ODEs is then linearized around an equilibrium point ( w , x , u , ˙ u ) of thesystem. The resulting linearized system can then be written as a set of coupled linear ODEs A ˙ w + H w + ( E + C ) ˙ x + G x = ¨ u + D ˙ u + K u = P w (cid:101) Kx = (cid:101) K c u . (3)where the fluid mesh variable x has been eliminated through x = (cid:101) Su = (cid:101) K − (cid:101) K c u and the linearizedoperators are A = A ( x ) , H = ∂ φ ∂ w ( w , x , ˙ x ) , E = (cid:20) ∂ A ∂ x ( x ) w (cid:21) (cid:101) S , C = (cid:20) ∂ φ ∂ ˙ x ( w , x , ˙ x ) (cid:21) (cid:101) S , G = (cid:20) ∂ φ ∂ x ( w , x , ˙ x ) (cid:21) (cid:101) S , D = ∂ f int ∂ ˙ u ( u , ˙ u ) , K = ∂ f int ∂ u ( u , ˙ u ) − ∂ f ext ∂ u ( w , u ) , P = ∂ f ext ∂ w ( w , u ) . (4)The variables ( w , u , ˙ u ) denote now perturbations around the reference state ( w , u , ˙ u ) . For sim-plicity, the subscripts are dropped in the remainder of this work. D. Structural Modal Reduction
The dimensionality of the structural subsystem is subsequently reduced by modal truncation.For that purpose, a set of m eigenpairs { ω i , z i } mi =1 satisfying the following generalized eigenvalueproblem is first computed: Kz = ω Mz . (5)The eigenpairs are gathered in a diagonal matrix Ω = diag ( ω , · · · , ω m ) of eigenvalues and a modalbasis matrix Z = [ z , · · · , z m ] .The structural displacements and velocities are then approximated as a linear combination ofthe eigenmodes as u ( t ) ≈ Zu m ( t ) , v ( t ) ≈ Zv m ( t ) . (6)The structural equations are then reduced by Galerkin projection as I m ¨ u m + Ω u m = f m ( t ) (7)4here f m ( t ) is the reduced vector of linearized aerodynamic forces defined as f m ( t ) = Z T f ( t ) = Z T Pw ( t ) = P m w ( t ) . (8) E. Frequency Domain Formulation
The coupled equations resulting from the CFD-based ALE formulation can be written as I ¨ u m + Ω u m = P m wA ˙ w + Hw + ( E + C ) ˙ u m + Gu m = (9)The aerodynamic and aeroelastic linear operators H , E , C , G and P are computed for specificstructural configurations and free-stream boundary conditions such as free-stream Mach number M ∞ , angle of attack α as well as flight altitude h . Let s ∈ C and assume that w ( t ) = w e st , u m ( t ) = u m e st . (10)Note that s can have a complex part. Then their derivatives are ˙ w ( t ) = s w e st , ˙ u m ( t ) = s u m e st , (11)and the fluid equation becomes ( s A + H ) w + [ s ( E + C ) + G ] u m = . (12)Hence w = − ( s A + H ) − [ s ( E + C ) + G ] u m . (13)Plugging this expression for w into the structural equation leads to I ¨ u m + Ω u m = − P ( s A + H ) − [ s ( E + C ) + G ] u m (14)Defining the matrix of generalized aerodynamic forces A ( s ) = P ( s A + H ) − [ s ( E + C ) + G ] , thisleads to I ¨ u m + (cid:0) Ω + A ( s ) (cid:1) u m = . (15)Note that this equation is also valid for values of s with non-zero real part.5 II. Aeroelastic Eigen AnalysisA. Nonlinear Aeroelastic Eigenproblem
Consider now a motion u m ( t ) = u h e st (16)where s = ω ( i + γ ) = iω + ωγ . Then, ¨ u m = s u h e st = s u m and (15) becomes (cid:2) s I + Ω + A ( s ) (cid:3) u h = . (17)This equation will now be transformed into state-space form. For that purpose, consider thevelocity v m ( t ) = v h e st . Then v h = s u h . The operator A ( s ) can also be decomposed into its realand imaginary parts as A ( s ) = A R ( s ) + i A I ( s ) . (18)Let s be also decomposed into its real and imaginary parts s = s R + is I (note that s R = ωγ and s I = ω ). Then A ( s ) = A R ( s ) + i A I ( s )= A R ( s ) − s R A I ( s ) s I + s R A I ( s ) s I + is I A I ( s ) s I = (cid:18) A R ( s ) − s R A I ( s ) s I (cid:19) + s A I ( s ) s I . (19)Equation (17) can then be written as s I ( s u h ) + Ω u h + (cid:18) A R ( s ) − s R A I ( s ) s I (cid:19) u h + A I ( s ) s I ( s u h ) = (20)which, combined with v h = s u h leads to the system s Iv h + (cid:104) Ω + A R ( s ) − s R s I A I ( s ) (cid:105) u h + A I ( s ) s I v h = , v h = s u h (21)that can be written in matrix form as − (cid:18) Ω + A R ( s ) − s R s I A I ( s ) (cid:19) − A I ( s ) s I u h v h = s u h v h . (22)This leads to the following non-linear eigenvalue problem A ( s ) u h v h = s u h v h (23)6here A ( s ) is the real-valued matrix A ( s ) = − (cid:18) Ω + A R ( s ) − s R s I A I ( s ) (cid:19) − A I ( s ) s I . (24)Note that there are as many solutions to this eigenvalue problem as there are generalized degrees offreedom in the structural system. The eigenvalues come in pairs s j = s Rj ± is Ij , j = 1 , · · · , m .There are two main approaches that can be followed to solve the non-linear eigenvalue prob-lem (23).The first class of approaches is based on a fixed point iteration approach. In the context of the p-k method for solving the flutter problem, that approach was first introduced in [2]. At each iteration,a linear eigenvalue system is solved and the algorithm proceeds until convergence of the iterates toan eigenvalue. The advantage of that approach is that it relies only on a simple linear eigenvaluesolver and is therefore straightforward to implement. A drawback of this approach is the potentialslow convergence requiring many iterations [7]. [7] presents an analysis of the convergence proposesan extension based on Newton’s method to address the slow convergence issue. That extensionhowever requires the computational of the sensitivities of the aeroelastic operators with respect tothe reduced frequency.The second class of approaches is based on continuation to solve the non-linear eigenvalue prob-lem. Such an approach was proposed in [8] in the context of the p-k method. The main advantageof this method is its robustness that guarantees a computation of all the eigenvalues. However, thisapproach is expensive as it requires computing flutter predictions for an entire trajectory from highaltitude to the desired altitude. Further, computing the flutter path in a robust fashion may re-quire small increments that increase the computational cost. Finally, sensitivities of the aeroelasticoperator are required for a continuation approach.In the remainder of this paper, the fixed point iteration approach is chosen for its simplicity ofimplementation. It is described in Section III B.7 . Fixed-point Iteration Algorithm The algorithm proceeds, for each of the m structural eigenvalues, by iteratively computing onesolution of the eigenvalue problem (23). For that purpose, an initial guess is provided by the dryeigenvalues and associated eigenvectors. The operator A is then evaluated at the dry mode andall of its eigenvalues computed. The eigenvalue are sorted according to their imaginary parts andthe eigenvalue associated with the index of interest retained. The operator A is then evaluatedfor that eigenvalue. The procedure is repeated until convergence. A pseudo-code is provided inAlgorithm 1. If convergence is not achieved after a maximum number of allowed iterations, thecontinuation procedure mentioned above is triggered.The main difference with the p-k method developed in [2] is the fact that the real part of theeigenvalues is retained during the iterations, leading to a possible convergence of the procedure toan aeroelastic eigenvalue with non-zero real part. This means that aeroelastic eigenvalue can becomputed for non flutter configurations as well. This is due to the fact that the linearized CFD-based framework developed in Section II is equally valid at flutter and non-flutter configurations,unlike the p-k method that is based on the doublet lattice method, valid only for harmonic motions.Once the eigenvalues { s j = s Rj + is Ij } mj =1 of the nonlinear eigenproblem (23) are computed, theassociated aeroelastic frequencies and damping ratios can be readily computed as (cid:98) ω j = s Ij , η j = − s Rj (cid:113) ( s Rj ) + ( s Ij ) , j = 1 , · · · , m. (25) C. Complexity analysis
The main cost associated with the proposed procedure originates from the computation ofthe linear operator A ( s ) for values of s arising in the procedure. Computing A ( s ) requires thecomputation of the matrix of generalized forces A ( s ) = P ( s A + H ) − [ s ( E + C ) + G ] . In practice,computing the term ( s A + H ) − [ s ( E + C ) + G ] is the most expensive step. It relies on the solutionof m large-scale, sparse systems of equations ( s A + H ) x i = [ s ( E + C ) + G ] e i , i = 1 , · · · , m . Thesesolutions are typically computed by a Krylov-subspace iterative technique. In the present work, thepre-conditioned Generalized Minimal Residual technique (GMRES) [9] is used.Assuming, for simplicity, that the computation of each of the m aeroelastic eigenvalues requires8 lgorithm 1 Iterative computation of the eigenvalues of the nonlinear eigenproblem (23)
Require:
Dry structural frequencies { ω j } mj =1 , tolerance (cid:15) Ensure:
Eigenvalues { s j = s Rj + is Ij } mj =1 , for j = 1 , · · · , m do Let s j = iω j and k = 0 while k = 0 or s kj − s k − j > (cid:15) do k = k + 1 Compute all the eigenvalues { λ l = λ Rl ± iλ Il } mj =1 of A ( s k − j ) and order them by frequency such that λ I < · · · < λ Im Let s kj = λ Rj + iλ Ij if k ≥ k max then Call the continuation procedure end if end while end for N iter iterations, N iter m such solutions of linear systems are required. D. Computation of the Flutter Speed Index
The Flutter Speed Index (FSI) is a non-dimensional quantity that characterizes the flutterproperties of an aeroelastic system. It is computed for a fixed free-stream Mach number M ∞ anddensity ρ ∞ . The free-stream pressure p ∞ is increased until the onset of flutter. This is equivalent,for a perfect gas, to increasing the free-stream velocity V ∞ as V ∞ = (cid:114) γp ∞ ρ ∞ M ∞ (26)with γ = 1 . . The FSI is then defined for the critical value of flutter velocity V cr ∞ asFSI = V cr ∞ b S ω α √ µ (27)where b S is a reference length such as the semi-chord of the wing at the root, ω α is the frequencyof the first torsional dry mode of the structural system. The parameter µ is defined as µ = m S ρ ∞ (cid:98) V (28)9here m S is the structural mass of the system and (cid:98) V is the volume of a conical frustum with lowerand upper based diameters being the stream-wise chords at the root and the tip of the wing andwith height equal to the semi-span of the wing of the system.Flutter can be in practice detected from the outputs of Algorithm 1 by the presence of aneigenvalue s j of the operator A with positive real part. In practice, Algorithm 1 is run for increasingvalues of the free-stream pressure until flutter is detected. IV. ApplicationsA. Flutter Analysis of the Agard Wing 445.6
1. Agard Wing 445.6 model
The weakened Agard 445.6 model 3 research wing model is described in [10]. The dimensions ofthe wing are reported in Figure 1. A finite element model (FEM) consisting of shell elements,corresponding to degrees of freedom (dofs) is created and the eigenfrequencies associated withits first m = 4 structural modes compared in Table 1 to the ones measured experimentally in [10].Good agreements are observed between the two sets of eigenfrequencies.A finite volume mesh with more than , nodes is then created, leading to more than , dofs in the CFD model. Mode Frequency (Hz) Frequency (Hz) Relativenumber (experimental model in [10]) (present FEM model) discrepancy .
60 9 .
84 2 . .
17 39 .
54 3 . .
35 50 .
50 4 . .
54 96 .
95 5 . Table 1 AGARD Wing 445.6 retained modes for the structure
2. Flutter Analysis
A flutter analysis is first carried on for three different Mach numbers M ∞ ∈ { . , . , . } andthree altitudes h ∈ { , , } ft. For every combination of flight conditions, the proposed10 ig. 1 Geometry of the Agard wing iterative approach is applied to predict the aeroelastic eigen-characteristics of the system. Hence,for each case, four sets of real and imaginary parts for the eigenvalues are computed, as well as thecorresponding damping ratio, as defined in Eq. (25). The densities and pressures associated witheach altitude are given in Table 2. The results are graphically depicted in Figure 2(a) for aeroelasticpredictions at sea level, in Figure 2(b) for h = 10 , ft and in Figure 2(c) for h = 20 , ft.In addition, the results are compared to predictions in the time domain for which damping andfrequencies characteristics are extracted by fitting a canonical response to the transient lift. Thissimple approach only allows the extraction of a single eigenvalue, as opposed to the approachproposed in this paper that allows the extraction of all the eigenvalues. The comparisons withthe proposed iterative approach show that the time domain are able to accurately predict the11haracteristics of one of the aeroelastic eigenmodes, but that eigenmode is not always associatedwith the smallest damping ratio. As a result, the time domain approach ofter overestimates theminimal damping ratio of the system, as observed in Figure 3 where the damping ratios computed byeach approach are reported. Furthermore, the time-domain mispredicts flutter for M ∞ = 0 . and h = 10 k ft. On the other hand, the proposed iterative approach returns all the eigenvalues and hencethe one associated with smallest damping. The aeroelastic analysis for M ∞ = 0 . and h = 20 , ft takes 4.25 minutes on 8 processors for a tolerance (cid:15) = 10 − . Computing the steady-state takes1.1 minute and the iterative approach takes 3.15 minutes. Altitude (ft) Pressure (lbf.in ) Density (lbf.s .in − ) (sea level) . . × − ,
000 10 .
11 0 . × − ,
000 6 .
76 0 . × − Table 2 Free-stream pressures and densities for the altitudes considered
Next, FSIs are computed for a range of Mach numbers varying from M ∞ = . to M ∞ = 1 . using the proposed approach. For that purpose, the free-stream density is first fixed to its sea levelvalue. The results are reported in Figure 4. One can observe a characteristic flutter dip around M ∞ = 0 . and a range of flutter conditions for M ∞ ∈ [0 . , . . B. Analysis and Design Optimization Based on the ARW-2 Wing
1. Aeroelastic Research Wing ARW-2 Model
The second model is a more realistic wing. It is the Aerodynamic Research Wing Number 2(ARW-2) which was designed for flutter suppression and gust load alleviation [11, 12]. The wingis based on a supercritical airfoil, has an aspect ratio of 10.3 and a leading-edge sweepback angleof . ◦ . The geometry of the wing is depicted in Figure 5. A FE model of the wing consists of shell, beam and truss elements, resulting in degrees of freedom. The structural FE modelis depicted in Figure 6. m = 8 structural eigenmodes are then selected to model the structuralsubsystem. The corresponding eigenfrequencies are reported in Table 3.12 − − Real Part I m ag i n a r y p a r t M ∞ = 0 . Frequency Domain (sea level)Time Domain (sea level)Frequency Domain (10k ft)Time Domain (10k ft)Frequency Domain (20k ft)Time Domain (20 kft) (a) M ∞ = 0 . − − − −
10 0100200300400500600
Real Part I m ag i n a r y p a r t M ∞ = 0 . Frequency Domain (sea level)Time Domain (sea level)Frequency Domain (10k ft)Time Domain (10k ft)Frequency Domain (20k ft)Time Domain (20 kft) (b) M ∞ = 0 . − − − − − − Real Part I m ag i n a r y p a r t M ∞ = 1 . Frequency Domain (sea level)Time Domain (sea level)Frequency Domain (10k ft)Time Domain (10k ft)Frequency Domain (20k ft)Time Domain (20 kft) (c) M ∞ = 1 . Fig. 2 Comparison of the aeroelastic eigenvalues predicted by the frequency and time domainmethods for the Agard wing
The fluid computational domain consists of an unstructured grid with , vertices corre-sponding to more than , dofs. Part of the fluid mesh is shown in Figure 7.
2. Flutter Analysis
A flutter analysis is carried on for the three Mach numbers M ∞ ∈ { . , . , . } and altitudes h ∈ { , , } ft. For each combination of flight conditions, the proposed approach is appliedto predict the aeroelastic eigen-characteristics of the system.Hence, for each case, eight sets of real and imaginary parts for the eigenvalues are computed.The eigenvalues are reported in Figure 8(a) for aeroelastic predictions at sea level, in Figure 8(b) for h = 10 , ft and in Figure 8(c) for h = 20 , ft, together with predictions in the time domain.One can observe that the time domain approach is again associated with large approximation errors13 Altitude D a m p i n g R a t i o M ∞ = 0 . Frequency DomainTime Domain (a) M ∞ = 0 . − Altitude D a m p i n g R a t i o M ∞ = 0 . Frequency DomainTime Domain (b) M ∞ = 0 . Altitude D a m p i n g R a t i o M ∞ = 1 . Frequency DomainTime Domain (c) M ∞ = 1 . Fig. 3 Comparison of the smallest aeroelastic damping ratios predicted by the frequency andtime domain methods for the Agard wing of the aeroelastic eigenvalues, especially at M ∞ = 0 . and M ∞ = 1 . . In turn, the minimaldamping ratio is often mispredicted by this time-domain approach,as observed in Figure 9. Theaeroelastic analysis for M ∞ = 0 . at sea level takes approximatively . minutes on 32 processors fora tolerance (cid:15) = 10 − . Computing the steady-state takes . minute and the iterative eigenanalysisapproach takes . minutes.Finally, the flutter speed indices associated with the wing at sea level altitude are computed andreported in Figure 10. One can observe that flutter is never reached, even at sea level, confirmingthe proper flutter-free design of the ARW-2 wing.14 .5 0.6 0.7 0.8 0.9 1 1.10.250.30.350.40.450.50.550.60.65 M ∞ F S I Wing FSISea Level Condition
Fig. 4 Variations with the Mach number of the Flutter Speed Index and comparison with thesea level conditions for the Agard wingFig. 5 Geometry of the ARW-2 wing
3. Design Optimization
Next, the proposed aeroelastic analysis approach is included in a design optimization procedurein order to design an aircraft system under specific dynamic and static aeroelastic constraints. Forthat purpose, the ARW-2 wing is considered as a baseline design at M ∞ = 0 . and sea levelaltitude. In this study the first m = 6 structural modes are considered.The following optimization15 ode number Frequency (Hz) .
922 28 .
833 33 .
274 57 .
045 73 .
706 86 .
087 93 .
918 110 . Table 3 ARW2 retained modes for the structureFig. 6 Finite element model of the ARW-2 wing (a) Surface and symmetry plane (b) Fluid surface mesh
Fig. 7 ARW-2 wing fluid mesh − − − − − − − − Real Part I m ag i n a r y p a r t M ∞ = 0 . Frequency Domain (sea level)Time Domain (sea level)Frequency Domain (10k ft)Time Domain (10k ft)Frequency Domain (20k ft)Time Domain (20 kft) (a) M ∞ = 0 . − − − − − Real Part I m ag i n a r y p a r t M ∞ = 0 . Frequency Domain (sea level)Time Domain (sea level)Frequency Domain (10k ft)Time Domain (10k ft)Frequency Domain (20k ft)Time Domain (20 kft) (b) M ∞ = 0 . − − − − − − Real Part I m ag i n a r y p a r t M ∞ = 1 . Frequency Domain (sea level)Time Domain (sea level)Frequency Domain (10k ft)Time Domain (10k ft)Frequency Domain (20k ft)Time Domain (20 kft) (c) M ∞ = 1 . Fig. 8 Comparison of the aeroelastic eigenvalues predicted by the frequency and time domainmethods for the ARW-2 wing problem is considered max µ ∈D L ( µ ) D ( µ ) s.t W ( µ ) ≤ W max σ ( µ ) ≤ σ max η ( µ ) ≥ η min , (29)where µ ∈ R N µ denotes a vector of N µ design variables and D ⊂ R N µ is the design space. L ( µ ) and D ( µ ) respectively denote the lift and drag associated with the static aeroelastic equilibrium for thedesign configuration µ . W ( µ ) denotes the weight associated with the design µ . The vector σ ( µ ) contains the nodal values of the von Mises stress in the structure at that configuration. Finally,17 Altitude D a m p i n g R a t i o M ∞ = 0 . Frequency DomainTime Domain (a) M ∞ = 0 . Altitude D a m p i n g R a t i o M ∞ = 0 . Frequency DomainTime Domain (b) M ∞ = 0 . Altitude D a m p i n g R a t i o M ∞ = 1 . Frequency DomainTime Domain (c) M ∞ = 1 . Fig. 9 Comparison of the smallest aeroelastic damping ratios predicted by the frequency andtime domain methods for the ARW-2 wing η ( µ ) is the vector of m modal damping ratios associated with the configuration µ for the flightcondition of interest. W max , σ max and η min provide respective bounds on the weight, stress anddamping ratios.Tthere are N µ = 6 design variables. The first three are associated with the shape of the wingand respectively control the sweep, twist and dihedral angles. The last three variables control thethickness of three distinct subsets of spar and rib elements in the wing. The design space D is abox domain defined in Table 4. Note that all six design variables are non dimensional and scaled.The nonlinear program (29) is solved with Matlab’s fmincon active set method with convergencetolerances set to − . The gradients are computed analytically except for the gradient of thedamping ratio that is computed by finite differences. This requires a total of N µ + 1 = 7 calls tothe iterative procedure introduced in this paper for each design considered.18 .4 0.6 0.8 1 1.20.250.30.350.40.450.50.550.60.650.7 M ∞ F S I Wing FSISea Level Condition
Fig. 10 Variations with the Mach number of the Flutter Speed Index and comparison withthe sea level conditions for the ARW-2 wing
Design variable Nature Lower bound Upper bound s Sweep angle − . . s Twist angle − . . s Dihedral angle − . . s Element thickness − . . s Element thickness − . . s Element thickness − . . Table 4 Design domain for the ARW-2 wing
The initial design is the one associated with the original ARW-2 wing. Its physical propertiesare reported in Table 5. The bounds in (29) are set to W max = 400 lbs, σ max = 2 . × and η min = 2 × − . As such, the dynamic aeroelastic constraint is violated for the initial design. Theoptimized design is found after 13 major iterations and reported in Table 5 and the two designconfigurations depicted in Figure 11. Unlike the original design, the damping ratio constraint issatisfied and active for the optimized design. The lift/drag ratio has increased as well as the weight.The maximum von Mises stress has decreased.The overall design procedure takes hours on 32 CPUs. About hour is associated with the19 roperty Initial design Optimized design s − . s − . s . s . s . s . L ( µ ) D ( µ ) 1 .
283 1 . W ( µ ) (lbs) . . σ ( µ )) (psi) . × . × min( η ( µ )) 1 . × − . × − Table 5 Properties of the initial and optimized designs computation of the multiple static constraint and its derivatives and hours to compute the dynamicconstraints and their sensitivities.The variation of the design over the optimization iterations is reported in Figure 12 and thecorresponding objective function in Figure 13. V. Conclusions
This paper presents a practical and accurate approach for determining the aeroelastic charac-teristic of an aircraft configuration. The approach relies on a linearized CFD-based formulationtogether with fixed-point iteration procedure. The application of the proposed approach to theaeroelastic analysis and design optimization of two realistic wing configurations illustrates its capa-bility to enables fast and accurate aeroelastic predictions in the subsonic, transonic and supersonicdomains.
VI. Bibliography [1] Dowell, E. H.,
A Modern Course in Aeroelasticity , Springer, 2014.[2] Hassig, H. J., “An approximate true damping solution of the flutter equation by determinant iteration.”
Journal of aircraft , Vol. 8, No. 11, Nov. 1971, pp. 885–889. ig. 11 Comparison of the initial design (in red) and the optimized design (in blue) [3] Ashley, H., “Piston Theory-A New Aerodynamic Tool for the Aeroelastician,” Journal of the Aeronau-tical Sciences (Institute of the Aeronautical Sciences) , Vol. 23, No. 12, Dec. 1956, pp. 1109–1118.[4] Geuzaine, P., Brown, G., Harris, C., and Farhat, C., “Aeroelastic Dynamic Analysis of a Full F-16Configuration for Various Flight Conditions,”
AIAA Journal , Vol. 41, No. 3, 2003, pp. 363–371.[5] Juang, J.-N. and Pappa, R. S., “An eigensystem realization algorithm for modal parameter identificationand model reduction,”
Journal of Guidance, Control, and Dynamics , Vol. 8, No. 5, Sept. 1985, pp. 620–627.[6] Lesoinne, M., Sarkis, M., Hetmaniuk, U., and Farhat, C., “A linearized method for the frequency anal-ysis of three-dimensional fluid/structure interaction problems in all flow regimes,”
Computer Methodsin Applied Mechanics and Engineering , Vol. 190, 2001, pp. 3121–3146.[7] Back, P. and Ringertz, U., “Convergence of Methods for Nonlinear Eigenvalue Problems,”
AIAA Jour-nal , Vol. 35, No. 6, June 1997, pp. 1084–1087.[8] Meyer, E. E., “Application of a new continuation method to flutter equations,” , 1988.[9] Saad, Y. and Schultz, M. H., “GMRES: A generalized minimal residual algorithm for solving nonsym-metric linear systems,”
SIAM Journal on Scientific and Statistical Computing , Vol. 7, No. 3, 1986, − − − − − Iteration s − − − − − Iteration s Iteration s Iteration s Iteration s − − Iteration s Fig. 12 Variations with the optimization iteration of the design variables pp. 856–869.[10] Yates, E. C., “AGARD Standard Aeroelastic Configurations For Dynamic Response - 1 - Wing 445.6,”NASA, 1987.[11] Seidel, D., Sandford, M., and Eckstrom, C., “Measured unsteady transonic aerodynamic characteristicsof an elastic supercritical wing,”
Journal of aircraft , Vol. 24, No. 4, April 1987, pp. 225–230.[12] Sandford, M., Seidel, D., and Eckstrom, C., “Geometrical and structural properties of an aeroelasticresearch wing (ARW-2),”