Aircraft Design 2 (1999) 207}230
Sensitivity analysis method for aeroelastic aircraft models Anthony A. Giunta*,1 National Research Council/NASA Langley Research Center, 18d West Taylor Street, Mail Stop 139, Hampton, Virginia 23681-2199, USA
Abstract A novel method has been developed for calculating gradients of aerodynamic force and moment coe$cients for an aeroelastic aircraft model. This method is intended for use in preliminary-level aircraft design which typically involves computationally expensive aerodynamic and structural analyses. This method uses the global sensitivity equations (GSE) to express the aero-structural coupling in an aircraft model. In addition, a reduced-order modal analysis approach is employed to condense the coupling bandwidth between the aerodynamic and structural models. Coarse-grained parallel computing is applied to reduce the wall-clock computational time of the expensive aerodynamic analysis needed in this sensitivity analysis method. A supersonic transport aircraft model is examined in this study, subject to Mach 2.4 cruise #ight conditions. Aerodynamic analysis is performed using a NASA-developed Euler/Navier-Stokes solver, and structural analysis is performed using commercial "nite element analysis software. The GSE/modal analysis method is used to compute the sensitivity of the aerodynamic performance of the aircraft subject to perturbations in the angle-of-attack, wing sweep angle, and wing thickness. Good agreement is obtained between gradients computed with the GSE/modal analysis approach and the same quantities computed using a traditional, computationally expensive, "nite di!erence approach. A cost analysis demonstrates that the GSE/modal analysis method is more computationally e$cient than the traditional approach if gradients are needed for two or more aircraft design parameters. Published by Elsevier Science Ltd.
1. Introduction The focus of this work is the development of a computational method that permits the rapid evaluation of sensitivity derivatives for a coupled aero-structural system. This method is intended
1 Postdoctoral Fellow, National Research Council (1997}1999), currently with Sandia National Laboratories, PO Box 5800, Mail Stop 0847, Albuquerque, NM 87185-0847, USA. * Tel.: 505/844-4280; fax: 505/844-9297. E-mail address:
[email protected] (A.A. Giunta) 1369-8869/99/$ - see front matter Published by Elsevier Science Ltd. PII: S 1 3 6 9 - 8 8 6 9 ( 9 9 ) 0 0 0 1 6 - 6
208
A.A. Giunta / Aircraft Design 2 (1999) 207}230
for use with expensive computational #uid dynamics (CFD) and computational structural mechanics (CSM) analysis tools that are often employed in the preliminary phase of aircraft design. In addition, the method developed in this work employs the CFD and CSM solvers as black boxes. Thus, virtually any CFD or CSM code may be used with this sensitivity analysis technique. Even without attempting to calculate sensitivity derivatives, a coupled aero-structural analysis is cumbersome and expensive using high-"delity CFD and CSM models. A typical aero-structural analysis problem may be stated as follows: given an unloaded aircraft "nite element model (i.e., the jig shape), "nd the "nal deformed shape of the model subject to the aerodynamic loads at the #ight conditions of interest. The aero-structural coupling arises due the interaction between the structural de#ections and the aerodynamic loads. This coupling is depicted in Fig. 1. A typical strategy to solve this coupled system is to perform a CFD analysis, pass aerodynamic loads to the CSM model, execute the CSM software, and pass the structural de#ections back to the CFD model. This process is then repeated until steady-state conditions are obtained for the aerodynamic loads and structural de#ections. Sensitivity analysis of this coupled aero-structural system is needed when attempting to quantify variations in parameters of interest (e.g., aerodynamic performance coe$cients) due to perturbations in the system inputs (e.g., aircraft shape parameters, cruise #ight conditions). The traditional approach to sensitivity analysis uses "nite di!erence approximations in which each input parameter is perturbed by a small amount and then a separate iterative solution is performed for the coupled aero-structural system. Clearly this approach is not attractive when there are many design parameters, and when the CFD and CSM evaluations are expensive. An alternative to the traditional "nite di!erence approach is demonstrated in this study which employs Sobieski's global sensitivity equations (GSE) [1] along with modal analysis methods from structural dynamics [2]. The GSE method provides a framework for decomposing the total sensitivity of a general coupled system into component sensitivity derivatives. For the aerostructural system considered here, the total sensitivity is comprised of sensitivity terms for the aerodynamic model, structural model, and aero-structural interaction. The application of modal analysis is used to simplify the aero-structural interaction terms in the GSE. Another important aspect of the GSE/modal analysis method is that the component sensitivity derivatives are independent and may be evaluated in parallel. Signi"cant wall-clock time savings are realized by applying coarse-grained parallel computing to the CFD analyses needed in calculating some of the GSE terms. The result of this study is a sensitivity analysis method that is signi"cantly less computationally expensive than the traditional "nite di!erence approach. In particular, it is shown that for more than two independent design parameters, the GSE/modal analysis method is superior to the traditional "nite di!erence approach. The GSE/modal analysis method is demonstrated on a simple cantilever beam problem, and on a realistic supersonic transport aircraft model. For the cantilever beam problem, sensitivity derivatives calculated using the GSE/modal analysis method are in good agreement with exact, analytical derivatives. For the aero-structural aircraft problem, sensitivity derivatives computed with the GSE/modal analysis method are in good agreement with derivatives calculated with the traditional "nite di!erence approach. Additional information on the GSE method and on mode shape superposition techniques is provided in Section 2. The formulation of the GSE/modal analysis approach is described in Section 3. A simple cantilever beam problem is presented in Section 4 to demonstrate the GSE/modal analysis method. Section 5 covers the modeling and aeroelastic analysis of a
A.A. Giunta / Aircraft Design 2 (1999) 207}230
209
supersonic transport aircraft. Section 6 contains the results of this study which includes sensitivity derivative data for the aircraft model computed using the GSE/modal analysis method. In Section 7 there is a comparison between the computational expense of the GSE/modal analysis method and the traditional "nite di!erence approach. Section 8 contains a summary of this work.
2. Background This study builds on the past research e!orts of Barthelemy et al. [3] and Dovi et al. [4] who employed the global sensitivity equations in design optimization of a supersonic transport aircraft. Related research by Kapania et al. [5] and Eldred et al. [6] used a variation of the GSE method for the sensitivity analysis of a simple forward-swept wing model. These studies employed inexpensive low-"delity analysis methods such as linear aerodynamics (panel codes) and equivalent plate structural models (ELAPS [7,8]). These inexpensive methods permitted the use of "nite di!erence approximations to evaluate the interdisciplinary coupling derivatives in the GSE (described in Section 3). The use of computationally expensive high-"delity analysis software in this study greatly increases the number of interdisciplinary coupling derivatives in the GSE. That is, the number of coupling terms is related to the amount of force and de#ection data exchanged between the aerodynamic and structural models. With the high-"delity CFD and CSM tools employed in this work, the number of coupling terms in the GSE is O(102}103). Thus, a pure "nite difference approach for approximating the partial derivative coupling terms is not computationally a!ordable. Fortunately, there has been considerable research in the aeroelasticity community in developing methods that simplify the coupling between aerodynamic and structural models. Numerous reduced-order modeling approaches are described in the survey papers by Friedmann [9], Livne [10], and Karpel [11]. Recent examples of the application of reduced-order modeling methods are found in the work of Raveh and Karpel [12] and Cohen and Kapania [13]. The reduced-order modeling approach followed in this study employs a linear superposition of basis vectors to approximate the structural de#ections of an aircraft "nite element model. Here, the basis vectors are supplied by a normal modes (eigenvalue) analysis of the "nite element model. The number of basis vectors used in this study is O(101). One advantage of this reduced-order approach is that it supplies analytic expressions for some of the interdisciplinary coupling terms in the GSE. Another advantage of this approach is that it becomes computationally a!ordable to compute the remaining partial derivatives in the GSE using "nite di!erence methods. The derivation and application of this GSE/modal analysis approach is described below.
3. Sensitivity analysis 3.1. Global sensitivity equations Consider a generic coupled aero-structural system depicted in Fig. 1. The coupling between the aerodynamic and structural models is accomplished by calculating external loads, F, on the aircraft
210
A.A. Giunta / Aircraft Design 2 (1999) 207}230
Fig. 1. A model of a coupled aero-structural system.
structure along with the resulting structural de#ections, D. The aerodynamic and structural data needed to compute F and D are obtained from a variety of computer programs along with suitable pre- and post-processing software to transfer the load and de#ection data between the models. The input into the aero-structural system is a vector of independent parameters X. This vector contains aerodynamic design parameters, such as wing planform and shape parameters, as well as structural design parameters, such as internal rib and spar thicknesses. The output quantities from the coupled system include the aerodynamic load distribution on the de#ected structure, F, along with the force and moment coe$cients C , C , and C 0. In addition, the output includes the "nal L D M structural de#ections, D, and the internal stresses in the structure. For an aircraft design application, quantities of interest such as the aerodynamic force and moment coe$cients are needed to ensure that the aircraft will satisfy the speci"ed mission requirements. In addition, the gradients of these quantities are needed to determine the sensitivity of the aerodynamic performance to small perturbations in the design parameters. Analytic expressions for these gradients usually are not available and they are estimated using "nite di!erence approximations. This requires a separate solution of the coupled system for each perturbation of X; a potentially prohibitive undertaking if either the aerodynamic model or the structural model is computationally expensive to evaluate. Following the notation employed by Sobieski [1], the loads and de#ection transfer in the aero-structural system are represented in functional form as F"F (X, D)
(1)
D"D(X, F).
(2)
and
A.A. Giunta / Aircraft Design 2 (1999) 207}230
211
Di!erentiating Eqs. (1) and (2) with respect to the vector of independent parameters, X, yields dF RF RF dD " # dX RX RD dX
(3)
dD RD RD dF " # . dX RX RF dX
(4)
and
Eqs. (3) and (4) are coupled and may be rearranged into a linear system of the form
C
I
RD ! RF
RF ! RD I
DG H G H dF dX
"
dD dX
RF RX RD RX
.
(5)
Eq. (5) is known as the global sensitivity equation (GSE), which provides a convenient framework for grouping related terms in a system of coupled sensitivity equations. Olds [14] provides a useful lexicon for describing the components of the GSE. Following Olds' approach, the matrix on the left-hand side of Eq. (5) is termed the local sensitivity matrix (LSM) and it contains the partial derivatives of the aerodynamic loads and structural de#ections with respect to each other. Note that the vector of independent parameters, X, does not appear in the LSM. The vector on the right-hand side of Eq. (5) is termed the local sensitivity vector (LSV). This vector contains the partial derivatives of the loads and de#ections with respect to the independent parameters. In the LSV there are no coupling derivatives. The vector on the left side of Eq. (5) is termed the total sensitivity vector (TSV) and it contains the total sensitivity derivatives which are the unknown quantities in this system of equations. While Eq. (5) provides a simple expression for the interdisciplinary coupling of an aero-structural system, the partial derivative terms in Eq. (5) may be particularly di$cult to obtain. For example, consider the term RF/RD in Eq. (5). The vector of aerodynamic forces applied to the structural model, F, is computed using an expensive CFD code (assume for this example that suitable force calculation and interpolation methods are available). If one attempts to calculate RF/RD using a "nite di!erence approximation, then a CFD code evaluation is required for each perturbation of the vector D. Even with parallel computing, this approach is not attractive when the length of D is large, e.g., O(102}103), as is typical for "nite element models used in the aircraft industry. For this reason there is motivation to explore methods that may reduce the computational expense of computing the term RF/RD. One approach to this problem is to represent the nodal de#ections using a linear superposition of a set of basis functions. A version of such a reduced basis approach is employed in this study, where the basis functions are the mode shapes of the structural "nite element model. An important aspect of this work is that the reduced basis method used here is employed while treating the CFD and CSM codes as black boxes. That is, access to the source code of the CFD and CSM analysis software is not needed. While this approach may appear restrictive, it is a realistic model of the aircraft design industry in which both commercial software and legacy software are used in the aircraft analysis and design process. An attractive outcome of this restriction is that
212
A.A. Giunta / Aircraft Design 2 (1999) 207}230
successful methods developed by treating the CFD and CSM codes as black boxes will be broadly applicable to aircraft design practices in industry. 3.2. Modal analysis In linear "nite element structural analysis the vector of structural displacements, D, is found by solving the equation F"KD,
(6)
where F is the vector of applied loads and K is the sti!ness matrix assembled from the individual shape functions of the "nite elements. Another common computation is a modal analysis of the "nite element model. This analysis provides the natural frequencies of the structure along with the associated mode shapes. The modal analysis is performed by solving the eigenvalue problem K/"M/j,
(7)
where M is the mass matrix of the "nite element model, and / is an eigenvector (mode shape) corresponding to a particular eigenvalue, j, that satis"es the linear system. Eq. (7) is solved to obtain the desired number of eigenvalues and eigenvectors for the "nite element model. These n eigenvalues are expressed using a diagonal matrix of the form K"[j , j ,2, j ]I, 1 2 n where I is an n]n identity matrix. The n eigenvectors are grouped using the matrix U as
(8)
U"[/ / ,2, / ], (9) 1 2 n where the n columns of U correspond to the n eigenvalues in K. Note that in this study the eigenvectors are scaled to meet the K-orthogonality and M-orthonormality conditions described by Bathe [2], such that UT KU"K
(10)
UT MU"I.
(11)
and
This scaling is performed internally in the CSM solver GENESIS [15] used in this study. The mode shapes in U are used to approximate the structural de#ections, D, where D"Uq
(12)
and q is a vector of unknown scale factors. The vector q is computed using a least-squares approach where q"(UTU)!1UTD.
(13)
Note that it is typical to use 10}30 mode shapes in this approximation. Thus, the length of q can be much smaller than the length of the vector of structural de#ections D.
A.A. Giunta / Aircraft Design 2 (1999) 207}230
213
To take advantage of this reduced basis approach involving q, the Chain Rule is applied to the term RF/RD as follows: RF RF Rq " , RD Rq RD
(14)
where, from Eq. (13) Rq "(UTU)!1 UT. RD
(15)
Thus, Eq. (14) becomes RF RF T !1 T " (U U) U . RD Rq
(16)
The advantage of this approach is that q is O(101) whereas D typically is O(102}103). It is much more reasonable to compute RF/Rq than RF/RD if using "nite di!erence approximations. The other partial derivative term in the local sensitivity matrix of Eq. (5) is RD/RF. The exact value of this derivative is K!1, however, this matrix is not available to the user in many "nite element analysis codes, particularly those that are commercially developed. For this reason, an approximation for RD/RF is needed. Using the modal analysis approach described above, it is possible to compose an explicit expression relating D and F. This is described below. Substituting Eq. (12) into Eq. (6) yields F"KUq.
(17)
Pre- and post-multiplying this expression by UT and U, respectively, gives UT FU"UT KUqU.
(18)
The K-orthogonality condition in Eq. (10) is used to simplify Eq. (18) to UT FU"KqU,
(19)
which reduces to UT F"Kq.
(20)
Rearranging this equation yields an expression for q where q"K~1UT F.
(21)
Finally, substituting Eq. (21) into Eq. (12) gives D"UK~1 UT F.
(22)
With this explicit relationship between D and F, the analytic expression for RD/RF is RD "UK~1 UT RF
(23)
214
A.A. Giunta / Aircraft Design 2 (1999) 207}230
and the local sensitivity matrix in Eq. (5) can be rewritten as
C
D
RF ! (UTU)~1UT Rq
I
LSM"
T
!UK~1U
I
.
(24)
The remaining unknown terms in the global sensitivity equation are RF/Rq in the LSM, along with RF/RX and RD/RX in the LSV. These terms are needed before solving the GSE for the total sensitivity vector. For simplicity, in this study the unknown partial derivative terms are evaluated using "rst order, forward step "nite di!erence approximations. Details on the "nite di!erence approximations are provided in the examples below.
4. Beam example A simple example problem is shown below to demonstrate the GSE/modal analysis method. Consider a cantilever beam of square cross section as shown in Fig. 2. The vertical force on the tip of the beam acts in the positive z-direction with a magnitude that depends on the amount of de#ection. The magnitude of the load, F, in units of pounds, is F"900 lb!(2000 lb/ft)d.
(25)
Using basic mechanics of materials, the vertical de#ection at the tip of the beam is F¸3 d" , 3EI
(26)
where the length of the beam, ¸, is 15 ft, Young's modulus, E, is 2.3]109 lb/ft2 (for titanium alloy Ti-6Al-4V), and the area moment of inertia, I, is 0.00521 ft4 (h"w"0.5 ft). The solution for F and d to this set of coupled equations is F"757.725 lb and d"0.071137 ft.
Fig. 2. Sample problem of a cantilever beam subject to a tip load, where the load magnitude depends on the tip displacement.
A.A. Giunta / Aircraft Design 2 (1999) 207}230
215
4.1. Exact gradients The total derivatives of F and d with respect to an independent variable, e.g., beam length, may be found analytically as dF !1.62]107EI¸2 " +!23.9567 lb/ft d¸ (3EI#2000¸3)2
(27)
8100EI¸2 dd " +0.011978 ft/ft. d¸ (3EI#2000¸3)2
(28)
and
The exact gradients also may be found using the GSE formulation. The GSE for this coupled system is
C
I
RF ! Rd
Rd ! RF
I
DG H G H
dF RF d¸ R¸ " . dd Rd d¸ R¸
(29)
Substitute the following partial derivatives into the GSE: RF "!2000, Rd
(30)
¸3 Rd " , RF 3EI
(31)
RF "0 R¸
(32)
Rd F¸2 " . EI R¸
(33)
and
Solving for the unknown total sensitivity vector terms gives !6000F¸2 dF " +!23.9567 lb/ft d¸ 3EI#2000¸3
(34)
dd 3F¸2 " +0.011978 ft/ft, d¸ 3EI#2000¸3
(35)
and
which, as expected, are identical to the total derivative values computed explicitly.
216
A.A. Giunta / Aircraft Design 2 (1999) 207}230
Table 1 Natural frequencies, eigenvalues, and eigenvectors for the four node beam "nite element model
Frequency Eigenvalue Eigenvector Components Node 1 Node 2 Node 3 Node 4
Mode 1
Mode 2
Mode 3
5.86 Hz 1354.0
36.77 Hz 53388.0
103.69 Hz 424454.6
0.0 0.058156 0.192160 0.351349
0.0 !0.208162 !0.149674 0.352866
0.0 0.262108 !0.230501 0.351239
4.2. Approximate gradients using GSE/modal analysis Using the GSE/modal analysis approach described in Section 3, a four-node, three-element "nite element model was created for the cantilever beam using the commercial CSM code GENESIS [15]. In this model, node points 1}4 are placed at 0.0, 5.0, 10.0 and 15.0 ft, respectively, with the applied load at node d4. The "rst three natural frequencies and mode shapes are listed in Table 1. The partial derivatives in the GSE were calculated using nodes 2}4, since these were free to de#ect. Note that the subscripts used below indicate node number. The partial derivatives of force with respect to de#ection are RF i "0, for i"2,2, 4 and j"2, 3 Rd j
(36)
and RF 4 "!2000. Rd 4 The remaining partial derivatives in the GSE are
(37)
Rd "UK~1U, RF
(38)
RF 2~4 "M0, 0, 0N R¸
(39)
Rd 2~4 "M0.00212893, 0.007618, 0.0143702N. R¸
(40)
and
Note that the term Rd/R¸ was computed using a "rst-order, forward-step, "nite di!erence approximation with a step size of 1.0% (*¸"0.15).
A.A. Giunta / Aircraft Design 2 (1999) 207}230
217
Solving the GSE for the terms in the total sensitivity vector yields dF 4 +!24.2028 lb/ft d¸
(41)
dd 4 +0.012101 ft/ft, d¸
(42)
and
These estimates for the total derivatives dF/d¸ and dd/d¸ agree to within 3.0% of the exact values for the total derivatives.
5. Aircraft analysis example The GSE/modal analysis method is applied to an aircraft analysis and design problem which involves computationally expensive CFD and CSM codes. High-"delity static aeroelastic analysis is performed for a supersonic transport aircraft at Mach 2.4, 1.0g, cruise conditions, at an angle-of-attack, a, of 3.53. The application of the GSE/modal analysis methods enables the e$cient evaluation of gradients of the aerodynamic lift, drag, and pitching moment coe$cients (C , C , L D and C 0 ) with respect to any of the design parameters. These gradients take into account the M aero-structural interaction in the aeroelastic system. 5.1. Aircraft parametric model A parametric wing/fuselage model of a supersonic transport was developed for this study [16] and contains 64 parameters which de"ne the outer mold line (OML) of the wing (i.e., planform and shape parameters). Six of the parameters which de"ne the wing planform are shown in Fig. 3. Other parameters de"ne the thickness, camber, and dihedral at various spanwise locations on the wing. In addition to the 64 OML parameters, there are 40 parameters that de"ne the structural model, including the thickness of the wing skin panels and the sizes of various rib and spar structural elements. 5.2. Aircraft CFD model Surface and volume grids for CFD analysis are generated based on the OML parameters using the code CSCMDO (coordinate and sensitivity calculator for multi-disciplinary design optimization) [17]. The output from CSCMDO is a structured grid with a C}O discretization scheme and dimensions of 121]41]61, in the streamwise, circumferential, and surface normal directions, respectively. The volume grid is divided into two blocks with the interface splitting the wing into upper and lower surfaces. There are approximately 300 000 grid points in the model (see Fig. 4). A CFD analysis of the HSCT model is performed using CFL3D [18] which is a time-dependent, three-dimensional Reynolds-averaged, thin-layer Navier}Stokes #ow solver. For this HSCT application, CFL3D is employed to solve the steady state inviscid #uid #ow equations, at nominal
218
A.A. Giunta / Aircraft Design 2 (1999) 207}230
Fig. 3. Planform variables for the wing of the supersonic transport model.
Fig. 4. A view of one block in the aerodynamic model of the supersonic transport. This grid shows the starboard wing, the x}z plane of symmetry, and the exit plane.
A.A. Giunta / Aircraft Design 2 (1999) 207}230
219
cruise #ight conditions of Mach 2.4, 3.53 angle-of-attack (with respect to the fuselage centerline), and an altitude of 63 000 ft. An initial #ow solution with CFL3D starting from a uniform #ow "eld requires approximately 60 CPU minutes on an SGI workstation with a 250 MHz, IP27 R10000 processor. Subsequent analyses require approximately 30 CPU min through the use of the restart capability in CFL3D. 5.3. Aircraft CSM model A "nite element model of the aircraft structure is generated using the 64 OML parameters and the 40 structural parameters with the aid of software developed by Balabanov [19] for related research involving supersonic transport con"gurations. The wing/fuselage model of the aircraft comprised a "xed number and topology of spar and rib elements, along with wing skin elements (Fig. 5). The layout of the structural elements is based on the OML parameters, whereas the size (e.g., thickness and area) of the "nite elements is speci"ed through the 40 structural parameters. The "nite element model of the supersonic transport con"guration has 226 nodes and 1130 elements with a total of 1254 degrees-of-freedom. Note that due to structural symmetry only the starboard portion of the model is constructed. The FE model contains triangular membrane elements for the fuselage and wing skins, along with rod elements for the spar cap and rib caps. The
Fig. 5. The structural model of the aircraft showing the wing skin elements (port) and the rib/spar elements (starboard).
220
A.A. Giunta / Aircraft Design 2 (1999) 207}230
spar and rib webs are modeled with a combination of shear panels and rod elements. The material for all structural elements is titanium alloy Ti-6Al-4V. The CSM solver GENESIS [15] is used in this study to perform linear structural analysis and modal analysis of the aircraft model. The computational expense for a single GENESIS analysis is approximately one CPU minute on an SGI workstation. 5.4. Loads transfer and structural deformation External aerodynamic loads for the GENESIS model are obtained using the CFD}CSM loads transfer software code FASIT (#uids and structures interface toolkit) developed by Smith et al. [20,21], and currently maintained by the Air Force Research Laboratory at Wright-Patterson Air Force Base. FASIT provides a suite of interpolation schemes for transferring the aerodynamic loads from the CFD model to the CSM model. The thin-plate spline method of Duchon [22] was used in this study, as recommended in the FASIT User's Manual [21]. The loads transferred from FASIT to the CSM surface grid are written in the NASTRAN bulk data format. Since this NASTRAN format is compatible with the GENESIS input format, no translation was needed between FASIT and GENESIS. The computational expense of running FASIT is negligible (less than 10 CPU s). Displacement of the nodes in the CSM model and the "rst 16 mode shapes are computed using GENESIS. A series of Mathematica [23] programs are used to calculate the mode shape scale factors (see Eq. (13)) needed to represent the node displacements as a superposition of the 16 mode shapes. The rationale for choosing 16 mode shapes is described in Section 5.6 below. Additional Mathematica programs calculate the changes in the 64 OML parameters due to structural deformation. A new CFD grid is then generated from the updated list of 64 OML parameters. 5.5. Aeroelastic analysis Aeroelastic analysis is performed by coupling the CFD, CSM, and interpolation codes as depicted in Fig. 6. The box labeled `Geometry Manipulatora in Fig. 6 contains the software used to generate the CFD and CSM parametric models, including the various Mathematica programs described above. The CFD, CSM, and interpolation software used in the study is loosely coupled using UNIX shell scripts and some simple "le manipulation codes. More detail on the software coupling methods is provided in Ref. [16]. Due to the nonlinearity of the aero/structural interaction, the aeroelastic analysis involves an iterative procedure whereby the aerodynamic and structural analyses are performed repeatedly until both the aerodynamic loads and the structural de#ections reach convergence. The convergence criterion used in this loop is based on the z-direction displacement of the leading edge of the wing tip. If the di!erence in this z-displacement value between two successive passes through the loop is less than 0.05 ft, then the aeroelastic analysis is considered to be converged. To reduce the oscillations that occur in this underdamped system, a constant factor underrelaxation method is used to accelerate convergence. This under-relaxation scheme follows the approach of Chipman et al. [24] and Tzong et al. [25]. A value of 0.7 for the relaxation parameter is used for all aeroelastic calculations conducted in this study. The e!ect of this damping parameter is shown in the convergence history plot shown in Fig. 7. This "gure shows the displacement history
A.A. Giunta / Aircraft Design 2 (1999) 207}230
221
Fig. 6. The arrangement of software used to perform static aeroelastic analysis.
Fig. 7. Convergence history of the aeroelastic analysis with and without relaxation.
of the wing tip versus the number of passes through the aeroelastic loop in Fig. 6. Typically, convergence is obtained in 6}10 iterations (2.7 CPU h) when under-relaxation is used as compared to 19 iterations (6.1 CPU h) without the relaxation method. Results for an aeroelastic analysis at the nominal Mach 2.4, 1.0g cruise conditions are shown in Figs. 8 and 9. In Fig. 8 the "nal deformed wing is shown in comparison to the solid outline of
222
A.A. Giunta / Aircraft Design 2 (1999) 207}230
Fig. 8. Orthographic view of the deformed wing (mesh) and the undeformed wing (solid outline). Note the X : > : Z scaling of 1 : 1 : 2 used to show the wing deformation.
Fig. 9. Airfoil sections at the (top-bottom) wing tip, wing break, and wing root for the initial undeformed (dashed) and "nal deformed (solid) wing shapes.
the undeformed wing. Note that the z-axis in this "gure has been scaled to exaggerate the wing de#ection for easier viewing. Fig. 9 shows the aeroelastic deformation of the aircraft wing at the wing tip, wing break (y"33.7 ft), and wing root. At the 1.0g cruise conditions the z-direction (upward) wing tip de#ection is approximately 1.0 ft.
A.A. Giunta / Aircraft Design 2 (1999) 207}230
223
Fig. 10. Convergence study showing the e!ect of varying the number of mode shapes used to approximate the structural deformation of the aircraft "nite element model.
5.6. Mode shape selection criteria For this study 16 mode shapes were found to be more than su$cient to represent the structural deformation of the aircraft "nite element model. This was demonstrated in a convergence study in which the number of modes used in the aeroelastic analysis procedure was varied from 1 to 16. Fig. 10 shows these results as compared to an aeroelastic analysis performed using the exact structural deformations, i.e., without introducing the superposition of mode shapes. The correct "nal wing position was obtained if at least four mode shapes were used in the linear superposition. All the aeroelastic analysis cases shown in Fig. 10 were started from the unde#ected aircraft shape with identical initial conditions.
6. Aircraft sensitivity analysis Once a converged aeroelastic analysis was obtained for the aircraft model, the GSE method was used to estimate gradients of the coupled system with respect to perturbations in the independent design parameters. For this study, gradients of C , C , and C 0 were computed with respect to L D M three di!erent design parameters: (1) X"a, i.e., the angle-of-attack at cruise, (2) X"(t/c) , the "3%!, thickness-to-chord ratio at the wing leading edge break location (see Fig. 3), and (3) X"" , LE the inboard leading edge wing sweep angle. Calculating gradients for this coupled aero-structural system involves computing the terms in the global sensitivity equation with the modi"cation to the local sensitivity matrix shown in Eq. (24). Solving the GSE yields the total derivatives dF/dX and dD/dX, which estimate the change in the aerodynamic loads and structural de#ections due to a perturbation in the parameter X. These gradient vectors are then used to compute the gradients for C , C , and C 0 . This is possible L D M
224
A.A. Giunta / Aircraft Design 2 (1999) 207}230
since the aerodynamic coe$cients have a functional form identical to that of the aerodynamic loads. That is, C "C (X, D), (43) L L along with similar expressions for C and C 0 . Di!erentiating Eq. (43) with respect to the D M parameter X, yields RC RC dD dC L" L# L . RX RD dX dX
(44)
Following the approach shown in Eqs. (14)}(16) yields dD RC RC dC L " L # L (UTU)~1UT . dX RX Rq dX
(45)
6.1. Interdisciplinary coupling terms Computing the interdisciplinary coupling terms in Eq. (24) involves the mode shape and eigenvalue data from the "nite element solver GENESIS. Some of these data are listed in Table 2 which includes the "rst 16 natural frequencies and eigenvalues for the supersonic transport model, along with the mode shape scale factors computed using Eq. (13). Note that these eigenvalue analysis data are obtained from the undeformed "nite element model of the aircraft, i.e., without application of any loads. The matrix of mode shapes, U, has 226 rows and 16 columns. That is, each row corresponds to one of the nodes in the "nite element model, and each column corresponds to a di!erent mode shape. Table 2 Natural frequencies, eigenvalues, and scale factors for the "rst 16 modes of the aircraft "nite element model Mode
Frequency (Hz)
Eigenvalue
Scale factor
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1.75 4.37 7.05 9.21 10.94 14.82 16.52 19.67 24.42 25.33 28.71 30.93 33.48 34.47 37.59 41.37
121.20 754.61 1961.95 3350.76 4727.34 8675.36 10768.64 15277.35 23533.04 25332.50 32539.54 37765.35 44241.15 46894.51 55771.28 67567.14
6.7740953 !1.9055302 0.5120722 0.0584607 0.0303337 !0.0675731 !0.1006403 0.0534937 !0.0464125 0.0233427 !0.0528983 !0.0113462 0.0265745 !0.0036646 !0.0068915 !0.0055217
A.A. Giunta / Aircraft Design 2 (1999) 207}230
225
The partial derivative term RF/Rq in the local sensitivity matrix is calculated by perturbing each mode shape coe$cient by 20%, and then performing a CFD analysis for each of the new aircraft shapes. The 16 CFD analyses are performed using coarse-grained parallel computing on an SGI Origin 2000 computer. That is, the CFD analyses are executed simultaneously, with each CFD analysis assigned to a separate processor on the parallel computer. The SGI Origin 2000 computers at NASA Langley Research Center and NASA Ames Research Center were used in this study. The performance of this parallel computing approach is shown in Fig. 11, where `speedupa is de"ned as the ratio of the total serial computational time to the parallel execution time. The discrepancy between ideal and actual speedup is due primarily to "le transfer operations between the processors and disk storage. With 16 simultaneous CFD analyses on the NASA Langley Origin 2000, an aggregate parallel performance of 380 MFLOPS was achieved. By using parallel computing, 16 CFD analyses can be completed in about the same wall-clock time as is needed for a single CFD analysis. 6.2. Local sensitivity terms The partial derivative terms in the local sensitivity vector on the right-hand side of Eq. (5) must be computed for each design parameter X, before solving the GSE. Some CFD and CSM solvers may provide these partial derivatives based on analytic di!erentiation of the underlying CFD and CSM state equations. If analytic expressions for the partial derivatives are unavailable, "nite di!erence approximations may be used to estimate the partial derivatives. The "nite di!erence approach is used in this study. 6.3. Global sensitivity calculations 6.3.1. Angle-of-attack sensitivity For the angle-of-attack sensitivity calculations the partial derivative term RF/RX"RF/Ra was computed using "nite di!erences. The aerodynamic loads were obtained from the initial aeroelastic analysis results at a"3.53, and from one additional CFD analysis at a"4.03.
Fig. 11. Speedup obtained using coarse-grained parallel execution of CFL3D on the NASA Langley SGI Origin 2000.
226
A.A. Giunta / Aircraft Design 2 (1999) 207}230
The nodal displacements from the CSM code are not explicitly dependent on angle-of-attack. Therefore, RD RD " "M0N. RX Ra
(46)
The GSE was then solved for the unknown total derivatives dF/da and dD/da. Gradients for C , C , and C 0 were computed using Eq. (45), with the terms RC /Rq, RC /Rq, and RC 0 /Rq L D M L D M computed using the results from the 16 CFD analyses for the mode shape perturbations. The gradients of the aerodynamic coe$cients are listed in Table 3. For comparison, the same gradients were computed using "nite di!erences, with an additional expensive aeroelastic analysis at a"4.03. Good agreement is obtained between the GSE-based gradients and the "nite di!erence gradients. 6.3.2. t/c Ratio sensitivity For the thickness-to-chord ratio sensitivity calculations the nominal value of the t/c ratio at the wing break was 2.36%, and the perturbed value of the t/c ratio was 2.60%. The partial derivative terms RF/R(t/c) and RD/R(t/c) were computed using "nite di!erence approximations. The aerodynamic loads for t/c"2.60% were obtained with a single additional CFD analysis. Similarly, the structural de#ections for t/c"2.60% were obtained from one additional CSM analysis. Since the terms in the local sensitivity matrix do not change with respect to perturbations in the design variables there was no need to recompute the mode shape partial derivative term RF/Rq. Thus, the previously computed LSM was reused when solving the GSE for this t/c sensitivity analysis case. The total derivatives dF/d(t/c) and dD/d(t/c) were computed by solving the GSE, and gradients for the aerodynamic coe$cients were estimated using Eq. (45). These values are listed in Table 4 and show good agreement with gradients computed using a "nite di!erence approximation based on the expensive aeroelastic analysis results. 6.3.3. Wing sweep angle sensitivity The sensitivity of the aerodynamic coe$cients to perturbations in the inboard wing sweep angle were computed using the GSE method. The nominal wing sweep for the aircraft was 74.03, and a perturbation of #2.03 was used in this study. A single CFD analysis for a wing sweep of 76.03 was used to create a "nite di!erence approximation for RF/R" . Similarly, the results from a single LE Table 3 Gradients of C , C , and C o due to perturbations in cruise angle-of-attack (X"a). The "nite di!erence values are L D M computed from aeroelastic analyses at a"3.53 (nominal condition) and a"4.03. The gradients are in units of deg~1
dC /dX L dC /dX D dC o /dX M
GSE
Finite di!erence
0.02432 0.002486 !0.02763
0.02404 0.002408 !0.02707
A.A. Giunta / Aircraft Design 2 (1999) 207}230
227
Table 4 Gradients of C , C , and C o due to perturbations in t/c ratio at the wing break. The "nite di!erence values are L D M computed from aeroelastic analyses at t/c"2.36% (nominal value) and t/c"2.60%. The gradients are nondimensional
dC /dX L dC /dX D dC o /dX M
GSE
Finite di!erence
6.317]10~4 3.835]10~4 !0.929]10~3
9.958]10~4 4.033]10~4 !1.446]10~3
Table 5 Gradients of C , C , and C o due to perturbations in the leading edge sweep angle (X"" ). The "nite di!erence L D M LE values are computed from aeroelastic analyses at " "74.03 (nominal condition) and " "76.03. The gradients are in LE LE units of deg~1
dC /dX L dC /dX D dC o /dX M
GSE
Finite di!erence
!0.003385 !0.000224 0.003359
!0.004111 !0.000266 0.004355
CSM analysis were used to estimate RD/R" . As before, the LSM does not change with respect to LE the aircraft design parameters so it is reused in this sensitivity analysis. After solving the GSE for the total sensitivity terms dF/d" and dD/d" , gradients of the LE LE aerodynamic coe$cients were estimated using Eq. (45). These gradients are listed in Table 5.
7. Computational savings with the GSE approach The advantage of computing gradients with the GSE/modal analysis approach rather than a pure "nite di!erence approach is illustrated in Fig. 12. This plot shows the number of CFD evaluations needed to compute gradients of the coupled aero-structural system, versus the number of independent design parameters. Here, the number of CFD evaluation is used as a cost metric since a CFD evaluation is about 30 times more expensive than any other portion of the aeroelastic analysis process. In the pure "nite di!erence approach, computing the gradients for each new parameter requires about 10 CFD analyses. This cost stems from the need to perform an aeroelastic analysis for each perturbation of a design parameter. Thus, the cost of the pure "nite di!erence approach is linear with respect to the number of independent variables and is given by the expression COST "10n , (47) FD 7 where n is the number of independent variables. 7 In the GSE/modal analysis approach there is an initial cost of 16 CFD evaluations, i.e., one CFD evaluation for each mode shape coe$cient perturbation in the partial derivative term RF/Rq.
228
A.A. Giunta / Aircraft Design 2 (1999) 207}230
Fig. 12. Cost of computing total sensitivity derivatives using "nite di!erences (solid) and the GSE/modal analysis approach (dashed).
However, after this initial cost only one new CFD evaluation is needed for each design parameter. Thus, the GSE/modal analysis approach quickly becomes more attractive, from a computational standpoint, if gradients are needed for more than two design parameters. The cost of the GSE/modal analysis approach is COST "n #16. GSE@M0$!7
(48)
Note that the GSE approach without modal analysis would require a separate CFD evaluation for each of the elements of the partial derivative term RF/RD. With the models used in this aeroelastic analysis, the cost of a GSE-alone approach would be "n #226, COST GSEv!-0/% 7
(49)
where 226 is the number of nodes in the "nite element model. For larger, more realistic "nite element models the value 226 would grow to be O(103}104). In theory, the application of coarse-grained parallel computing renders the wall-clock time identical for both the GSE/modal analysis method and the GSE-alone method. However, the burden of "le management and serial "le input/output make the GSE-alone approach unattractive, even when a parallel computer having hundreds or thousands of processors is available. The GSE/modal analysis approach o!ers signi"cant opportunities for multi-level parallelization. For example, "ne-grained parallel computing could be used to perform a CFD analysis on an aircraft model having several million grid points. With such a large grid the computational domain may be subdivided into hundreds of zones, each of which is assigned to a separate processor on a parallel computer. The coarse grained GSE/modal analysis approach would then provide an additional level of parallelism. As an illustration of this multi-level paradigm, consider a large CFD grid having 100 zones, each of which executes simultaneously on a separate processor of a parallel computer. The GSE/modal analysis approach would replicate this 100-zone grid 16 times, i.e., once for each perturbation of a mode shape coe$cient. Such a strategy would e$ciently utilize 1600 processors. Furthermore, this multi-level strategy readily accommodates increases in the number of zones ("ne-grained scalability) and in the number of mode shapes (coarse-grained scalability).
A.A. Giunta / Aircraft Design 2 (1999) 207}230
229
8. Conclusions A method based on the global sensitivity equations and modal analysis has been developed to calculate gradients of aerodynamic force and moment coe$cients for an aeroelastic aircraft model. The global sensitivity equations capture the the aero-structural coupling in the supersonic transport aircraft model examined in this study. Modal analysis is employed to reduce the coupling bandwidth between the aerodynamic and structural models. Coarse-grained parallel computing is used with the high-"delity computational #uid dynamics solver, CFL3D, for the e$cient calculation of partial derivative terms in the GSE/modal analysis method. Sensitivity analysis was performed for the supersonic transport aircraft model at nominal Mach 2.4 cruise conditions. The GSE/modal analysis approach was used to estimate the gradients of C , C and C 0 with respect to variations in three of the aircraft design parameters. Good L D M agreement was obtained between the GSE-based gradients and "nite di!erence-based gradients of the aerodynamic coe$cients. A comparison of the computational expense for the GSE/modal analysis method and for the traditional "nite di!erence method demonstrates the advantage of using the GSE/modal analysis approach if sensitivity information is needed for two or more aircraft design parameters. Thus, the initial cost of the GSE/modal analysis approach is quickly recovered in a realistic aircraft sensitivity analysis which may involve tens or hundreds of independent parameters.
Acknowledgements This study was performed while the author held a National Research Council Research Associateship at NASA Langley Research Center in Hampton, Virginia. Funding for this work was provided by the Computational Aerosciences Team at NASA Langley, under the direction of Dr. Jaroslaw Sobieski.
References [1] Sobieszczanski}Sobieski J. Sensitivity of complex, internally coupled systems. AIAA Journal 1990;28(1):153}60. [2] Bathe K-J. Finite element procedures. Upper Saddle River, NJ: Prentice-Hall, Inc., 1996, p. 838}45. [3] Barthelemy J-FM, Wrenn GA, Dovi AR, Coen PG, Hall LE. Supersonic transport wing minimum weight design integrating aerodynamics and structures. Journal of Aircraft 1994;31(2):330}8. [4] Dovi AR, Wrenn GA, Barthelemy J-FM, Coen PG, Hall LE. Multidisciplinary design integration methodology for a supersonic transport aircraft. Journal of Aircraft 1995;32(2):290}6. [5] Kapania RK, Eldred LB. Sensitivity analysis of wing aeroelastic response. Journal of Aircraft 1993;30(4):496}504. [6] Eldred LB, Kapania RK. Sensitivity analysis of aeroelastic response of a wing using piecewise pressure representation. Journal of Aircraft 1996;33(4):803}7. [7] Giles GL. Equivalent plate analysis of aircraft wing box structures with general planform geometry. Journal of Aircraft 1986;23(11):859}64. [8] Giles GL. Further generalization of an equivalent plate representation for aircraft structural analysis. Journal of Aircraft 1989;26(1):67}74. [9] Friedmann PP. Renaissance of aeroelasticity and its future. Journal of Aircraft 1999;36(1):105}21. [10] Livne E. Integrated aeroservoelastic optimization: status and directions. Journal of Aircraft 1999;36(1):122}45.
230
A.A. Giunta / Aircraft Design 2 (1999) 207}230
[11] Karpel M. Reduced-order models for integrated aeroservoelastic optimization. Journal of Aircraft 1999;36(1):146}55. [12] Raveh DE, Karpel M. Structural optimization of #ight vehicles with cfd-based maneuver loads. 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, St. Louis, MO, AIAA Paper 98-4832, 1998. [13] Cohen DE, Kapania RK. Trim angle of attack of #exible wings using non-linear aerodynamics. 7th AIAA/ USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, St. Louis, MO, AIAA Paper 98-4833, 1998. [14] Olds J. System sensitivity analysis applied to the conceptual design of a dual-fuel rocket ssto. 5th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, Panama City Beach, FL, AIAA Paper 94-4339, 1994. [15] Vanderplaats, Miura and Associates, Inc., GENESIS User's Manual, Version 4.0. Colorado Springs, CO, 1997. [16] Giunta AA, Sobieszczanski}Sobieski J. Progress toward using sensitivity derivatives in a high-"delity aeroelastic analysis of a supersonic transport. In Proceedings of the 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Design, St. Louis, MO, AIAA Paper 98-4763, 1998, p. 441}53. [17] Jones WT, Samareh JA. A grid generation system for multi-disciplinary design optimization. In Proceedings of the 12th AIAA Computational Fluid Dynamics Conference, San Diego, CA, AIAA Paper 95-1689, 1995, p. 657}69. [18] Krist SL, Biedron RT, Rumsey CL. CFL3D user's manual, version 5.0. NASA/TM-1998-208444, 1998. [19] Balabanov VO. Development of approximations for HSCT wing bending material weight using response surface methodology. Ph.D. Thesis, Virginia Polytechnic Institute and State University, Blacksburg, VA, 1997. [20] Smith MJ, Hodges D, Cesnik C. An evaluation of computational algorithms to interface between CFD and CSD methodologies. Flight Dynamics Directorate, Wright Laboratory, Wright-Patterson Air Force Base, OH, Report WL-TR-96-3055, 1995. [21] Smith MJ, Hodges D, Cesnik C. Fluids and Structures Interface Toolkit (FASIT), Version 1.0., Georgia Tech Research Institute, Atlanta, GA, Report GTRI A-9812-200, 1996. [22] Duchon J. In: Schempp W, Zeller K. editors. Splines Minimizing Rotation-Invariant Semi-Norms in Sobolev Spaces. Springer, Berlin: 1977, p. 85}100. [23] Wolfram S. The mathematica book, 3rd edition. Champaign, IL, 1996. [24] Chipman R, Waters C, MacKenzie D. Numerical computation of aeroelastically corrected transonic loads. In Proceedings of the 20th AIAA/ASME/ASCE/AHS Structures, Structural Dynamics, and Materials Conference, St. Louis, MO, AIAA Paper 79-0766, 1979, p. 178}84. [25] Tzong G, Chen HH, Chang KC, Wu T, Cebeci T. A general method for calculating aero-structure interaction on aircraft con"gurations. In Proceedings of the 6th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, Bellevue, WA, AIAA Paper 96-3982, 1996, p. 14}24.