AIAA 2000-0419 Multidisciplinary High-Fidelity Analysis and Optimization of Aerospace Vehicles, Part 2: Preliminary Results J. L. Walsh, R. P. Weston, J. A. Samareh, B. H. Mason, L. L. Green, and R. T. Biedron NASA Langley Research Center Hampton, VA 23681 38th Aerospace Sciences Meeting & Exhibit 10-13 January 2000 / Reno, NV For permission to copy or republish, contact the American Institute of Aeronautics and Astronautics 1801 Alexander Bell Drive, Suite 500, Reston VA 20191-4344 AIAA-2000-0419 MULTIDISCIPLINARY HIGH-FIDELITY ANALYSIS AND OPTIMIZATION OF AEROSPACE VEHICLES, PART 2: PRELIMINARY RESULTS J. L. Walsh,* R. P. Weston,† J. A. Samareh,‡ B. H. Mason,§ L. L. Green,¶ and R. T. Biedron** NASA Langley Research Center, Hampton, VA 23681 Abstract tational AeroSciences (CAS) team began a multi- disciplinary analysis and optimization software An objective of the High Performance Computing development project. Initially, the focus of the CAS and Communication Program at the NASA Langley project was on the software integration system, or Research Center is to demonstrate multidisciplinary framework, that was used to integrate fast analyses on a shape and sizing optimization of a complete aerospace simplified design application. The sample application vehicle configuration by using high-fidelity finite- has been the High-Speed Civil Transport (HSCT, Fig. 1). element structural analysis and computational fluid Over the years, the CAS project has focused on dynamics aerodynamic analysis in a distributed, progressively more complex engineering applications, heterogeneous computing environment that includes with the application in the present study known as high performance parallel computing. A software HSCT4.0. A companion paper1 summarizes two system has been designed and implemented to integrate previous applications, known as HSCT2.12 and a set of existing discipline analysis codes, some of them HSCT3.5,3 and presents the HSCT4.0 formulation. The computationally intensive, into a distributed computa- HSCT has also been the focus of other research studies tional environment for the design of a high-speed civil (see Refs. 4–10). transport configuration. The paper describes both the preliminary results from implementing and validating In 1997, the sample application11 shifted to more the multidisciplinary analysis and the results from an realistic models and higher fidelity analysis codes. aerodynamic optimization. The discipline codes are Preliminary results from implementing the HSCT4.0 integrated by using the Java programming language and application are presented in this paper. The HSCT4.0 a Common Object Request Broker Architecture com- application objective is to demonstrate simultaneous pliant software product. A companion paper describes multidisciplinary shape and sizing optimization of a the formulation of the multidisciplinary analysis and complete aerospace vehicle configuration by using high- optimization system. fidelity finite-element structural analysis and computational fluid dynamics (CFD) aerodynamic analysis in a distributed, heterogeneous computing Introduction environment that includes high performance parallel computing. To this end, an integrated set of discipline An objective of the High Performance Computing and analysis codes and interface codes has been formulated Communications Program (HPCCP) at the NASA as a distributed computational environment for the Langley Research Center (LaRC) has been to promote design of an HSCT configuration. The analysis part of the use of advanced computing techniques to rapidly the design loop has been implemented into a software solve the problem of multidisciplinary optimization of integration system that is known as CORBA Java aerospace vehicles. In 1992, the LaRC HPCCP Compu- Optimization (CJOpt)12,13 and is based on a Common Object Request Broker Architecture (CORBA)14 compliant software product and the Java™ programming *Senior Research Engineer, Senior Member AIAA language. †Senior Research Engineer ‡Research Scientist, Senior Member AIAA The formulation of the multidisciplinary design §Aerospace Engineer, Member AIAA optimization problem and of its component analyses is ¶Research Scientist, Member AIAA presented in the companion paper.1 The present paper **Research Scientist is focused on the results obtained up to this point. First, an overview of the HSCT4.0 multidisciplinary opti- Copyright©2000 by the American Institute of Aeronautics and mization application is given and the validation test Astronautics, Inc. No copyright is asserted in the United States under cases are described. Next, the overall analysis process Title 17, U.S. Code. The U.S. Government has a royalty-free license is summarized, and results from the two validation to exercise all rights under the copyright claimed herein for Governmental purposes. All other rights are reserved by the cases are presented for many of the processes included copyright owner. in the complete analysis. There follows a discussion of 1 American Institute of Aeronautics and Astronautics the present state of the sensitivity analysis and selected include constraints on fuel volume, ply-mixture ratio, sensitivity analysis results. Finally, the optimization airfoil interior thickness, takeoff ground scrape angle, demonstration problem description and results—based and landing scrape angle. The structural constraints on a nonlinear aerodynamic analysis—are presented. include buckling and stress constraints. The perfor- mance constraints include constraints on range, takeoff field length, landing field length, approach speed, a Overview time-to-climb-to-cruise requirement, and noise. The weight constraints include constraints on various HSCT4.0 Model weight components. The total number of constraints is The HSCT4.0 application considers a realistic on the order of 32,000 with the majority being aircraft concept and is a multidisciplinary application structural constraints. More detail on the constraints is that integrates high-fidelity analyses representing included in the companion paper. aerodynamics, structures, and performance. For the HSCT4.0 application, a realistic model* that was a The HSCT4.0 application has 271 design variables candidate for development as a commercial HSCT is for optimization—244 structural thickness variables used. This model was originally presented in Ref. 15. and 27 shape variables. To limit the number of inde- Other researchers are also investigating the use of pendent structural design variables, the optimization multidisciplinary analyses, but with simple generic model is divided into 61 design variable zones, as HSCT models.4–10 shown in Fig. 3. Each zone consists of several finite elements. Thirty-nine zones are located on the fuselage Figure 2 shows both the linear aerodynamics grid and 22 zones are located on the wing (11 on the upper and the structural finite-element model (FEM) for half surface and 11 on the lower surface). Within each of the symmetric baseline HSCT4.0 model. A surface zone, four structural design variables are used. These grid with approximately 1100 grid points for a linear structural design variables consist of three ply-thickness code (USSAERO)16 and a volume grid with variables (a 0o fiber variable, a 90o fiber variable, and a approximately 600,000 grid points for a nonlinear code variable that sizes the 45o and –45o fibers) and a core (CFL3D)17 are used in combination for the aerodynamic thickness variable. The ply orientations and composite analyses. For efficiency, only the wing and fuselage are laminate stacking sequence are shown in Fig. 4. modeled in the CFL3D calculation, but the corres- ponding linear aerodynamics model includes the tail The 27 shape design variables (see Fig. 5) consist of surfaces; neither aerodynamic model includes the two sets. The first set contains the nine planform engines. An FEM with approximately 40,000 degrees variables shown in Fig. 5a—the root chord C, the r of freedom (DOFs) is used with the structural analysis outboard break chord C, the tip chord C, the semispan 2 3 code (GENESIS®, a product of VMA Engineering†).18 distance to the outboard break B, the leading edge 2 The engines are modeled as masses on beams. Eight sweep of the two outer wing panels SLE and SLE , the 2 3 laterally symmetric load conditions are used—one total projected area of the three wing panels A, and the t cruise load condition, six arising from maneuver fuselage nose and tail lengths L and L. Note that the n t conditions at +2.5g and - 1g, and one representing a taxi root chord also sets the length of the center fuselage condition. The Flight Optimization System (FLOPS)19 section and that the wing semispan variable B is 3 code is used for modeling the aircraft performance. dependent on other planform variables, including the total area. Optimization Problem Description The objective function of the HSCT4.0 optimization The second set of shape design variables (see Fig. problem is to minimize the aircraft gross takeoff weight 5b) consists of control points that define the wing (GTOW) subject to geometry, structural, performance, camber, thickness, twist, and shear at a set of airfoil and weight constraints. The geometry constraints shape definition points. The definition point locations for camber and thickness are identical, and those for *The computational model for this example has been supplied by the shear and twist are identical. The 18 airfoil shape variables for HSCT4.0 are the vertical (z) perturbations Boeing Company and the results are presented without absolute scales in this paper under the conditions of a NASA Langley Property Loan of the camber, thickness, and shear from the wing Agreement, Loan Control Number I922931. baseline shape and the wing twist perturbation from the baseline shape in constant y planes. Note that the †Any use of trademarks or names of manufacturers in this report is for airfoil camber and thickness perturbations are smooth accurate reporting and does not constitute an official endorsement, globally, while the twist and shear perturbations are either expressed or implied, of such products or manufacturers by the National Aeronautics and Space Administration. linear between the line definition points. 2 American Institute of Aeronautics and Astronautics Optimization Formulation checks that the results obtained are reasonable, based on the experience of the discipline experts. Optimization Process - Figure 6 illustrates the optimization procedure, which consists of a multi- The first set of design variable values, known as the disciplinary analysis (Analysis), gradient calculations “Baseline” configuration, is based on the baseline (Sensitivity Analysis), and a gradient-based optimizer FEM. It has zeros for all planform and airfoil shape (Gradient-Based Optimizer). In the diagram, circles are variables, i.e. for all changes from the baseline shape. used to indicate processes (or functions) and arrows When this set of design variable values is used, the show the data that are passed between processes. By baseline shape and structural thicknesses are convention, this paper uses italics for process names. reproduced. The second set of design variable values, Not all data passed between processes are explicitly known as the “Higher Aspect Ratio” (HAR) configu- shown—only enough data to indicate the required ration, results in a planform shape with a higher aspect sequencing among processes. The outer loop in Fig. 6 ratio than the Baseline configuration and has structural represents one design “cycle,” defined as analysis, design variable values that are increased based on the evaluation of the objective function and constraints, following arbitrary schema. If a design variable repre- sensitivity analysis, and optimization. The Analysis is senting a 0o ply or a 45o ply is within 10% of its lower summarized in a later section and described in detail in bound value, the value is increased by approximately the companion paper. 23%. If a design variable representing a 90o ply is within 10% of its lower bound value, the value is Sensitivity Analysis Process - The Sensitivity Analysis increased by approximately 145%. If a design variable process provides the derivatives of the constraints and representing the core is within 10% of its lower bound the objective function. This process is still being value, the value is increased by approximately 355%. formulated. The other design variables are unchanged. The HAR configuration is expected to have an increased weight Gradient-Based Optimizer Process - The Gradient- as well as changes in stress and buckling responses. Based Optimizer process,20 based on a sequential linear programming (SLP) technique, consists of a general- Implementation Status purpose optimization program (CONMIN),21 an approx- The analysis part of the HSCT4.0 multidisciplinary imate analysis that is used to reduce the number of full optimization process is fully operational within the analyses during the optimization procedure, and some CJOpt framework; it is executed on a heterogeneous set minor process steps. The approximate analysis is used of computers linked by a local area network. All the to extrapolate the objective function and constraints component processes and the complete analysis have with linear Taylor series expansions. This extrapolation been validated as described above. The Sensitivity is accomplished by using derivatives of the objective Analysis process for the complete system is still under function and constraints (from the Sensitivity Analysis development. It will link sensitivity derivatives from process) computed from the analysis at the beginning of the component processes. The completion of the each design cycle. Move limits are imposed on the optimization application awaits the sensitivity analysis, design variables during the Gradient-Based Optimizer although the Optimization process has been validated process to control any errors introduced by the linearity by using an aerodynamic analysis (described below) in assumption. place of the complete multidisciplinary analysis. Validation Test Cases Two sets of initial design variables are used to HSCT4.0 Analysis Process Validation validate that the disciplinary analysis processes des- cribed in the companion paper are integrated correctly. The HSCT4.0 Analysis process is formulated as the The term “integrated correctly” means that the values of sequence of processes shown in the data flow diagram the design variables and all quantities derived from the in Fig. 7 and summarized below. More details of the design variable values are passed from one process to formulation are provided in the companion paper. another process correctly. The method of validation is to execute the Analysis process by using the CJOpt Analysis Process Summary system and to compare the output from each process The Analysis process begins at the top of the data with the output from the original stand-alone code for flow in the figure where the design variable values are that process and the same input for that process. The prescribed. First, the Geometry process derives standalone discipline codes had been validated updated geometries and grids, from baseline geometries previously by comparisons with other engineering and grids, for use by later processes. Next, the Weights results. The integrated system validations include process uses the derived FEM and section properties to 3 American Institute of Aeronautics and Astronautics calculate detailed weights and the center of gravity processes is effectively hidden by executing them in (c.g.) locations for specified flight conditions. The parallel with the Displacements, Loads Convergence, weights data are needed before the remaining processes and Stress & Buckling processes. The Nonlinear can be executed. Next, the Nonlinear Corrections Corrections process remains as the most time- process can be executed. Note that the flow lines to this consuming component process. For this reason, it is process are dashed; the dashed lines indicate that the not executed in every design cycle. Without the Nonlinear Corrections process may not be executed in Nonlinear Corrections process, the parallelized some design cycles due to the high computational time Analysis process executes in just 5 hours. requirements. When this process is not run, the most recent nonlinear corrections continue to be used until an The following sections describe the validation of update is available. several component processes of the Analysis process. No results are presented for the Rigid Trim, Next the Rigid Trim process is executed for the Displacements, or Ground Scrape processes. cruise condition to determine the configuration angle of attack and the tail deflection angle that combine to yield Geometry Results a lift equal to the weight, with no net pitching moment. The MASSOUD22 shape parameterization metho- Once the Rigid Trim process has completed, the left dology was developed for parameterizing changes from branch in Fig. 7, comprising the Polars, Performance, a baseline aircraft shape. Benefits of this approach are and Ground Scrape processes, can proceed in parallel the ease of implementation for the parameterization of with the right branch, comprising the Displacements, complex existing analysis models and grids, the Loads Convergence, and Stress & Buckling processes; relatively few shape design variables required, and the the processes in each branch, however, must proceed consistency of the resulting parameterization across all sequentially. disciplines. The MASSOUD method has been successfully demonstrated for aerodynamic shape The Polars process calculates the drag polars that are optimization, including analytical sensitivity derivative used by the Performance process. The Performance computations, with a structured CFD grid23 and with an process uses the FLOPS code to calculate the mission unstructured CFD grid.24 performance metrics. The Ground Scrape process provides constraints so that the aircraft tail will not Figure 8 shows the FEM for the Baseline scrape the ground on takeoff or landing. configuration on the left-hand side and the HAR configuration, as modified by using the MASSOUD The Displacements process uses a finite-element method, on the right-hand side. The heavy solid lines analysis to generate linear static structural displace- represent the baseline locations of the hexahedral solid ments due to applied aerodynamic and inertial loads for elements that control the planform variation. Note that the cruise condition. These displacements are saved as the large number of elements in the FEM are smoothly a reference set for use in the Loads Convergence perturbed by changes in the 27 shape design variables. process. That process performs the aeroelastic trim calculations for the six noncruise load conditions, Weights Results producing the aeroelastically converged loads on the The as-built weight of a component includes the as- aircraft. Lastly, in the Stress & Buckling process, stress built structural weight, plus nonstructural, systems, and buckling constraints are computed for all elements payload, and fuel weights, of which only the structural contained in the 61 design zones on the fuselage and and fuel weights change with the design variables in the wing. current formulation. The as-built structural weight of a component includes both the theoretical FEM structural Table 1 shows typical execution times and computer weight and structural weight increments for production types used for the Analysis process and the component splices, local pad-ups, side-of-body joints, adhesives, processes within it. The Nonlinear Corrections, Polars, paints, materials for damage tolerance, sealants, and and Loads Convergence processes are the most time- fasteners essential in building the aircraft. Non- consuming; the time for each of the other component structural weight items include windows, landing gear processes is very small. If executed serially, the whole doors, access doors, seat tracks, fuel tank baffles, and process would take almost 17 hours of wall clock time. passenger doors, and system attachment fittings.25 The By using coarse-grain parallelism within the Loads theoretical FEM structural weights are calculated by the Convergence process, executing three load conditions GENESIS® code as applied nodal forces due to a (ten aeroelastic iterations each) simultaneously, the time gravity load vector. The remaining weights are calcu- for that process is reduced from 12 to 4 hours. The time lated by empirical weight estimation methods. for the Polars, Performance, and Ground Scrape 4 American Institute of Aeronautics and Astronautics In order to obtain the as-built weights, the FEM is decreases along with the fuel weight, although the divided into the 14 weight calculation regions shown in structural weights increase. Fig. 9. The as-built weight distributions for HSCT4.0 are defined by a set of 42 files, each of which Nonlinear Corrections Results corresponds to one of the 14 regions. Some regions For efficiency, a Nonlinear Corrections process is have more than one associated file. For example, the used that requires at most one computationally four files associated with region 1 (part of the inboard intensive, nonlinear CFD calculation per load condition wing) represent the structural weight distribution, the during each design cycle. The result is used to nonstructural weight distribution, the systems and calculate a nonlinear correction relative to the equipment weight distribution, and the wing main corresponding linear aerodynamics calculation. The landing gear weight distribution. Nonlinear Corrections process uses the CFL3D code in the Euler (inviscid) mode to capture nonlinear Only the theoretical FEM structural weights of aerodynamics. The computed pressure distribution for regions 1 through 5 are changed directly by the design each load condition is transferred to the panels of the variables. The theoretical FEM structural weights of linear aerodynamics grid while maintaining the same regions 8 through 11 change as these regions stretch total normal force and pitching moment. The nonlinear and/or contract to remain consistent with geometric correction is the panelwise difference between the changes in regions 1 and 2. Under the current nonlinear aerodynamics pressure distribution and the formulation, the structural weight does not change in corresponding linear aerodynamics pressure distri- regions 6, 7, 12, and 13, and none of the nonstructural, bution from the USSAERO code. This correction is systems, or payload weights change. The inboard wing applied many times during the iterations of the Loads fuel weight changes in proportion to the available fuel Convergence process. volume as the geometry changes. Also, the fuel weight files are the only ones that differ between the two mass Figure 11 shows the upper surface distribution of the cases (cruise and GTOW) considered in HSCT4.0. nonlinear corrections, expressed as a correction These two mass cases are obtained by adding the pressure coefficient, for the Baseline and HAR payload weight and the appropriate files of fuel weight configuration cruise conditions and the load factor to the operational empty weight (OEW). extremes of 2.5g and –1g. Because the primary aerodynamic changes are accounted for by the linear Figure 10 shows the percent change of the HAR aerodynamic calculations, the differences between the configuration relative to the Baseline configuration for correction distributions are modest between the two several as-built weight items. For comparison, the configurations, despite significant wing shape changes. figure also shows percent changes in three geometric items that are generally significant in empirical weight For the first Analysis process cycle, the nonlinear estimation methods: the wing average thickness-to- aerodynamics corrections are set to zero. Therefore, all chord ratio, taper ratio, and aspect ratio. Note that the the following results use zero nonlinear corrections. outboard wing as-built structural weight, which has the largest relative weight increase (about 19%), is well Polars Results correlated with the wing aspect ratio increase; all other In the Polars process, the cruise lift and drag from weights change by 5% or less. the Rigid Trim process are augmented by calculating a table of drag polars for a range of Mach numbers and of It is important to realize that a change in the aircraft lift coefficients; this table provides input to the structural weight has a relatively small effect on the calculations in the Performance process. The lift- aircraft GTOW, which is to be minimized in the dependent drag at each specified Mach number is optimization. The entire as-built structural weight is determined by interpolating between USSAERO only about 12% of GTOW and 27% of the OEW, which calculations for a set of angles of attack and tail excludes the payload and fuel weights. Almost half of deflection angles to determine the trimmed drag the structural weight is the as-built structural weight coefficients at the specified lift coefficients. This lift- increment, so it is important that this increment is induced drag contribution is then combined with the included in the optimization. The payload, systems, lift-independent drag contributions of skin friction drag, and non-structural weight increments are held constant wave drag,26 and other miscellaneous drag increments in the current formulation. The fuel is the most to obtain the drag polars. Nonlinear corrections are not significant contributor to weight because it is about calculated for the polars because of the computational 30% of the cruise weight and about 46% of the GTOW. time that would be required. The 1g cruise shape is This effect can be seen in Fig. 10, where the cruise used for all of the Polars process calculations. weight is almost unchanged and the GTOW actually 5 American Institute of Aeronautics and Astronautics Figure 12 shows the polar curves for the Baseline Figure 14 shows the convergence history through ten and the HAR configurations. These curves require 320 aeroelastic iterations in terms of the scaled angle of USSAERO executions in addition to the lift- attack and the scaled tail deflection angle for the independent calculations of skin friction, wave, and Baseline and the HAR configurations. As shown in miscellaneous drags. The 1–2% decreases in drag Figs. 14a and 14c, both the angle of attack and the tail values for the HAR configuration are primarily a result deflection angle of the Baseline configuration at the of lower induced drag due to the higher wing aspect 2.5g load condition have essentially converged after ratio. five aeroelastic iterations. For the –1g load condition, more iterations are required. However, for the HAR Performance Results configuration, more than 10 aeroelastic iterations will The Performance process uses the FLOPS code with be required in the current formulation. The number of the input geometry, weight, and drag polar data for the iterations required can be reduced by applying a current geometric configuration to calculate a variety of relaxation factor. mission performance metrics. These metrics include the range, the takeoff and landing field lengths, the Figure 15 shows aeroelastically displaced wing aircraft noise, the excess thrust, and the takeoff and shapes for both the Baseline and the HAR landing speeds (used by the Ground Scrape process). configurations under the same two load conditions. These results correspond to iteration 10 of Fig. 14. The Figure 13 shows the percent change of the HAR left half of each part of the figure shows the displaced configuration relative to the Baseline configuration for FEM grid that results from application of the several of the performance metrics. The excess thrust aerodynamic loads: the right half shows the shown is needed so that the aircraft has more available corresponding displaced linear aerodynamics surface thrust than drag at all times during the climb to cruise. grid that is used to compute the aerodynamic loads. In this example, only the approach and climb-out excess This figure demonstrates that the structural and thrust, second-segment excess thrust, and thrust-to- aerodynamic shapes of the wing are consistent. This weight ratios increase with the change from the consistency indicates that the aeroelastic iteration is Baseline configuration to the HAR configuration. The converged. (The fuselage shapes are not consistent maximum change occurs for the second-segment excess because the fuselage aerodynamic surface is not thrust, which increases by about 3.5%; all other changes displaced in the current formulation, while the FEM are about 1.5% or less. The decrease in range is well fuselage is displaced; the fuselage shapes will be made correlated with the decrease in fuel volume, as consistent in the next phase of the project.) expected; the decrease in landing field length is well correlated with the decreases in the approach speed and Stress & Buckling Results the GTOW (from Fig. 10). The changes in the wing In the Stress & Buckling process, a fuselage cabin loading and thrust-to-weight ratio are as expected, pressure is added to each of the six converged load given the decrease in GTOW (Fig. 10), while the wing conditions from the Loads Convergence process, and reference area and thrust are held constant. the total is multiplied by a 1.5 factor of safety to produce a set of augmented loads. The GENESIS® Loads Convergence Results code uses this set of augmented loads and a taxi load to In the Loads Convergence process, the trimmed compute stress failure indices and stress resultants for aerodynamic loads for the six noncruise load each element in the 61 design zones. conditions, excluding the taxi condition, are determined from an iterative aeroelastic analysis. An aeroelastic Stress and buckling results for both the Baseline and analysis is required because the structural displace- the HAR configurations are presented in Figs. 16–19 ments depend on the aerodynamic pressure loading and for the 2.5g and –1g load conditions. Elements shown the aerodynamic pressures depend on the displaced in white in the plots are not sized. Stress result plots shape of the aircraft. Because of this mutual depen- represent the Hoffman stress failure index (SFI)18 dency between aerodynamic pressures and structural values. The GENESIS® code computes one SFI value displacements, an iterative process is used to determine for each ply in a sized element (a total of eight values the converged aerodynamic loads at each load per element) for each load condition. Only the condition. Results are presented for the two of the six maximum layerwise SFI value in each element is load conditions that represent the load factor extremes retained. Buckling result plots represent a normalized of 2.5g and –1g. Of the six load conditions, these two buckling load factor (BLF), which is defined in the load conditions are the slowest to converge. companion paper. One BLF value is computed per sized finite element per load condition. Retaining one SFI and one BLF value for each element for each of 6 American Institute of Aeronautics and Astronautics seven load conditions would result in 31,640 structural differentiated with the ADIFOR tool. These derivatives constraints. For simplicity in this paper, stress and have been verified for accuracy by comparisons with buckling results are presented as the largest SFI and the finite-difference approximations. Thus, sensitivity largest BLF per design zone for each load condition. In derivatives for the majority of the codes used in the Figs. 16 to 19, SFI and BLF values above the critical Sensitivity Analysis process are expected to be readily level (1.0) indicate that the constraint has been violated. available. In Fig. 16a, only one design zone on the upper wing The GENESIS® source code is not available for surface of the Baseline configuration has a stress ADIFOR processing, leading to the major difficulty in constraint violation for the 2.5g load condition. As obtaining derivatives for the HSCT4.0 application: shown in Figs. 16b and 17b, no stress constraint choosing a method to obtain the total stress and buckling violations occurred in the Baseline or the HAR constraint derivatives. The total stress and buckling configurations for the –1g load condition. Three zones constraint derivatives depend on component derivatives on the wing upper surface and one on the lower surface obtained by differentiating the equilibrium equation for at the inboard-outboard break had moderate stress linear static analysis with respect to the design variable violations for the HAR configuration for the 2.5g load vector V: condition, according to Fig. 17a. ¶ K u+K¶ u = ¶ f (1) According to Fig. 18a, two design zones on the side ¶ V ¶ V ¶ V and one on the lower surface of the fuselage have moderate buckling constraint violations, and one zone where K is the stiffness matrix, u is the vector of nodal on the lower surface of the fuselage has a severe displacements, and f is the applied load vector. buckling constraint violation for the Baseline Normally, in structural optimization, it is assumed that configuration for the 2.5g load condition. In Fig. 18b, constant loads are used, so ¶ f/¶ V = 0, and methods exist three buckling constraint violations occur in the in the GENESIS® code for obtaining the stress and Baseline configuration on the fuselage for the –1g load buckling constraint derivatives based on that assumption. condition. No buckling constraints are violated for the The plan for the HSCT4.0 project is not to assume HAR configuration for either load condition, as shown constant loads, because the trimmed aeroelastic loads in Fig. 19. from the Loads Convergence process are expected to vary with the shape design variables. One method is to This section has summarized the status of the obtain ¶ f/¶ V for the stress and buckling constraint Analysis process validation. The next section discusses derivatives by finite differences; this method can be the status of the Sensitivity Analysis process. computationally intensive for 271 design variables. An alternate, approximate method to incorporate non-zero ¶ f/¶ V is to exploit the modal approach described in Ref. Sensitivity Analysis Results 9. Sensitivity Analysis Process Sample sensitivity results from codes that may be The Sensitivity Analysis process that provides the used in the Sensitivity Analysis process are presented derivatives of the constraints and the objective function next. The analytical sensitivity derivatives have been is currently under development. Because not every verified by comparisons with finite-difference results. analysis is a direct function of the design variables, it is necessary to obtain the constraint and/or objective Geometry Sensitivities function derivatives by chain-ruling the component The Geometry process calculates shape, structural, derivatives. The plan is to use analytical derivatives and miscellaneous geometries. For the shape geometry, whenever possible, either by hand differentiating the exact sensitivity derivatives required by a gradient- equations or by using the automatic differentiation tools based optimizer are defined as the partial derivatives of ADIFOR27–29 and ADIC,30 to obtain the component the geometry model or grid-point coordinates with derivatives from any analysis for which source code is respect to a design variable. The ADIC tool was available. For example, the geometry tools MASSOUD applied to MASSOUD to augment it with these analytic and CSCMDO31 have been differentiated with the sensitivity derivatives. Figure 20 shows the sensitivity ADIC tool. Similarly, representative versions of the derivative, x/ C, of the x coordinate of the FEM r FLOPS code (used in Weights and Performance nodes with respect to root chord. Because the aft processes), the USSAERO code (used for linear fuselage is fixed, x/ C is small near the inboard r aerodynamics calculations), and the CFL3D code (used trailing edge. The maximum x/ C is near the wing r for nonlinear aerodynamics calculations) have been 7 American Institute of Aeronautics and Astronautics tip. Because in the HSCT4.0 formulation1 the total area Nonlinear Analysis Sensitivities is a design variable, the root chord change has a direct The CFL3D code has been differentiated by using impact on the span of the outboard wing panel. Figure the ADIFOR tool. The resulting derivatives have been 21 shows the volume grid and x/ C for the CFD grid. validated by using finite differences.23 This r The behavior for x/ C is the same as for the FEM. differentiated code is utilized in a following section for r an aerodynamic optimization of the HSCT4.0 The structural geometry sensitivities (ply mixture configuration. and airfoil interior thickness sensitivities) are computed analytically. It is expected that the remaining geometry Performance Sensitivities codes can also be augmented to produce analytical The differentiated FLOPS code that is used for derivatives. weights sensitivity is also used to calculate a variety of exact derivatives of performance metrics. Table 3 Weights Sensitivities shows a sample of sensitivity derivatives of normalized The FLOPS code used as part of the Weights process performance metrics. The performance derivatives in has been differentiated with the ADIFOR tool. The this table are normalized by the same method as the resulting code can compute a variety of exact weights derivatives in Table 2. The independent derivatives of empirical weights metrics with respect to parameters (GTOW, AR, and SW) are defined in the a variety of independent parameters including planform preceding Weights Sensitivities section. The dependent and wing section geometry inputs and weight inputs. variables in Table 3 are a combined noise figure of merit (COMBND), the aircraft range (RANGE), the Table 2 shows a sample of normalized weights approach speed (VAPP), the takeoff field length sensitivity derivatives computed by the ADIFOR- (FAROFF), the landing field length (FARLND), the generated FLOPS code. The derivatives are normalized approach/climb-out excess thrust (AMFOR), and the as (x/y)*(¶ y/¶ x), where x is the value of independent second-segment excess thrust (SSFOR). parameter, and y is the value of the dependent variable. The independent input parameters are the fuselage As can be observed in Table 3, the normalized length (XL), the width of the fuselage (WF), the depth derivatives in the sample are generally within a few of the fuselage (DF), the GTOW, the aspect ratio (AR), orders of magnitude, with both positive and negative the reference wing area (SW), and the inboard wing derivatives represented. Though spanning fewer orders leading edge sweep angle (SWEEP). The dependent of magnitude than the computed weights derivatives variables are the wing weight (WING), the fuselage above, the actual computed performance derivatives weight (FUSELAGE), the total structural weight again span several orders of magnitude. (STRUTOT), the total systems weight (SYSTOT), and OEW. Stress & Buckling Sensitivities The codes used to perform structural sensitivity As shown in the table, the normalized derivatives are analyses and post-process the results are currently being generally of the same order of magnitude, with both tested by comparisons with sensitivities computed by positive and negative derivatives represented. The finite differences. Although the GENESIS® code can actual computed derivatives span several orders of been used to generate stress and stress-resultant magnitude. Computation of a derivative matrix for 115 sensitivities for ¶ f/¶ V = 0 (i.e., constant load vector; see independent by 33 dependent variables requires about 4 Eq. 1), the current HSCT4.0 plan is not to use a minutes on a Silicon Graphics R10000™ workstation. constant load vector. Only the sensitivities of the For comparison, an analysis-only execution of the laminate with the largest SFI value for each element are FLOPS code requires about one minute on the same post-processed by the stress sensitivity code. The BLF workstation. Thus in this case, the use of ADIFOR formulas in the companion paper are differentiated by provides a thousand-fold speedup over the use of finite the chain rule with respect to the shape and structural differences. design variables to form the buckling sensitivities. The resulting BLF sensitivity equations include terms for The GENESIS® code, the other major part of the stress resultant sensitivities; the BLF sensitivity code Weights process, does not have the capability to uses the stress-resultant sensitivities generated by the generate the theoretical FEM nodal weight sensitivities. GENESIS® code for these terms. But again, these are A method for calculating the theoretical FEM nodal based on a constant load vector. The constant load weight sensitivities is being formulated. vector assumption will be tested, and if necessary, an alternative sensitivity method, such as that described in Ref. 9, will be developed. 8 American Institute of Aeronautics and Astronautics Aerodynamic Shape Optimization Results and by the ADIFOR tool for the Nonlinear CFD Sensitivity Analysis process. The drag coefficient One aspect of any optimization formulation is the gradient with respect to each of the 27 design variables proper choice for the upper and lower bounds on the was computed within the CFL3D code and passed on to design variable values. For the HSCT4.0 application the optimizer. structural design variables, these upper and lower bounds are easy to determine. However, the Gradient-Based Optimizer Process - The Gradient- appropriate upper and lower bounds for the shape Based Optimizer process is based on the same design variable values are not so readily determined. sequential linear programming process described above Therefore, as a way to help in specifying upper and for the complete HSCT4.0 application. lower values for the shape design variables, as well as to test the Optimization process (Fig. 6) with a CFD Aerodynamic Optimization Results analysis, an aerodynamic shape optimization problem The specific version of the nonlinear aerodynamic was formulated. code used in HSCT4.0 is CFL3Dv4.1hp. This version of the code has been ported to parallel computer Aerodynamic Shape Optimization Formulation architectures, where it has been demonstrated to scale This optimization problem uses the HSCT4.0 set of well with the number of processors in both the function 27 shape design variables in the same Optimization and derivative modes.23 Figure 24 illustrates the process as the complete HSCT4.0 application (Fig. 6), speedup obtained by using multiple processors for but with a small subset of the Analysis processes shown gradient calculations with respect to the 27 shape in Fig. 7. No structures or performance calculations are design variables considered in HSCT4.0. For an used. increasing number of processors, at least up to the maximum of 32 used for the scaling study, a nearly The objective is to minimize the pressure drag on the linear speedup is observed. With 32 processors, the Baseline configuration (again, wing-body only) at the function calculation shows a superlinear speedup. This cruise Mach number of 2.4. The lift is held fixed at the superlinear speedup is not unusual for large baseline cruise level. This information, along with the problems—as the per-processor memory requirements Analysis and Sensitivity Analysis results, is fed to the decrease, the problem fits more completely into cache, optimizer, which determines a new set of design so that the per-processor floating point operation rate variables that is used to start the process again. The improves over that for fewer processors. However, process is repeated until no significant drag reduction is based on past experience, it is expected that for a observed in subsequent optimization cycles. Other than sufficiently large number of processors the speedup upper and lower bounds on the design variables, no will degrade to sublinear as communication time constraints are applied. becomes a larger part of the overall execution time. The speedup obtained by using multiple processors is Analysis Process - Fig. 22 shows the nonlinear crucial to successful use of high-fidelity CFD methods aerodynamic-only Analysis process that replaces the in the optimization process. multidisciplinary Analysis process of Fig. 7. The Geometry process here consists only of the MASSOUD Figure 25 shows the design cycle history for aircraft shape parameterization code combined with the drag, as measured relative to the baseline values. It can CSCMDO31 volume-grid deformation code. The be seen that the drag has been reduced by Nonlinear CFD Analysis process uses the CFL3D code. approximately 7.5% relative to the baseline. With CFL3D executing in parallel on 32 processors, each The inviscid aerodynamic calculations are done on a design cycle required approximately 1 hour of CPU 129 X 177 X 25 grid split into 32-equal sized blocks time, the bulk of which was the CFL3D computation of and run on 32 processors of a Silicon Graphics Origin the 27 gradients. Although the optimizer has not fully 2000™. In this case, lift coefficient was held fixed by converged for this case, the convergence history from using an option in the CFL3D code that adjusts the 20 design cycles suggests that little additional reduction angle of attack to generate a specified lift coefficient. in drag would be obtained from additional design The drag computed within CFL3D was passed on to the cycles. optimizer. Figure 26 shows a comparison of the baseline and Sensitivity Analysis Process - Figure 23 shows the final surface pressures on both the upper and lower aerodynamic shape optimization Sensitivity Analysis surfaces. The planform changes that occurred between process, which uses analytic sensitivity derivatives the initial and final design cycles are also evident. The generated by the ADIC tool for the Geometry process primary effect on the planform has been to increase the 9 American Institute of Aeronautics and Astronautics