Eng 2

Reliability-Based Design Optimization for Aerospace and Automotive Structures Roberto d’Ippolito ([email protected]), Stijn Donders ([email protected]), Luc Hermans ([email protected]), Herman Van der Auweraer ([email protected]) LMS International, Interleuvenlaan 68, B-3001 Leuven, Belgium Dirk Vandepitte ([email protected]) Katholieke Universiteit Leuven, Department of Mechanical Engineering, Division PMA, Kasteelpark Arenberg 41, B-3001 Heverlee, Belgium

Abstract The need for improvements in engineering designs especially for aerospace and automotive structures is nowadays becoming a major industry request. Deterministic approaches are unable to take into account all the variabilities that characterize design input properties without leading to oversized structures. This paper intends to give a description of the most commonly used methods for reliability based design optimization using the performance measure approach to point out the advantages of the application of this method in the design process. Furthermore, an outline of the methodology will be given that could serve as a guideline to develop a more efficient and optimized design process that takes into account the input parameters variability. For this purpose two Finite Element models have been created, one for an industrial composite wing structure and one for a vehicle knuckle. Realistic variability has been assigned to material properties and static load conditions. A range of analysis cases, with gradually increasing complexity, has then been defined. This allows the assessment of the reported methods in terms of accuracy, computation time and applicability in conjunction with Finite Element models.

1 INTRODUCTION Currently, aerospace and automotive industries are dedicating a lot of attention to improve product quality and reliability already in a virtual simulation environment. The aerospace industry is nowadays using composite materials more and more for both aeronautic and space applications. This trend is shifting the traditional design paradigm dealing with metal structures to a modified paradigm involving hybrid metal/polymer structures. The new design process that is needed to build these structures requires also a shift from the traditional

design approach to a new approach that integrates all the variabilities and uncertainties that hybrid designs introduce, mainly due to the use of composite materials. Also the automotive industry is facing new challenges, especially in the field of durability and reliability of structural components, in order to achieve better performances together with improved safety. Unfortunately, current deterministic approaches are unable to take into account all the variabilities that characterize aerospace and automotive designs, without oversizing structures or assuming a too pessimistic view on the actual material properties and on the operating conditions. Traditionally, safety factors have been applied to the product specifications and dimensions in order to make sure that the structure complies with the internal and external requirements, while being also robust to variability on the design parameters. Unfortunately, these safety factors are often rather arbitrary and typically overestimate the actual variability. This might result in a design with a much higher level of robustness than what is actually required. Especially for (composite) aerospace materials, this could lead to structures that are “heavier than needed”, thus producing too conservative designs that don’t fully exploit the real benefits that lightweight composite materials can bring to the aerospace design process. A more promising approach is to include the property variability in the design flow and to use computer simulation methods to guarantee the design reliability. This paper gives a description of the most commonly used methods to improve the reliability of aerospace and automotive structures. Particularly, a Reliability-based Design Optimization (RBDO) procedure will be illustrated, that uses reliability methods not only to assess the reliability of a given design, but sets a step further: it allows computing a more reliable design for the given product variability.

When accurate stochastic descriptions of physical properties are available, these can be integrated into the design process of aerospace and automotive structures and a procedure management and optimization tool can be developed to • assess the inherent risk of the preliminary design • identify the most relevant parameters • optimize the structure to meet targets on structural performance and safety

2 TEST CASES DESCRIPTION For demonstration purposes, an aerospace and an automotive analysis projects have been created, to improve the reliability of a composite wing subject to gust loads and the reliability of a vehicle knuckle fatigue life, respectively. Finite Element (FE) models are used to represent these structures in a simulation environment and realistic variability has been assigned to material properties and static load conditions. The RBDO procedure has then been used to improve the design of these structures.

2.1 The composite wing A static analysis is performed on a real prototype composite wing in order to optimize the structure for a high altitude, long endurance (HALE) unmanned air vehicle (UAV), for which the characteristics are briefly

approximation is not always true in reality because the realistic characterization of the correlation coefficients is often very difficult and therefore considered impractical for most engineering applications. The composite wing has been examined under gust load conditions, with the wing clamped in the section that connects to the fuselage as boundary condition. The Pratt formula has been used to evaluate the gust load. To characterize a worst-case scenario, a stronger gust situation is considered for most of the presented structural analyses. To stiffen the structure, stringers are distributed over the entire wing. They are considered as onedimensional elements in the FE model, so only the fiber modulus along the fiber direction has been used (Stringer E1, SE1). As stated before, a Finite Element model of the wing has been created to perform the analysis and a particular framework for the necessary computations has been set up. More details on the analysis implementation are given in section 4.1.

2.2 The vehicle knuckle To demonstrate the advantages of the probabilistic methodology applied to fatigue cases, a typical automotive analysis case is considered as well. The vehicle lifetime is highly determined by the fatigue life of its components; in turn, variability in material properties may affect the fatigue life. This is demonstrated for the vehicle knuckle structure shown in Figure 2. The structure is made of steel, and has properties (with inherent variabilities) as specified in Table 2. A fatigue analysis is performed based on a validated FE model; the aim is to increase the fatigue life in the presence of variability.

Figure 1. Main composite wing characteristics and geometry summarized in Figure 1. The composite wing has been made of AS4 12k/3502 unitape. Many properties of this material have already been statistically characterized in the Military Handbook 17-2E and can therefore immediately be used in a probabilistic formulation . The considered input parameters are summarized in Table 1: Property Fiber E1 [Pa] Fiber E2 [Pa] Fiber G12 [Pa] Density [Kg/m3] Stringer Fiber E1 [Pa]

Distrib. Normal Normal Normal Normal Normal

Mean 1.3307•1011 9.3079•109 3.7438•109 1575 1.3307•1011

Variance 7.6381•109 3.9652•108 1.9318•108 2.5 7.6381•109

Table 1. Statistically characterized material properties at 75º F (23.89º C)

be

In this paper, all the input variables are considered to statistically independent. Unfortunately this

(a)

Figure 2: Knuckle finite element model Property

Distrib.

Mean

Elastic Modulus Tensile Strength

Normal Normal

2·105 [MPa] 800 [MPa]

Coefficient of Variance (C.O.V.) 1% 5%

Table 2: Statistically characterized material properties The loads applied on the structure have been obtained from a multibody simulation of an assembly consisting of the structure and the wheels, the steering mechanism and the chassis . All these contributions have been simulated over a time period of 630s and then applied on the different interface points of the knuckle. Then the fatigue effect of the cyclic application of the load histories has been evaluated. The types of loads applied are here summarized: • Wheels load components along the X, Y and Z axes + momentums around axes X and Y • Steering load components along the X, Y and Z axes • Upper chassis connection load components along the X, Y and Z axes • Lower chassis connection load components along the X, Y and Z axes Examples of some of the load time histories are shown in Figure 3.

(b) Figure 3: Wheel load conditions in Newton for X and Y directions (a) and for the Z direction (b) over a time interval of 630s. The analysis case for the knuckle consists of a reliability-based design optimization in order to meet a safety requirement in terms of fatigue life. Details on the procedures used for the optimization will be given in Section 4.2. In general, fatigue is the degradation of a material due to repeated cyclic loading. For metals, this typically means that there is an initiation of small cracks from active slip bands in grains on the free surface of a specimen, component or structure, which eventually link to form large cracks that either break or severely degrade the performance of a component. Two basic approaches have been developed to estimate the crack initiation life of components and to address fatigue analysis and design for durability: the stress-life approach and the strain-life approach. The goal of both approaches is to estimate the crack initiation life of structures, by combining a mechanics-based analysis of the stresses or strains in the structure with the results of basic material property tests. The type of fatigue analysis approach to use will depend on the type of failure mode that is expected for the structure. When repeated fatigue life tests are performed, one does not find exactly repeated results, but rather a scatter of results. This scatter comes from several sources. There are natural variations in material properties, component dimensions, customer service loads and manufacturing tolerances. The largest variation is typically found in the customer usage. One can use probabilistic distributions to

characterize the variation of these variables, and include this variability in the design process to better explain and predict the scatter obtained in fatigue test results. Probabilistic methods offer a powerful tool to assess the inherent risk of designs for durability through the characterization of the inherent variabilities, by allowing the design engineer to minimize the weight and costs requirements, without losing reliability with respect to durability requirements.

3 THE RELIABILITY PROBLEM STATEMENT As stated in Section 1, the variabilities of an engineering design can be characterized by the variations of a random system parameter set x=[ xi ]T

withi =1...n

The probability distribution of each xi is described either by its cumulative distribution function (CDF, Fxi ( xi ) ) or by its probability density function

(PDF, f xi ( xi ) ) and is often bounded by tolerance limits on the system parameter values. The system performance is described by Performance Functions (PF, Gj ( x) with j = 1...m ) that, for a structural design, are usually the selected failure criteria. Thus, with Gj ( x) one of the m system PFs, the system is con-

sidered to fail if Gj ( x) < 0 for at least one index j. The probability of failure of the system for every jth performance function is then

(

)

FGj ( 0) = P Gj ( x) < 0 = ∫ ∫ f X ( x) dx1 dxn Ωj

{

}

The event space Ω j is the region of the stochastic space where only failure events occur. Thus, the integral of the joint probability density function of the system response over Ω j yields the probability of failure Pf , j of the structure for the jth criterion. In general, given a fixed performance index (or measure) g j of the system, Eq (1) can be generalized

( )

(

) ∫ ∫ f X ( x) dx1 dxn

{

ΩG j

n where ΩGj : x ∈ ℜ : Gj ( x) < g j

}

3.1 The Reliability-Based Design Optimization (RBDO) Model In engineering design, the traditional deterministic design optimization model has been successfully applied to systematically improve the system design process, yielding a reduction of the costs and an improvement of the final quality of the products. However, uncertainties in either engineering simulations and/or manufacturing processes exist. This calls for different optimization models that can yield not only an improvement in the design, but also a higher level of confidence. Thus, a reliability-based design optimization (RBDO) model for robust and costeffective designs can be defined using mean values of the random system parameters as design variables and optimizing the cost subject to prescribed probabilistic constraints (e.g. a maximum on the allowed probability of failure) by solving a mathematically nonlinear programming problem. As a result, the RBDO solution provides not only an improved design but also a higher level of confidence in the design. The general RBDO model can be defined as

[ Cost( d) ]  min  d  s.t. Pf , j = P Gj ( x) < 0 ≤ Pf , j with j = 1 m

(

)

(3)

where the cost function can be any function of the design variable d = [ di ]T ≡ [ µi ] T and a probabilistic constraint Pf , j can be defined for each failure mode. Thus, for each

(1)

n where Ω j : x ∈ ℜ : Gj ( x) < 0

FGj g j = P Gj ( x) < g j =

In this case, the event space ΩGj represents the region of the stochastic space where the performance of the structure is below a prescribed quantity g j .

(2)

iteration of the optimization loop, an estimation of the probabilistic constraint has to be computed. For this purpose, different methods exist. In general, each prescribed probability of failure Pf , j can be represented in terms of the target reliability index βt, j as

βt, j = −Φ −1( Pf , j ) ⇒ Pf , j = Φ( − βt, j )

(4)

also defined by the more general relation

( )

(

)

(

FGj g j = P Gj ( x) < g j = Φ − β Gj

)

(5)

where Φ( •) is the standard normal CDF (zero mean and standard deviation 1). Using Eq. (5), the second condition of Eq. (3) can be rewritten as

(

)

(

)

Pf , j = FGj ( 0) = Φ − β s, j ≤ Φ − βt, j = Pf , j

(6)

⇒ β s, j ≥ βt, j

where β s, j is traditionally called “reliability index” of the structure for the specified failure mode j. Through inverse transformations, Eq (6) can be written as

( (

))

g*j = FG−j1 Φ − β t, j ≥ 0

(7)

* where g j is named “target probabilistic performance

measure” and represents the value of the performance function “equivalent” to the target reliability index βt, j .

3.2 Performance Measure Approach The formulation of the general RBDO model of Eq. (3) that uses Eq. (6) to describe the probabilistic constraint is called Reliability Index Approach (RIA) . An alternative and more effective way of formulating the general RBDO model is to use the inverse formulation of the probabilistic constraint of Eq. (7), known as Performance Measure Approach (PMA) . The RBDO model can be re-written as

[ Cost( d) ]  min d  s.t. g*j = FG−1j Φ − β t, j ≥ 0 with j = 1 m  

( (

))

(8)

[ ] [ ] T

T

In this case, at a given design d (k) = di(k) ≡ µi(k) , the evaluation of the target probabilistic performance * measure g j is performed using

    g*j d (k) = FG−1j   f X ( x) dx1 dxn   ΩG   j 

( )

∫ ∫

(9)

For each iteration of the optimization process the following problem has to be solved

[

]

 Gj ( u) with j = 1 m min  u  s.t. u = t,j

(

)

which depends on the design points, has to be estimated within each iteration and the target of the optimization is to find the design point that can satisfy all the probabilistic constraints in terms of the performance measure. This approach yields an improved convergence rate of the PMA methodology with respect to the RIA formulation. In fact, PMA minimizes a complex cost function subject to a simple constraint, which gives rise to less numerical problems and therefore a better algorithm performance, while the RIA algorithm minimizes a simple cost function subject to a complex constraint.

3.3 Limit State Approximations (LSA) In most industrial cases, the integral of Eq. (9) cannot be evaluated in closed form. Only for simple academic cases, fX can be defined analytically. In such a case, the performance function Gj ( x) for each specific failure mode of the structure is available so that Ω j can be defined. The boundary of the domain Ω j is also called

limit-state surface and the function Gj ( x) is called the limit-state function (LSF) . In the most widely used reliability method, approximations are made in the space of standard uncorrelated normal variates, Y, obtained from a transformation of the basic variables Y = T( X )

(12)

where the transformation T is expressed in terms of the distributions of X. For the sake of simplicity, only the case with one failure criterion (j=1) will be considered in this paper. Therefore, from here on, the j index is removed from all the quantities that depend on j. In the space Y, denoted as the standard normal space, approximations of the probability of failure Pf , j of Eq. (3) are obtained by replacing the limit state surface with first or second order approximating surfaces. These surfaces are fitted to the limit state surface at points with minimal distance to the origin (design points).

(10)

The approximate solution of this problem yields an estimation of the performance measure, given by

g*j ≈ G u*β =βt , j

* With this formulation, the performance measure g j ,

(11)

3.4 First and Second Order Reliability Methods In the first order reliability method (FORM), the limit state surface is approximated by its tangent hyperplane in the standard normal space Y at the global minimum distance point (Most Probable Point, MPP) on the limit state function. With this approximation, the reliability index s coincides with the distance of the MPP to the ori-

gin. In the second order reliability method (SORM) the limit state surface is approximated by a paraboloid, defined with the same principal curvatures as estimated for the limit state surface at the design point (MPP). Determination of these curvatures essentially requires computing the second-order derivative matrix of the limit state surface and solving an associated eigenvalues problem. The determination of the MPP involves the solution of the following constrained optimization problem:  1   min f ( y) = yT y  y  2   s.t. G( y) = 0

with y ∈ Y

(13)

To achieve the optimal solution, the Karush-KuhnTucker (KKT) conditions have to be satisfied. For only one equality constraint, the Lagrangian expression and the KKT conditions can be reduced to: L( y, λ ) = f ( y) + λG( y)

(14)

 ∇ y L y , λ = 0  *  G( y ) = 0

(15)

(

*

*

)

*) is the solution point. where (y*, Many algorithms are available in literature to solve this problem, usually combining a direction search algorithm and a line search procedure along that direction. A well-established algorithm for the reliability index approach (RIA) is the modified HL-RF algorithm. For the Performance Measure Approach (PMA), the Hybrid Mean Value (HMV) scheme exhibits good performance. Although they are well established and widely used, LSA approximations such as FORM and SORM exhibit some limitations when the degree of nonlinearity in the LSF increases. Since the FORM approach is a linear approximation, it is easy to understand that it ignores all nonlinearities in the LSF. Furthermore, as the number of dimensions of the problem increases, also the error of the approximation for all LSAs increases accordingly, without the possibility to estimate it within the same reliability calculation (“curse of dimensionality”). In practice, this means that for problems with a number of dimensions n>15, special care should be paid to the degree of nonlinearity of the performance function(s) because the error of the approximation strongly depends on this nonlinearity.

4 ANALYSIS CASES IMPLEMENTATION In order to perform the analysis the methodology reported in Section 3.2 has been implemented in a MATLAB code and used to perform the analysis. The Optimization process is managed by an OPTIMUS thread (Figure 4

and upper part of Figure 6) that determines the necessary steps to reach the optimum design point that satisfies the probabilistic constraint of Eq. (16). The optimization algorithm selected for all the analysis cases is the Sequential Quadratic Programming (SQP) using a tolerance of 10-3 and forward finite difference gradient estimation. Note that the optimization process is a deterministic process with a probabilistic constraint. Thus the result is an optimum point that is also robust with respect to the required failure probability.

4.1 The composite wing For the composite wing structure proposed in Section 2.1, a Finite Element model has been created and the external load configuration has been applied. With the variability introduced in the input parameters, the performance behavior of the FE model has been observed and the RBDO process reported in Section 3 has been applied to meet the following probabilistic constraint:

βt = 6 ⇒ Pf = 9.86-10

(16)

Figure 4. OPTIMUS project for Reliability-based Design Optimization process

Figure 5. OPTIMUS project for Process automation For each iteration of the optimization process, the algorithm requires the evaluation of the performance function at specific values of the input parameters, so that a FE static analysis (for the present case the Nastran solver has been used) has to be performed for each evaluation. To accomplish this task, another process management thread (Figure 5) using OPTIMUS has been created: when the MATLAB algorithm (for the estimation of the performance measure) requires a performance function evaluation for a given parameter vector, OPTIMUS launches a FE computation and submits the results back to MATLAB, where the LSF criterion is evaluated.

4.2 The vehicle knuckle For the testcase proposed in Section 2.1, a similar procedure has been applied to optimize the fatigue life response of the knuckle FE model. Due to the computational effort required for each fatigue life evaluation, a hybrid meta-model/FE strategy has been used. In fact, a full factorial, multi-level Design of Experiments (DOE) plan has been designed in order to explore the fatigue life response within the boundaries of interest. Subsequently, a quadratic Response Surface (RS) model has been computed based on the DOE results. The obtained meta-model has then been used for a first RBDO process using a PMA methodology. The DOE strategy used to build the meta-model starts with a full factorial, 3 levels design, in the range [-3 , +3 ] . The DOE plan has then been expanded to a 5 levels design in the range [-6 , +6 ] and then up to a 7 levels design in the range [-9 , +9 ] . This choice has been made to better explore the possible nonlinear behavior of the fatigue life response and to validate the meta-model on a wider interval of analysis. Thus a total of 49 experiments have been run and the maximum absolute error of the meta-model showed to be less than 0.1 when compared with FE simulation results. The outline of the procedure used for the analysis is 1.

Compute the response of the structure with 3 DOE plans a. Compute the 3 levels DOE plan in the range [-3 +3 ] b. Compute the 5 levels DOE plan in the range [-6 +6 ] , re-using previous step results c. Compute the 7 levels DOE plan in the range [-9 +9 ] , re-using previous step results 2. Compute the quadratic least squares response model RS (using Taylor Approximation) 3. Perform an RBDO loop using the RS model and obtain a first tentative optimum point 4. Validate (and refine if needed) the results of the previous step with FE computations. The implementation of the analysis for the vehicle knuckle differs from the previous one in the process automation loop. In fact, in order to evaluate the performance function at specific values of the input parameters, a sequence of FE static analysis and LMS FALANCS fatigue analysis has to be performed for each evaluation. The probabilistic constraint set for the RBDO process is:

βt = 6 ⇒ Pf = 9.86-10

(17)

This means that the design optimum will be a 6 design w.r.t. the PF selected in Section 5. OPTIMUS Optimization Loop

MATLAB Algorithms (HL-RF, HMV)

OPTIMUS Process Integration

Figure 6: Optimization loop and Process Integration with MATLAB and OPTIMUS As in the previous implementation, the optimization process is managed by an OPTIMUS thread (as shown in Figure 6) that determines the necessary steps to reach the optimum design point that satisfies the probabilistic constraint of Eq. (17).

5 PERFORMANCE FUNCTION CHOICE The selection of a proper performance function is a key task in the reliability analysis process. Of course it is strongly case dependent and even for the same structure, it can change depending on the particular aspect that has to be analyzed. In this section, a description of the specific performance function selection for the two analysis cases is reported.

5.1 The Composite Wing As the wing considered in this paper consists entirely of composite material, a special failure criterion should be chosen. The majority of failure criteria for composite materials are based on energetic theories, but few of them, like the Tsai-Wu criterion, have been developed and tested for industrial composite applications. It is common convention to consider that the composite structure breaks when one of its constituents breaks, as until now there are no criteria that take into account the presence of residual stress in the laminate. In the case considered in this paper for the composite wing, the Tsai-Hill failure criterion, derived from the more general Hill criterion for anisotropic materials, has been used as performance function. This criterion is defined by the following formula:

2 2  σ  2 σ σ σ  σ   G = 1−  1  − 1 22 +  2  +  12    F1  F1  F2   F12   

5.2 Vehicle Knuckle (18)

Before discussing the choice of the performance function for the present case, some observations are needed. Various failure modes exist in durability cases, and the identiThe maximum Tsai-Hill index over all the plies of all fication of the failure mode is important in both testing laminates (i.e. a first-ply failure criterion) represents the and analysis. Testing at load levels that are too high will performance of the structure. An example of the resulting change the failure mode and give the incorrect safety multidimensional surface, given the input parameter varifactors or may identify incorrect critical locations. The ability, is represented in Figure 7: the displayed values failure modes here considered are belong to a section of the multidimensional space referA static (non-fatigue) failure mode: this results in ring to the main hyper-planes passing through the origin. • large structural deformation and is controlled To prevent having to deal with an infinite design space, primarily by the resistance of the net section. the transformed standard normal space is bounded in a • Low cycle fatigue: this is usually the result of a hyper-cube between –10 and +10 and the optimiza large plastic zone at a notch, or stress concentration. tion range for each variable is between –4 and Low cycle fatigue behavior depends on the notch +4 with respect to the nominal starting value (these severity and the inelastic material response. are the initial mean values of the parameters exhibiting High cycle fatigue: this is representative of a variability). As can be seen in Figure 7, the surface is • situation in which there is little plasticity. Here the quite flat around the nominal values of the design paranotch severity, manufacturing processes and residual meters and becomes quite steep (w.r.t. the scatter around stresses play an important role. the nominal value) when approaching the LSF. For fatigue analysis, the strain-life approach is typic10 x 10 E ally used in situations where yielding of the material may E 4 be present at certain locations on the structure. This ap5 LSF 2 proach has been developed based on knowledge of more 0 0 detailed material behavior when subjected to stress levels LSF above the yield strength of the material. On the contrary, -2 -5 -10 the stress-life approach is typically more suitable in situ0 -10 -5 0 -10 5 ations where the stresses are below (or not much above) 10 10 -10 -5 0 5 10 E the elastic limit of the material. (a) (b) Following the previous considerations, the choice of Figure 7. Contour plot of the Limit State Function in the performance function is related to the Maximum the transformed space Damage sustained by the structure. For this case, a Maximum Damage value of –1 (in log10 scale) has been selecFour analysis sub-cases have been designed, using either ted. The performance function can than be written as the input variability in Table 1 or the nominal values. The sub-cases follow the description given in Table 3. MaxDamage = −1⇒ In this paper, only material properties have been con⇒ G = −( MaxDamage + 1) ⇒ (19) sidered. Note however that the methodology is not lim⇒ G > 0 for MaxDamage < −1 ited to this class of parameters as similar analyses can be performed for geometrical or strength properties. The increasing number of parameters in the sub-case sequence With this choice, the G function is positive where in Table 3 will be used to show how the addition of new MaxDamage is less than –1, which is in agreement with parameters influences the results of the optimization pro- the definition of the reliability problem reported in Seccedure (and the corresponding reliability index for tion 3. each iteration). In this paper, only material properties have been considered. Note however that the methodology is not limited to this class of parameters as similar analyses can be performed with variability in geometrical or strength properties, load conditions, … . 8

1

E2

2

1

Testcase 1 2 3 4

Variable 1 Fiber E1

Variable 2 Fiber E2

Variable 3 Mean Fiber G12

Fiber E1

Fiber E2

Fiber G12

Fiber E1 Fiber E1

Fiber E2 Fiber E2

Fiber G12 Fiber G12

Variable 4 Mean (Rho) Mean (Rho) (Rho) (Rho)

Variable 5 Mean Stringer E1 Mean Stringer E1 Mean Stringer E1 Stringer E1

Testcase

Fiber E1

Fiber E2

Fiber G12

1 2 3 4

-1.0388 -0.8677 -0.7614 0.31719

3.9200 4.0134 4 -0.3487

Nominal -1.3472 1.7929 -0.03564

(Rho) Nominal Nominal 4 0.81411

Stringer E1 Nominal Nominal Nominal 1.7215

Table 5. Optimization results in standard normal space (Y-Space) PMA Figure 8. Optimization process for sub-case 1 using 6 using RESULTS PMA The results of the optimization processes for both the analysis cases are reported in this section. The differences in A visualization of the optimization process for subthe procedures and the influence on the results of the difcase 1 is shown in Figure 8(a) and Figure 8(b). ferent hypotheses are stressed. For the other sub-cases, the optimum points are reported in Table 5. It is shown that the position of the optimum point is noticeably affected by the introduction of new 6.1 Composite wing results random parameters. These random parameters are represFor the composite wing case, Table 4 shows the optimiza- entative of the variability of the input parameters and are tion results using PMA. The optimization process has described by probability distributions. As a result, while been repeated for all 4 sub-cases, in order to demonstrate in the first sub-case Fiber E1 is decreased and Fiber E2 is how the results may change by introducing additional increased, in sub-case 4 the opposite occurs. variability. An indication of the total computational effort spent for each optimization procedure is given by the required number of Limit State Function (LSF) evaluations: more computations means a higher computational effort. 10 8 6

2

2076

1775

-2

-2

Performance Measure G(x) SC 1 SC 2 SC 3 SC 4 0.2514 0.2523 0.2523 -0.044 0.1182 0.1252 0.1117 0.0046 0.0441 0.1242 0.1036 0.0025 -0.0034 0.0615 0.0695 -0.0009 0.0003 0.0163 0.021 -0.0005 0.0001 0.0054 -0.0108 0.0001 0.0049 0.0031 0.0046 0.0023 0.0033 630

0

0 G

SQP Iteration 0 1 2 3 4 5 6 7 8 N. of LSF evaluations

2

X

E2

4

-4

-6 10

-6

-8 0

-8 -10 -10

-4

-10 -10

-5

0 Stringer E1

5

-5

0

10

5

10 -10

E2

Stringer E1

(a) (b) Figure 9. Optimization process for sub-case 4 and parameters E2 and Stringer E1

3605

Table 4. Performance Measure estimation in the optimization results for different testcases using PMA (results in the table are expressed in terms of performance function value in the parameter space)

The introduction of the last parameter has a strong influence on the position of the optimum point. This is due to the fact that the stringers are stiffening elements in the structure, thus variability in their properties influences the structure more than the variability in the 2nd, 3rd and 4th parameter. This can clearly be observed in Figure 9 that plots the LSF for the 2nd and 5th parameter.

6.2 Vehicle knuckle results 10 8

x 10

8 4 6 3 4 2 GX

E1

2 0

0

-2

-1

-4

-2 -10

-6

-5

-8 -10 -10

1

0

-5 5

-5

0 E2

(a)

5

10

10 10

5

E1

(b)

0 E2

-10

Referring to the outline of the procedure reported in Section 4, different steps have been considered for the vehicle knuckle. Due to the higher computational effort required for each fatigue analysis, a design of experiments strategy has been used to assess the response of the structure with a limited number of evaluations . The first step consisted in the computation of the fatigue life (in terms of Maximum Damage) using 3 different DOE plans. Table 6 reports the values of the Maximum Damage for the 3, 5 and 7 levels DOE plans, as defined in Section 4.

Using the results of the DOE plan of Table 6, a least squares quadratic Taylor Response Surface has been computed to approximate the FE model. A plot of the response surface is shown in Figure 10. The RS has an absolute maximum error of less than 0.1, thus can be used as meta-model in order to assess the reliability of the FE model and to compute the first part of the optimization E ↓\ TS→ [Mpa] 182000 188000 194000 200000 206000 212000 218000

44 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

In order to meet both targets, the optimization procedure outlined in Section 4.2 has been implemented for an objective that is a combination of the two constraints of Eq. (20). The new objective can be expressed as:

) min( β S − βt − MaxDamage

560

680

800

920

1040

1160

-0.41494 -0.42428 -0.42778 -0.43051 -0.43406 -0.43653 -0.43881

-1.05 -1.0689 -1.0885 -1.1044 -1.12 -1.1369 -1.1511

-1.4369 -1.4634 -1.4879 -1.5122 -1.5342 -1.556 -1.5774

-1.6716 -1.7043 -1.7315 -1.7626 -1.7892 -1.8143 -1.8389

-1.8451 -1.8811 -1.9122 -1.9438 -1.9737 -2.0007 -2.0275

-1.9948 -2.033 -2.07 -2.1039 -2.1356 -2.1656 -2.1952

(21)

(in log10 scale) for the 7 levels DOE plan. Different colors show the 3 and 5 levels DOE plans. procedure. Using this new objective, the optimization process has been carried out using the Performance Measure Approach and the HMV algorithm , looking for a robust optimum point in the range [-3 +3 ] for each parameter (i.e. E and TS). The results of the optimization are reported in Table 7 and Table 8 and shown in Figure 11.

10

8

6

0.5 4

0 OPTIMIZED POINT 2 TS, Tensile Strength

Maximum Damage

-0.5

-1

0

STARTING POINT

-2

-1.5

-10

-2

-4

-5 -2.5 -10

-6

0 -8

-6

-4

-2

0

5 2

4

-8

6

8

10

10

E

Tensile Strength

Figure 10: Response surface of the Maximum damage (log10 scale) as a function of the Elastic modulus and Tensile Strength.

-10 -10

-8

-6

-4

-2

0 E, Elastic Modulus

2

4

6

8

10

Figure 11: RBDO procedure using the PMA and the Hybrid Mean Value (HMV) algorithm.

Apparently, the optimization procedure seems simple, as there’s only a single equality constraint to be formally satisfied (see Section 4). A closer look at the problem reveals that the optimization loop should in fact meet two optimization objectives: the first is the equality constraint and the second is that the newly found design point should also minimize the Maximum Damage. It’s easy to figure out that the solution locus for the first constraint is constituted by all the points at a distance of 6 from the LSF. But among these points, there is only one point that minimizes the Maximum Damage of the structure. Thus the multi-objective procedure should satisfy two constraints, one equality constraint and one inequality constraint: β t = 6  ) minimize( MaxDamage

(20)

SQP Iteration 1 2 3 4 5 6 7

Elastic Modulus (MPa) 200000.00 200155.19 200277.41 200415.38 201100.98 204519.58 206000.00

Tensile Strength (MPa) 800.00 897.67 905.78 905.75 905.38 903.51 902.79

Table 7: RBDO results using PMA

Number of SQP Iterations 7 Number of LSF evaluations 240 Time needed for completion 26 seconds

Table 8: RBDO performance using PMA The optimized point found with the use of the response surface model has then been used as a starting point for a refinement optimization procedure involving only FE computations. Before stepping into the refinement process results, some considerations on the choice of using a Response Surface approximation are reported, to point out the advantages and disadvantages of this meta-modeling technique.

6.3 Considerations on the Response Surface model To better understand the differences between the Response Surface model and the FE computations, some relevant contours of the Maximum Damage values in the stochastic transformed space are presented. These contours have been computed with the FE model and with the DOE/RS model. Figure 12 compares a quadratic RS model with a contour based on a higher number of FE analyses than the DOE plans points. In the region where the behavior is quite linear, the accuracy of the terms with order higher than linear might not be sufficient. This gives rise to slight errors in the estimates with the RS model in this linear region. It can however be concluded that the RS model is a suitable first approximation model that can yield sufficiently accurate results with a very limited computational effort. Moreover, it saves many FE computations in the optimization process, by locating an initial point which is very close to the solution, so that the subsequent refinement process (that usually implies additional FE computations) becomes much faster. 10

10 8

8 6

6

4

2

2

TS, Tensile Strength

TS, Tensile Strength

4

0

LIMIT STATE FUNCTION -2

0

-2

-4

-4

-6

-6

-8

-8

-10 -10

LIMIT STATE FUNCTION

-8

-6

-4

-2 0 2 E1, Elastic Modulus

4

6

8

10

-10 -10

-8

-6

-4

-2

0 2 E, Elastic Modulus

4

6

8

10

(a) FE Computation (b) DOE/RS Model Figure 12: Differences between the FE and the RS model.

6.4 Refinement of the RBDO results The results of the refinement process using only FE computations are here reported. As it can be seen, even if the number of required optimization loops is less than the number required with the RS model, the time needed for the complete process to finish is some orders of magnitude greater. The results of the refinement process are reported in Table 9 and Table 10 and the comparison with the RS RBDO process in Table 11.

SQP Iteration 1 2 3

Elastic Modulus (MPa) 206000.00 206000.00 206000.00

Tensile Strength (MPa) 902.58 894.75 894.23

Table 9: RBDO refinement results for the vehicle knuckle using PMA

Number of SQP Iterations Number of PF evaluations Time needed for completion

3 289 80920 seconds (= 22h 28m 40s)

Table 10: FE RBDO refinement loop performance for the vehicle knuckle case using PMA

RS RBDO Time (seconds) 26

FE refinement RBDO Time (seconds) 80920

Ratio (FE/RS) 3112.3

Table 11: Performance comparison of the refinement process with the RS model first optimization loop. The results obtained from the refinement process show that the RS model is indeed conservative. Figure 12 shows that the “real” LSF (in terms of FE computations) is located slightly further than the LSF of the RS model. This means that the optimum point found with the RS model is at a distance a little larger than 6 . However, the difference between the Tensile Strength (TS) of the RS model and the TS of the FE refinement is sufficiently small when compared to the standard deviation of the TS property. This means that the FE refinement process can be considered too expensive (in terms of computational time) w.r.t. the advantage it brings. To conclude, it must be noted that the limited number of significant digits available in FE computations slows down the convergence of both the PMA iterations and the SQP optimizations. In order to limit this effect, a semiadaptive step for the finite difference gradient estimation has been used.

7 CONCLUSIONS In this paper, a composite wing and a vehicle knuckle have been analyzed under different load conditions and with variabilities in the input parameters. Several observations have been made. It has been shown that the methodology used to optimize the presented structures enables the designer to reach a higher level of confidence in the reliability of the design performance in the presence of realistic variabilities. The presented procedures are generally applicable and not limited to a specific engineering field. Moreover, the simplicity of integration in the design loop has been demonstrated in Section 4. Also different analysis strategies have been outlined, to demonstrate the flexibility that these procedures have, as they take advantage of existing and well-known modeling techniques. Thus, the methodology outlined in this paper for the analysis of so different cases constitutes an effective general-purpose tool to assess the inherent risk of a design. It allows the design risk to be reduced in accordance with safety regulations, to gain advantages in costs and time-to-market. The major outcome of this analysis effort is a more robust and safe design. For the specific applications considered in this paper, some other considerations arise. In particular, the effect of nonlinearities introduced by new variabilities in the input parameters should not be underestimated. One should carefully investigate what are the most relevant paramet-

ers to be included in the optimization loop through a preliminary assessment of the performance of the structure. For this purpose, one can use PMA. As an example, the results of the analysis of the composite wing are not only influenced by the probability that high- or low-quality materials are used, but also by the probability that a particular combination of different panel and stringer material properties occurs. Furthermore, it has been demonstrated that the dimensionality of the problem plays an important role as well. The introduction of the 5th parameter variability drastically changes the configuration of the failure surface defined by the Tsai-Hill criterion. In the case of more computational intensive analyses, a hybrid meta-model/FE approach can be used that proved to be effective to save analysis time and get an accurate result as well. In fact, the response surface model used for the first RBDO process in the analysis of the vehicle knuckle exhibits some conservatism w.r.t. the FE computations of the DOE plan, but can be successfully used to find a sufficiently accurate result. Therefore, the refinement process (using FE computations in place of the RS model) yields a small improvement in the solution that had been found with the RS RBDO procedure. Thus, the FE refinement can be considered to be not necessary if the RS model captures the behavior of the real FE model with sufficient accuracy. In the majority of the cases, this applies for performance functions that don’t exhibit strongly non-linear behavior in the region where the response surface is defined.

8 ACKNOWLEDGMENTS The presented methodologies are being studied and developed in the context of the EC research training network (RTN) MADUSE (Modeling product vAriability and Data Uncertainty in Structural dynamics Engineering), see http://www.maduse.org/, and the ongoing research project ‘Analysis Leads Design – Frontloading Digital Functional Performance Engineering’, which is supported by IWT Vlaanderen. The support of the European Commission, of IWT and of the Marie Curie Fellowship Association is gratefully acknowledged.

9 REFERENCES [1] MIL-HDBK-17, [2] [3]

Polymer Matrix Composites, Volume 1, Guidelines, February 1992. Long, M.W.; Narciso, J.D.: Probabilistic Design Methodology for Composite Aircraft Structures, Final Report DOT/FAA/AR-99/2, June 1999. d’Ippolito, R.; Donders, S.; Tzannetakis, N.; Van der Auweraer, H., Vandepitte, D., Integration of Probabilistic Methodology in the Aerospace Design Optimization Process, 46th

[4] [5]

[6] [7]

AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Austin, Texas, USA, Apr. 18-21, 2005, AIAA-2005-2137. LMS International, LMS Virtual.Lab, Rev. 5A, July 2005. Sudret, B.; Der Kiureghian, A.: Stochastic Finite Element Methods and Reliability, A state of the art report, Report N. UCB/SEMM-2000/08 Department of Civil and Environmental Engineering, University of California, Berkley, November 2000. Melchers, R.B.: Structural Reliability Analysis and Prediction, 2nd Edition, 1999. Youn, B.D.; Choi, K.K.; Park, Y.H.: Hybrid Analysis Method for Reliability-Based Design Optimization, Journal of Mechanical Design, vol 125, pp 221-232,

June 2003, ASME. [8] Schuëller, G.I., Pradlwarter, H.J., Koutsourelakis, P.S., A Critical appraisal of reliability estimation procedures for high dimensions, Probabilistic Engineering Mechanics, Vol. 19, pp. 463-474, 2004. [9] The MathWorks Inc., MATLAB Version 6.5, June 2002. [10]Noesis Solutions, OPTIMUS Rev. 5.1, July 2005. [11]MSC, MSC/NASTRAN 2001, February 2001. [12]Barkey, M., Hack, M., Speckert, M., Zingsheim, F., FALANCS User Manual, Kaiserlautern, 1998. [13]d’Ippolito, R., Donders, S., Tzannetakis, N., Van de Peer, J., Van der Auweraer, H., Design Improvements of composite aerospace structures using reliability analysis, Proc. IMAC XXIII, Orlando, Florida, USA, Jan. 31-Feb. 3, 2005.