Probabilistic Design of a Wind Tunnel Model to Match the Response of a Full-Scale Aircraft

Probabilistic Design of a Wind Tunnel Model to Match the

Response of a Full-Scale Aircraft

* † ‡ Brian H. Mason , W. Jefferson Stroud , T. Krishnamurthy*, and Charles V. Spain NASA Langley Research Center, Hampton, VA, 23681-2199 and § Ahmad S. Naser Lockheed Martin Space Operations, Hampton, VA, 23681-2199 An approach is presented for carrying out the reliability-based design of a plate-like wing that is part of a wind tunnel model. The goal is to design the wind tunnel model to match the stiffness characteristics of the wing box of a flight vehicle while satisfying strength-based risk/reliability requirements that prevents damage to the wind tunnel model and fixtures.

The flight vehicle is a modified F/A-18 aircraft. The design problem is solved using reliability-based optimization techniques. The objective function to be minimized is the difference between the displacements of the wind tunnel model and the corresponding displacements of the flight vehicle. The design variables control the thickness distribution of the wind tunnel model. Displacements of the wind tunnel model change with the thickness distribution, while displacements of the flight vehicle are a set of fixed data. The only constraint imposed is that the probability of failure is less than a specified value. Failure is assumed to occur if the stress caused by aerodynamic pressure loading is greater than the specified strength allowable. Two uncertain quantities are considered: the allowable stress and the thickness distribution of the wind tunnel model. Reliability is calculated using Monte Carlo simulation with response surfaces that provide approximate values of stresses.

The response surface equations are, in turn, computed from finite element analyses of the wind tunnel model at specified design points. Because the response surface approximations were fit over a small region centered about the current design, the response surfaces were refit periodically as the design variables changed. Coarse-grained parallelism was used to simultaneously perform multiple finite element analyses. Studies carried out in this paper demonstrate that this scheme of using moving response surfaces and coarse-grained computational parallelism reduce the execution time of the Monte Carlo simulation enough to make the design problem tractable. The results of the reliability-based designs performed in this paper show that large decreases in the probability of stress-based failure can be realized with only small sacrifices in the ability of the wind tunnel model to represent the displacements of the full-scale vehicle.

I. Introduction In the future, advanced aerospace vehicles will be designed using requirements that are expressed in terms of a probabilistic likelihood that the vehicle will achieve a desired performance level. In this paper, that likelihood is referred to as the reliability, and the associated design approach is referred to as reliability-based design (RBD).

Reliability is defined numerically as 1.0 minus the probability of failure to meet a prescribed performance metric.

In a traditional deterministic structural design approach, a factor of safety is applied to reduce the risk of failure.

In that approach, uncertainties are accounted for by designing the structure to carry loads that are substantially greater than the largest loads anticipated. In some cases, this approach can lead to overly conservative designs. Use * Aerospace Engineer, Computational Structures and Materials Branch, Senior Member, AIAA.

† Aerospace Engineer, Computational Structures and Materials Branch, Associate Fellow, AIAA.

‡ Senior Research Engineer, Aeroelasticity Branch.

§ Aerospace Engineer, Aerodynamics, Structures, and Materials Department.

American Institute of Aeronautics and Astronautics

of RBD promises improvements over traditional deterministic design methods by quantifying these uncertainties and accounting for them in the design process. Additionally, in RBD it is possible to determine the sensitivity of the reliability to design changes that can be linked to changes in manufacturing cost. Reliability can then be used as a measure of performance in trade studies. For the same manufacturing cost, it may be possible to design aerospace vehicles that are more reliable than vehicles designed with traditional deterministic methods, or to design vehicles with the same reliability at a lower cost. One of the disadvantages of RBD, however, is an increased computational cost over traditional deterministic methods.

The objective of the present study is to demonstrate methods for design of a plate-like wing of a wind tunnel model to match the stiffness characteristics of the wing box of a full-scale flight vehicle model while satisfying strength-based risk/reliability requirements. The flight vehicle is a modified F/A-18 aircraft. The objective of this study differs from typical aircraft designs (e.g. refs. 1-5), which, instead design for minimum weight. Also, unlike earlier, deterministic work (e.g., refs. 1-4), this design procedure considers several uncertainties and accounts for them by using a probabilistic approach. Stiffness matching of the two models is performed by matching displacements for a specified loading. The design variables in this study define the thickness distribution of an aluminum plate, which simulates a wing box in a wind tunnel model. To address the disadvantage of high computational cost of RBD, this paper describes computational techniques that help make RBD tractable. These techniques include methods referred to as response surface approximation and coarse-grained parallelism.

This paper is organized as follows. Section II provides a description of the wing design problem. In Section III, details of the approaches used in the studies are presented. In Section IV, the results of the deterministic and probabilistic design studies are presented. The deterministic design study was performed to obtain a reasonable starting point for the probabilistic design studies. Finally, a summary of the paper and the conclusions of the design studies are presented in Section V.

II. Description of Wing Design Problem In this section, the wing design problem is described. The finite element models and analyses that simulate the wind tunnel model and its response are discussed, followed by a description of the parameters that are used to change the finite element models in the design process. Then, calculation of the objective function is discussed.

Finally, the deterministic and probabilistic constraints used in the design process are defined.

A. Analysis models The objective of this study is to design a plate-like wing of a wind tunnel model. A discussion of the wind tunnel model (shown in Figure 1) and the experimental test setup is presented in ref. 6. An aluminum plate forms the core of the model, where the wing box would be in the flight vehicle. The control surfaces (4 flaps) are connected to the plate with actuators and hinges. In Figure 1, the aluminum plate (painted green) is covered with end-grain balsa wood to obtain the aerodynamic shape.

A finite element model of the wind tunnel model is shown in Figure 2. In Figure 2, the aluminum plate, the actuators, and the control surfaces are shown in green, blue, and white, respectively. The computer code MSC/NASTRAN (ref. 7) was used to conduct finite element analyses. The structural finite element grid consisted of 88195 nodes, 30246 shell (CTRIA3 and CQUAD4) elements, 46774 solid (CPENTA and CHEXA) elements, 453 beam (CBAR and CBEAM) elements, and 63 rigid (RBAR) elements. Three load cases were used for the analyses.

Load cases 1 and 2 were linear static conditions, and load case 3 was a static aeroelastic load condition. Load case 1 was used to compute bending displacements using the forces shown in Figure 3. Load case 2 was used to compute torsional displacements using the forces shown in Figure 4. Static aeroelastic analyses were performed in MSC/NASTRAN using the structural grid, the aerodynamic vortex lattice grid shown in Figure 5, and load condition 3 (given in Table 1).

Table 1. Static aeroelastic load conditions.

Load case Altitude (ft) Mach Number Dynamic Pressure (psi) 3 5000 0.95 1.607 American Institute of Aeronautics and Astronautics

Aluminum Core Aluminum Core

Balsa Wood Balsa Wood

Actuators Actuators

Figure 1. Wind tunnel model.

60” 60” Aluminum Core Aluminum Core 24” 24” Actuators Actuators Phantom Phantom Fuselage Shape Fuselage Shape Fuselage (beams Fuselage (beams & point masses) & point masses) Figure 2. Finite element model of wind tunnel model.

American Institute of Aeronautics and Astronautics F F Z Z Z F F Y Y Y F F X X X F F F F Clamped Clamped F F Boundary Boundary Condition F Condition F F F F F F F F F F = 108.4 lb F = 108.4 lb f F f F Figure 3. Bending loads applied to the wind tunnel finite element model.

F F F F Z Z Z Y Y Y F F X X X F F F F Clamped Clamped F F Boundary Boundary Condition F Condition F F F F F F F F F F = 108.4 lb F = 108.4 lb f F f F Figure 4. Torsion loads applied to the wind tunnel finite element model.

Figure 5. Aerodynamic grid of wind tunnel model.

American Institute of Aeronautics and Astronautics B. Design variables In the design studies performed in this investigation, the core plate thickness distribution is changed to minimize the objective function subject to deterministic or probabilistic stress constraints. The plate thickness distribution is defined at the finite element nodes by a piecewise linear formulation. In this formulation, 50 thickness design variables are defined in a 5 by 10 grid of points as shown in Figure 6. Nodal thicknesses are interpolated between these 50 grid points. The thickness design variables have upper and lower constraints. The lower constraint is a minimum gage of 0.05 in. The upper constraint is the thickness of the airfoil shape, which varies from 0.51 in. to 2.21 in. as shown in Figure 7.

50 50 49 49 40 40 48 48 47 47 30 30 39 39 46 46 38 38 45 45 20 20 29 29 37 37 44 44 43 43 28 36 28 36 10 19 10 19 42 42 35 35 41 27 41 27 18 18 9 34 9 34 26 26 33 33 17 17 8 25 8 25 32 32 16 16 24 24 31 31 7 7 23 23 15 15 6 6 22 22 14 14 21 21 5 5 13 13 4 12 4 12 3 11 3 11 2 2 1 1 Figure 6. Thickness design variable grid.

t (in.) t (in.)

2.250 2.250 1.375 1.375 0.500 0.500 Figure 7. Airfoil thickness distribution shown on core of wind tunnel model.

actuators and control surfaces removed for visualization C. Objective function The objective of the design studies is to match the stiffness of the aluminum core of the wind tunnel model to the stiffness of the wing box of the full-scale flight vehicle. Stiffness matching is defined in terms of matching displacements from the wind tunnel model to scaled displacements from the full-scale flight vehicle. The displacements on the full-scale flight vehicle were obtained from another finite element model shown in Figure 8.

Scaling is necessary for comparison of the displacements because the wind tunnel model and the full-scale flight vehicle are different in size, and the atmospheric flight conditions are different from the wind tunnel conditions.

Aeroelastic scaling techniques for this model are presented in ref. 6, and those techniques are used in this paper.

American Institute of Aeronautics and Astronautics Aeroelastic scaling must be applied to the bending and torsional forces and to the displacements so that a proper comparison can be made. To scale the displacements of the full-scale vehicle for comparison with the displacements of the smaller wind tunnel model, a length scale factor must be applied to the full-scale displacements. The length scale factor of 0.2609 from ref. 6 is used for this displacement scaling. The forces applied to the full-scale vehicle model are multiplied by a force scale factor of 0.05654 and applied to similar locations in the wind tunnel model as shown in Figures 3 and 4.

The design objective function ( f ) was computed using Equation 1.

F M F M N N − − w w w w i i i i + = (1) f

∑ ∑

F F

( ) ( ) w w max max

= = i i 1 1 i i T B where i =1 to the number of finite element nodes in the plate, w is the nodal displacement in the global z direction of node i , i M designates the wind tunnel model displacements, F designates the full-scale flight vehicle model displacements, B designates the bending load (load condition 1), and T designates the torsion load (load condition 2).

Inspection of Equation 1 shows that the objective function f is effectively the normalized difference in the nodal displacements. During optimization, the wind tunnel displacements are recomputed, and the full-scale flight vehicle displacements are constant. The full-scale displacements were scaled and fit to the wind tunnel finite element model using a piecewise linear interpolation technique, and are shown in Figure 9. The scaled displacements of the flight vehicle were a set of fixed data that remained constant throughout all the design studies.

D. Constraints A stress-based failure constraint is used in all design studies. For the deterministic design study, the constraint is that the von Mises stress (computed from load condition 3) at the centroid of each of the elements in the core must be less than the specified allowable. For the probabilistic design studies, the constraint is that the probability of system failure must be less than the specified value. System failure is assumed to occur if the von Mises stress in any element exceeds an allowable stress. The probability of failure is equal to a count of the number of simulations with system failure divided by the total number of simulations performed. The probabilistic design problem is an example of a series system in which there is a possible failure at multiple locations, any of which can cause the system to fail. Values used in the constraint allowables are discussed next in the design approach section.

Figure 8. Finite element model of the flight vehicle.

American Institute of Aeronautics and Astronautics

Displacement (normalized) Displacement (normalized) 1.0 1.0 0.5 0.5

Bending Bending

0.0 0.0 Displacement (normalized) Displacement (normalized) 1.0 1.0 0.5 0.5

Torsion Torsion

0.0 0.0 Figure 9. Displacements of the flight vehicle mapped onto the wind tunnel model.

III. Design Approaches The design approaches consisted of a combination of a finite element analysis tool (MSC/NASTRAN, refs. 7 and 8), an optimizer (Design Optimization Tools or DOT, ref. 9), and additional NASA-developed codes to perform data transfer, response surface approximations, and probabilistic analysis. This system of codes was combined within framework software called ModelCenter  (ref. 10). In the following subsections, the approaches for deterministic design and probabilistic design are described. The deterministic design is used to obtain a design that can be used as a reasonable starting point for a probabilistic design. Two probabilistic design approaches are discussed. The first probabilistic design approach uses only one random variable, the allowable stress. The second, more computationally expensive design approach uses randomness in the allowable stress and the 50 thickness design variables.

A. Deterministic design approach In the deterministic optimization, the thickness distribution in the core was adjusted to match the bending and torsional displacements of the full-scale model to the wind tunnel model subject to stress constraints. The material ultimate strength of the aluminum plate was 60 ksi. A factor of safety of 4 was used for the deterministic analysis.

The ultimate strength is divided by the factor of safety to produce an allowable stress constraint of 15 ksi. This factor of safety is much larger than values typically used in aircraft design because an additional margin of safety American Institute of Aeronautics and Astronautics was desired to ensure that the wind tunnel model would not suffer from a structural failure and risk damaging the wind tunnel.

Although this margin of safety may appear overly conservative, experience at Langley’s Transonic Dynamics Tunnel (TDT), where the model resulting from this study would presumably be tested, indicates that high safety factors are generally needed, and are not as conservative as they first appear. First of all, the criterion is based on static loading. For flexible models undergoing testing in the TDT, dynamic loading can substantially increase the stress experienced by the model. The dynamics result from naturally occurring turbulence, from interaction of the structures and aerodynamics, and from unpredictable behavior of the model at transonic tunnel conditions.

Although a static aeroelastic model must undergo a flutter clearance analysis, subcritical interactions can still occur, such as low-damped modal responses, which tend to increase model flexing and therefore stress. This increased flexing is especially true at transonic conditions where nonlinear analysis does not always predict static loading or flutter speed accurately.

An implementation of the deterministic optimization in the software framework is illustrated in Figure 10. In the flowchart presented in Figure 10, the Compute Responses phase uses MSC/NASTRAN to compute the objective function and constraints, and the Optimization Search phase uses the optimizer DOT. Each response computation (using load cases 1, 2, and 3) required 15 minutes of CPU time on a single node of a 2.0 GHz Pentium IV Linux cluster. Convergence of the design by having DOT call MSC/NASTRAN directly was very slow, so an alternative approach using approximate analysis techniques was used. This approximate analysis technique utilizes moving response surfaces and is discussed in detail in the Appendix.

START START Current Design Current Design Compute Responses Compute Responses Not Not done done Optimization Optimization Search Search Done Done Optimum Design Optimum Design STOP STOP Figure 10. Functional flow diagram for deterministic optimization.

B. Probabilistic design approach The probabilistic design problem is defined in much the same way as the deterministic design problem. The M objective function (Equation 1) is the same. The displacements w are calculated using the nominal thickness i distribution, with no uncertainty. The moving response surface strategy described in the Appendix is also used. The difference between the two design approaches is in the stress constraint. In the probabilistic approach, stresses are calculated using uncertain (random) values of the thickness in the finite element model rather than nominal values of the thickness. In addition, in the probabilistic approach, the allowable stress is also taken to be random.

Quantification of the randomness in these two quantities is discussed below, followed by a discussion of the two probabilistic studies that were performed.

The first set of random variables specifies the design thickness distribution. A limited amount of manufacturing data was available for the core thicknesses. The manufacturing data consisted of differences between the manufactured and the design thicknesses at 89 points on the core plate for the design in ref. 6. These data are shown as a probability density function (PDF) and a cumulative density function (CDF) in Figures 11 and 12, respectively.

In Figures 11 and 12, a negative thickness difference indicates that the manufactured thickness is less than the design thickness. The CDF shows that 20% of the time, the manufactured thickness is less than the design American Institute of Aeronautics and Astronautics thickness. As an approximation of the randomness observed from the manufacturing data, a combination of three uniform distributions is used in the probabilistic design process. The three uniform distributions are shown in red in Figures 11 and 12. In this approximation, a random thickness increment is added to the nominal thickness ( i.e. the thickness given by the design variables). As shown in Figures 11 and 12, 20% of the time, the random thickness increment is between -0.008 in. and 0.0 in, and another 20% of the time, the random thickness increment is between +0.006 in. and +0.040 in.

0.12 0.12 0.11 0.11 0.10 0.10 0.09 0.09 0.08 0.08 Manufacturing Data Manufacturing Data 0.07 0.07 Approximation (Formula) Approximation (Formula) 0.06 0.06 0.05 0.05 0.04 0.04 0.03 0.03 Probability Density Probability Density 0.02 0.02 0.01 0.01 0.00 0.00 -12 -8 -4 0 4 8 12 16 20 24 28 32 36 40 44 -12 -8 -4 0 4 8 12 16 20 24 28 32 36 40 44 -3 -3 Difference in Manufactured and Design Thickness (10 in.) Difference in Manufactured and Design Thickness (10 in.)

Figure 11. Probability density function for manufacturing data.

1.00 1.00 0.90 0.90 0.80 0.80 0.70 0.70 0.60 0.60 0.50 0.50 Manufacturing Data Manufacturing Data 0.40 0.40 Probability Probability Approximation (Formula) Approximation (Formula) 0.30 0.30 0.20 0.20 0.10 0.10 0.00 0.00 -12 -8 -4 0 4 8 12 16 20 24 28 32 36 40 44 -12 -8 -4 0 4 8 12 16 20 24 28 32 36 40 44 -3 -3 Difference in Manufactured and Design Thickness (10 in.) Difference in Manufactured and Design Thickness (10 in.)

Figure 12. Cumulative probability density function for manufacturing data.

The second random variable is the allowable stress. The allowable stress was assumed to have a normal distribution and a coefficient of variation of 10%. Normally, the factor of safety would be removed from the allowable stress for probabilistic design. However, for this study, the normal distribution represents only randomness in the material properties. The wind tunnel model has a large safety factor of 4 due to uncertainty in turbulence effects in the tunnel, and other factors as described above in the deterministic design approach.

Therefore, for this study, a factor of safety of 3 was used to account for the uncertainties that could not be quantified. Hence, the PDF of the allowable stress is a normal distribution with a mean of 20 ksi and a standard deviation of 2 ksi.

Two probabilistic design case studies are considered in this paper. In the first case, a probabilistic design study is performed using only one random variable, the allowable stress. For this case the system probability of failure (P ) was computed using Equation 2.

f American Institute of Aeronautics and Astronautics x m x   − 1 −   d 2   dx e P = ) Stress Maximum ( (2) f

∫

d 2 π x −∞ = where P (MaximumStress) is the probability of failure computed from the nominal maximum stress, f d is the standard deviation of the allowable stress, m is the mean of the allowable stress, and x is the nominal maximum stress.

The current nominal maximum stress used in Equation 2 is computed from the current nominal values of the design variables using response surfaces. An illustrative example of this approach is presented in Figure 13. In Figure 13, P is computed as the area bounded by the PDF curve of the allowable stress and the nominal maximum f -3 -4 -5 stress. Three cases were studied with different P limits, 10 , 10 , and 10 .

f 0.6 0.6 Change in Change in Computed Maximum Stress, Initial Computed Maximum Stress, Initial design variables design variables Computed Maximum Stress, Designed Computed Maximum Stress, Designed 0.5 0.5 PDF of Allowable Stress PDF of Allowable Stress 0.4 0.4 0.3 0.3 0.2 0.2 Probability Density Probability Density 0.1 0.1 Designed P Designed P f f Initial P Initial P f f 0.0 0.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0 Stress (ksi) Stress (ksi) Figure 13. Probabilistic design with allowable stress as a random variable.

For the second probabilistic design study, random thickness increments were applied to 50 thickness design variables. The random allowable stress was also used in this design study yielding a total of 51 random variables.

This approach is illustrated for a hypothetical design study in Figure 14. In Figure 14, the nominal maximum stress computed from the nominal values of the design variables are represented as vertical lines within the PDFs for the computed load. The PDF of the maximum stress changes with the design variable values during the design process.

The region of overlap between the PDF of the computed load and the PDF of the allowable stress is an indication of the P . The Monte Carlo simulation (MCS) used to compute P in this study is discussed in the next two paragraphs.

f f American Institute of Aeronautics and Astronautics

Change in Change in design variables design variables 1.2 1.2 PDF of Computed Load (Stress), Initial PDF of Computed Load (Stress), Initial PDF of Computed Load (Stress), Designed PDF of Computed Load (Stress), Designed 1.0 1.0 PDF of Allowable Stress PDF of Allowable Stress 0.8 0.8 0.6 0.6 0.4 0.4 Probability Density Probability Density 0.2 0.2 0.0 0.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0 Stress (ksi) Stress (ksi) Figure 14. Probabilistic design with thickness and allowable stress as random variables.

A large number of Monte Carlo simulations are more likely to reflect accurately the distribution of the underlying random variables. As a result, a large number of simulations are more likely to provide an accurate prediction of P . In addition, the optimizer requires derivatives of P with respect to the design variables. These f f derivatives are computed by forward differences using small changes in the design variables. In MCS, P is f computed from an integer count of the number of simulations that exceed the allowable stress divided by the number of simulations performed (N ). Therefore, the resolution of the P calculation is 1.0 divided by N . Small changes S f S in the design variables may not change the P calculation unless N is very large. In this paper, to maintain the f S ** required accuracy, one million MCS trials were used for each estimate of P . Because the responses computed f from the response surface are not accurate for designs outside of the response surface sampling zone, the move limit in Figure A-3 in the Appendix is set to 0.008 in. smaller than the size of the response surface sampling zone.

A Monte Carlo simulation with one million sampling points using a set of linear response surfaces for all 15789 responses required 25 minutes on a single node of a 2.0 GHz Pentium IV Linux cluster. Use of the parallelism software reduced the execution time to 3 minutes. Unfortunately for this design study, DOT required an average of 500 function evaluations for convergence. Thus, the execution time for an optimization within one sampling zone was around 25 hours. Because of this long execution time, only one Monte Carlo design was performed in this -4 paper. For this design study, the constraint was that the P was less than 10 .

f IV. Results of the Design Studies In this section, the results of several optimization studies are presented. In design study 1, a deterministic design was conducted, and its results were compared to the initial design. A probabilistic analysis at the deterministic optimum design was conducted to evaluate the effect of the randomness in the thicknesses on the stresses. In design study 2, the three probabilistic designs were conducted using only one random variable, the allowable stress. In design study 3, a Monte Carlo-based probabilistic design was conducted with 51 random variables.

A. Design study 1 – deterministic optimization For the deterministic optimization study, the design developed for the wind tunnel studies discussed in ref. 6 was used as the initial design. The results of the deterministic optimization are compared to the initial design in Table 2.

For this discussion the three elements with the largest stresses are shown in Figure 15 and referenced in Table 2.

Also, the nodes at the wing tip leading and trailing edge are identified in Figure 15 and referenced in Table 2.

Thickness distributions in the initial and the deterministic optimum designs are shown in Figure 16. The thickness of the core is largest at the root (1.34 inches) and decreases to minimum gage (0.050 inches) at the tip. Differences ** It is possible to estimate the error associated with a Monte Carlo simulation. For example, if a Monte Carlo simulation with one million trials indicates a reliability of 0.9900, there is a 95% probability that the true answer is between 0.9898 and 0.9902 (refs. 11,12).

American Institute of Aeronautics and Astronautics in the displacements of the full-scale and wind tunnel models for the initial and the deterministic optimum designs are shown in Figures 17 and 18. For the bending loads (Figure 3), the differences are shown in Figure 17. For the torsional loads (Figure 4), the differences are shown in Figure 18. Von mises stresses in both the initial and the deterministic optimum designs are shown in Figure 19.

Table 2. Responses computed for initial design and design study 1.

Quantity Initial Deterministic Optimum Stress (Element 1818), ksi 7.424 15.006 Stress (Element 5072), ksi 24.944 14.472 Stress (Element 5073), ksi 19.787 14.107 Objective function f (Equation 1) 609.551 80.128 Difference in Bending Displacement (Node 1), in. 4.420780 0.384090 Difference in Bending Displacement (Node 541), in. 2.691790 -0.068590 Difference in Torsion Displacement (Node 1), in. 1.495337 -0.037647 Difference in Torsion Displacement (Node 541), in. 0.221687 0.015796 Node 1 Node 1 Elem 1818 Elem 1818 Elem 5073 Elem 5073 Elem 5072 Elem 5072 Node 541 Node 541 Figure 15. Elements with high stresses and nodes with large displacements.

Initial Design (ref. 6) Initial Design (ref. 6)

Deterministic Optimum Deterministic Optimum

1.45 0.05 0.75 1.45 0.05 0.75 1.45 0.05 0.75

Thickness (in.) Thickness (in.) Thickness (in.)

Figure 16. Thickness distributions for initial (ref. 6) and deterministic optimum designs.

American Institute of Aeronautics and Astronautics

Initial Design (ref. 6) Initial Design (ref. 6)

Deterministic Optimum Deterministic Optimum

4.55 -0.35 2.10 4.55 -0.35 2.10 4.55 -0.35 2.10

Displacement Difference (in.) Displacement Difference (in.) Displacement Difference (in.)

Figure 17. Bending displacement differences for initial (ref. 6) and deterministic optimum designs.

Initial Design (ref. 6) Initial Design (ref. 6)

Deterministic Optimum Deterministic Optimum

1.30 -0.10 0.60 1.30 -0.10 0.60 1.30 -0.10 0.60

Displacement Difference (in.) Displacement Difference (in.) Displacement Difference (in.)

Figure 18. Torsion displacement differences for initial (ref. 6) and deterministic optimum designs.

American Institute of Aeronautics and Astronautics

Initial Design (ref. 6) Initial Design (ref. 6)

Deterministic Optimum Deterministic Optimum

15.00 1.00 8.00 15.00 1.00 8.00 15.00 1.00 8.00

Stress (ksi) Stress (ksi) Stress (ksi)

Figure 19. Von mises stresses for initial (ref. 6) and deterministic optimum designs.

B. Probabilistic analysis at deterministic optimum An evaluation of the randomness in the stress response due to randomness in the thickness distribution was performed using 100000 Monte Carlo simulations. The PDF and CDF of the maximum stress response at the deterministic optimum design are plotted in Figures 20 and 21, respectively. Note that the PDF of the maximum stress response is much more narrow than the PDF of the allowable stress. The standard deviation of the maximum stress response is 0.374 ksi, which translates into a 2.7% coefficient of variation. The probability of failure as computed by MCS at the deterministic optimum is 1.376e-3.

1.1 1.1 Nominal Stress Nominal Stress 1.0 1.0 PDF of Calculated Maximum Stress due to PDF of Calculated Maximum Stress due to 0.9 0.9 Randomness in Thickness Randomness in Thickness 0.8 PDF of Allowable Stress due to Material 0.8 PDF of Allowable Stress due to Material Randomness Randomness 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 Probability Density Probability Density 0.2 0.2 0.1 0.1 0.0 0.0 12 13 14 15 16 17 18 19 20 21 22 23 24 12 13 14 15 16 17 18 19 20 21 22 23 24 Stress (ksi) Stress (ksi) Figure 20. Probability density functions at deterministic optimum.

American Institute of Aeronautics and Astronautics 1.0 1.0 0.9 0.9 0.8 0.8 0.7 0.7 Nominal Nominal 0.6 0.6 Stress Stress 0.5 0.5 0.4 0.4 Probability Probability 0.3 0.3 0.2 0.2 0.1 0.1 0.0 0.0 12.0 12.5 13.0 13.5 14.0 14.5 15.0 15.5 16.0 12.0 12.5 13.0 13.5 14.0 14.5 15.0 15.5 16.0 Stress (ksi) Stress (ksi) Figure 21. Cumulative probability density function for stress due to randomness in thicknesses at deterministic optimum.

C. Design study 2 – probabilistic optimization with 1 random variable – allowable stress The deterministic design was used as the initial design for design study 2. For this study, Equation 2 was used to compute P using the nominal maximum stress. The allowable stress was a random variable that varied according to f a normal distribution. The moving response surface approach with parallelism (discussed in the Appendix) was used to compute the nominal maximum stress. Convergence required several response surface refits (see Figure A-3 in the Appendix). As the design approached convergence, the response surface size was reduced from 0.080 in. to 0.025 in. to improve the accuracy of the response surfaces. The stresses and objective function values predicted by the response surfaces at the final designs were verified by finite element analysis (FEA).

Three design cases were evaluated in this design study. For each design case, an optimization was performed using a specified P as a constraint. A different P constraint was used for each design case, as shown in Table 3. In f f Table 3, the three largest stresses and the objective function were computed from FEA results at each optimum design. Also, the P for maximum stress was computed using Equation 2 for each case. For comparison, a Monte f Carlo simulation was performed to compute the P using randomness in the 50 thickness variables and the allowable f stress. In all cases, the P computed by Monte Carlo simulation with 51 variables was lower than the P based only f f on the random allowable because the thickness PDF is not symmetric about the nominal value.

Table 3. Responses computed for design study 2.

Probability of Failure (P ) f Quantity Case 1 Case 2 Case 3 -3 -4 -5 10 10 10 P (Maximum Stress) 0.000998 0.000100 0.000011 f P (Monte Carlo analysis) 0.000539 0.000010 0.000002 f Stress (1818), ksi. 13.818 11.225 11.495 Stress (5072), ksi. 13.751 12.499 11.194 Stress (5073), ksi. 13.600 12.560 11.473 Objective function f (Equation 1) 84.227 87.305 89.172 Thickness distributions in the three probabilistic designs are shown in Figure 22. Most variables differed by only a small amount from one case to another. The largest differences in design variables occur at the variables controlling the thickness of elements with high stresses. Thickness design variables 23, 34, and 44 (see Figure 6) differed by more than 0.040 in. from case 1 to case 2. Thickness design variables 23, 34, and 44 (see Figure 6) differed by more than 0.060 in. from case 2 to case 3. Differences in nodal displacements between cases 1, 2, and 3 are very small and are hence not plotted. The value of the objective function increases as the P constraint is reduced f in value as shown in Table 3. Von mises stresses in the three probabilistic designs are shown in Figure 23.

American Institute of Aeronautics and Astronautics Reliability is defined numerically as 1.0 minus P . Therefore, as the reliability requirement is increased, the stresses f are reduced, as shown in Table 3.

-3 -3

P = 10 P = 10

f f -4 -4

P = 10 P = 10

f f -5 -5

P = 10 P = 10

f f

1.45 0.05 0.75 1.45 0.05 0.75 1.45 0.05 0.75

Thickness (in.) Thickness (in.) Thickness (in.)

Figure 22. Thickness distributions for 3 reliability designs with allowable stress as the only random variable.

American Institute of Aeronautics and Astronautics -3 -3

P = 10 P = 10

f f -4 -4

P = 10 P = 10

f f -5 -5

P = 10 P = 10

f f

15.00 1.00 8.00 15.00 1.00 8.00 15.00 1.00 8.00

Stress (ksi) Stress (ksi) Stress (ksi)

Figure 23. Von mises stresses for 3 reliability designs with allowable stress as the only random variable.

D. Design study 3 – probabilistic optimization with 51 random variables For design study 3, both the 50 thickness design variables and the allowable stress were random variables. The -4 constraint imposed for this probabilistic optimization study was that the P computed by MCS was less than 10 .

f From Table 3, the design for case 1 has a P from MCS of 5.39e-4, which is close to the target constraint for this f design study. Therefore, the case 1 design from study 2 was used as the starting point for design study 3.

In design study 3, parallelism was used in generating the response surfaces that the Monte Carlo code used to compute the objective function and the constraint. Convergence of the optimization process required several response surface refits (see Figure A-3 in the Appendix). As the design approached convergence, the response surface size was reduced from 0.080 in. to 0.025 in. to improve the accuracy of the response surfaces. The stresses and objective function values predicted by the response surfaces at the final designs were verified by FEA.

The results of design study 3 are shown in Table 4. In Table 4, the P constraint was computed using one million f Monte Carlo simulations. The three largest stresses and the objective function were computed from FEA results at the optimum design. Note that the designs from cases 1 and 2 in study 2 (Table 3) bracket the optimum design for design study 3. The thickness distribution for the optimum design is shown in Figure 24. Von mises stresses in the optimum design are shown in Figure 25.

American Institute of Aeronautics and Astronautics

Table 4. Responses computed for design study 3.

-4 Quantity Probability of Failure (P ) = 10 f P (Monte Carlo analysis) 0.000098 f Stress (1818), ksi. 11.540 Stress (5072), ksi. 12.747 Stress (5073), ksi. 12.773 Objective function f (Equation 1) 86.253 -4 -4

P = 10 P = 10

f f

1.45 0.05 0.75 1.45 0.05 0.75 1.45 0.05 0.75

Thickness (in.) Thickness (in.) Thickness (in.)

Figure 24. Thickness distributions for reliability design with 51 random variables.

-4 -4

P = 10 P = 10

f f

15.00 1.00 8.00 15.00 1.00 8.00 15.00 1.00 8.00

Stress (ksi) Stress (ksi) Stress (ksi)

Figure 25. Von mises stresses for reliability design with 51 random variables.

V. Summary and Concluding Remarks An approach for carrying out reliability-based design (RBD) of an aluminum plate that represents the wing box of a wind tunnel model is presented. The plate is designed to match the stiffness distribution of the wing box of a full-scale flight vehicle. Stiffness matching was accomplished by matching displacements for a given loading. The plate is also designed to meet stress requirements defined in terms of risk/reliability. Although the design studies in this paper do not address all the issues associated with designing this type of wind tunnel model, these studies can serve as the first step in developing a capability to account for uncertainties in the design of such a wind tunnel model. The design approach can also be used for preliminary design of wing structures for aircraft or spacecraft with thin wings.

The design problem that was studied can be defined as follows. Determine the thickness distribution such that the differences in displacements between the full-scale flight vehicle and the wind tunnel model are minimized and American Institute of Aeronautics and Astronautics

the probability of failure is less than a specified value. Failure is assumed to occur if the stress caused by an aeroelastic loading condition is greater than a specified allowable. There are two types of uncertainties: wing thickness and allowable stress. Data for uncertainties in the thickness distribution of the aluminum wing were obtained from detailed measurements of the experimental wind tunnel model and a comparison of those measurements with the thickness specified. Uncertainties in the allowable stress of the material were assumed to follow a normal distribution.

Three design studies were conducted, with an increasing complexity in the definition of the constraints. In study 1, a deterministic (no random variables) design optimization was conducted with a factor of safety on the allowable stress. The optimum from design study 1 was used as the starting design for design study 2. In study 2, RBD optimization was performed using randomness only on the allowable stress. The probability of failure (P ) in study f 2 was calculated from a computed nominal maximum stress using an equation for normal distribution of the -3 -4 -5 allowable stress. Three design cases with different constraint limits (P = 10 , 10 , and 10 ) were optimized in f study 2 in order to bracket the design space for the more complicated design study 3. In study 3, P was calculated f by Monte Carlo simulation (MCS) with one million trials and using two types of random quantities (allowable stress and thickness). MCS interacted with response surfaces that provided values of stresses. The response surface equations are, in turn, approximations of stresses computed from finite element analysis of the wind tunnel model at specified design points. Because the response surface approximations were fit over a small region centered about the current design, the response surfaces were refit periodically as the design variables changed. Coarse-grained parallelism was used successfully to speed up generation of the response surfaces in all design studies.

From the analyses conducted and presented in this report, the following conclusions were drawn: • The moving response surface approach combined with coarse-grained parallelism is very effective in decreasing the computation time in RBD. Use of response surfaces makes it possible to perform a million MCS in a few minutes.

• Large decreases in the probability of stress-based failure can be realized with only small sacrifices in the ability of the wind tunnel model to match the stiffness distribution of the full-scale vehicle.

• RBD with a small number of random variables can be used to bracket the design space for RBD with a larger number of random variables.

Appendix – Moving Response Surface Approach Response surface approximations were used to provide estimates of the responses that were used to compute the objective function and constraints. An implementation of response surface generation in the software framework is illustrated in Figure A-1. In Figure A-1, a table of sampling point designs was generated from an initial set of design variables by perturbing each variable by a prescribed sampling size, d . The initial design plus 50 thickness design variable perturbations produce a total of 51 sampling point designs. The only variables in the response surface equations are the 50 thickness design variables (see Figure 6). In Figure A-1, a finite element analysis (FEA) is performed at each sample point design using MSC/NASTRAN. In the Build Response Surface phase in Figure A-1, the responses are collected, and response surfaces are fit to the MSC/NASTRAN data. The responses considered included nodal bending and torsional displacements and von Mises stresses evaluated at the centroid of the core shell elements. Each time the Build Response Surface phase is executed, one response surface was generated for each nodal displacement in the z-direction (5358 nodes for 2 load conditions, bending and torsion) and for the stress in each element in the aluminum plate (5073 shell elements). Thus, the total number of sets of response surface approximations was 15789. A linear polynomial was fit to each response surface approximation.

Parallel computing was used to reduce the time required to build response surfaces. This reduction is possible because the ModelCenter framework can execute a NASA-developed code to perform coarse-grained parallel execution of the MSC/NASTRAN analyses. With this code, the user is able to specify the number of concurrent jobs to run simultaneously on different processors of a 2.0 GHz Pentium IV Linux cluster. By using parallelism, it was possible to perform MSC/NASTRAN analyses at all 51 sampling points in 82 minutes using 11 nodes in the cluster, as opposed to 750 minutes on a single node. Thus, by using 11 nodes, all of the MSC/NASTRAN analyses were completed 9.2 times faster than by running sequentially.

The sampling size used in the generation of a response surface has a large effect on the accuracy of the response surface approximation as illustrated in Figure A-2. In Figure A-2, the curve fitted to the larger sampling size (green curve) is less accurate than the curve fitted to the smaller sampling size (red curve). Initially, several sets of response surfaces with different sampling region sizes were generated to study the accuracy of the response surfaces.

Differences between the approximations and the actual values at the sampling points were too large to use the entire design space (minimum gage to maximum airfoil thickness) for a sampling region. As shown in Table A-1, it was American Institute of Aeronautics and Astronautics necessary to use response surface sampling zones no larger than 0.080 in. to have errors in the stress prediction less than 5%.

Table A-1. Response surface sampling size study.

Sampling Size ( d ), in 0.010 0.025 0.050 0.080 0.100 0.150 0.200 % Error in objective function 0.75 0.87 1.08 1.40 1.58 1.27 0.19 % Error in maximum stress 0.35 0.41 0.49 0.62 0.71 1.00 1.22 % Error in 1000 largest stresses 1.86 2.46 3.32 4.65 5.84 23.13 49.58 For the design studies in this paper, multiple response surfaces are used. When the optimizer selects a new design, move limits (upper and lower bounds on the design variables) are adjusted to match the region over which the response surfaces are accurate. The move limit strategy is illustrated in Figure A-3, which shows a hypothetical optimization search in 2-D design space. Once the move limits are reached, new response surfaces are generated, a new move direction is calculated, and new move limits are set. The response surface sampling region size may be increased or decreased at the discretion of the analyst to expedite convergence of the design or to improve the fidelity of the response surfaces, respectively.

START START Current Design Current Design Generate Sampling Point Designs Generate Sampling Point Designs i = 0 i = 0 Design i Design i Finite Element Analysis Finite Element Analysis Yes Yes i = i + 1 i = i + 1 i <= # Samples i <= # Samples No No Build Response Surfaces Build Response Surfaces Response Surface Approximations Response Surface Approximations STOP STOP Figure A-1. Functional flow diagram for response surface generation.

4 d Initial 4 d Initial Design Design d d Figure A-2. Effects of response surface sampling size on response accuracy.

American Institute of Aeronautics and Astronautics

Key: Key: th th th Create n response surface Create n response surface Create n response surface n n n th th th Move limits for n response surface Move limits for n response surface Move limits for n response surface Initial design, Initial design, Initial design, Deterministic optimum Deterministic optimum Deterministic optimum 1 1 1 2 2 2 Second Second 3 3 3 design design variable, variable, 4 4 4

t t

2 2 Final design Final design Final design 5 5 5

First design variable , t First design variable , t

1 1 Figure A-3. Move limit strategy during deterministic and probabilistic optimization.

References Rudisill, Carl S., and Bhatia, Kumar G., “Optimization of Complex Structures to Satisfy Flutter Requirements,” AIAA Journal , Vol. 9, pp. 1487-1491, 1971.

Stroud, W. Jefferson, Dexter, Cornelia B., and Stein, Manuel, “ Automated Preliminary Design of Simplified Wing Structures to Satisfy Strength and Flutter Requirements,” NASA TN D-6534, 1971.

McCullers, L. A., and Lynch, R. W., “Composite Wing Design for Aeroelastic Requirements,” Proceedings of the Conference on Fibrous Composites in Flight Vehicle Design, AFFDL-TR-72-130, 1972.

Haftka, Raphael T., “Automated Procedure for Design of Wing Structures to Satisfy Strength and Flutter Requirements,” NASA TN D-7264, 1973.

Stroud, W. Jefferson, Krishnamurthy, T., Mason, Brian H., Smith, Steven A., and Naser, Ahmad S., “Probabilistic Design of a Plate-Like Wing to Meet Flutter and Strength Requirements,” 43rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Denver, CO, AIAA Paper No. 2002-1464, April 22-25, 2002.

Heeg, J., Spain; C. V., and Rivera, J. “Wind Tunnel to Atmospheric Mapping for Static Aeroelastic Scaling.’ 45th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Palm Springs, CA, AIAA Paper No.

2004-2044, April 19-22, 2004.

Caffrey, John P., and Lee, John N., MSC/NASTRAN Linear Static Analysis, User's Guide , V68, The MacNeal-Schwendler Corporation, Los Angeles, CA, 1996.

Rodden, William P., and Johnson, Erwin H., MSC/NASTRAN Aeroelastic Analysis, User's Guide , V68, The MacNeal- Schwendler Corporation, Los Angeles, CA, 1994.

Anonymous, DOT, Design Optimization Tools, User’s Manual , Version 5.0, Vanderplaats Research & Development, Inc., Colorado Springs, CO 80906.

“Process Integration Using ModelCenter  ,” Technical White Paper, Phoenix Integration, Inc., Blacksburg, VA, 2000.

Ang, Alfredo H.-S., and Tang, Wilson H., Probability Concepts in Engineering Planning and Design. Vol. I, Basic Principles , John Wiley, 1975.

Ang, Alfredo H.-S., and Tang, Wilson H., Probability Concepts in Engineering Planning and Design. Vol. II, Decision, Risk, and Reliability , John Wiley, 1984.

American Institute of Aeronautics and Astronautics