AIAA 2004-4618

10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference 30 August - 1 September 2004, Albany, New York

A Convex Hull Approach for the Reliability-Based Design Optimization of Transient Dynamic Problems Palaniappan Ramu*, Samy Missoum♣, Raphael T. Haftka§ *§

University of Florida, Gainesville, USA Ecole Supérieure d’Ingénieurs Léonard de Vinci, Paris La Défense, France

♣

Abstract Nonlinear problems such as transient dynamic problems exhibit structural responses that can be discontinuous due to numerous bifurcations. This hinders gradient-based or response surface-based optimization. This paper proposes a novel approach to split the design space into regions where the response is continuous. This makes traditional optimization viable. A convex hull approach is adopted to isolate the points corresponding to unwanted bifurcations in the design space. The proposed approach is applied to a tube impacting a rigid wall representing a transient dynamic problem. Since nonlinear behavior is highly sensitive to small variations in design, reliability-based design optimization is performed. The proposed method provides the designer an optimal design with a prescribed dynamic behavior. I. Introduction Currently, large-scale structural optimization problems can be solved efficiently with commercially available Finite Element software. However, these problems are often limited to linear and simple nonlinear behaviors. There are no systematic and efficient methods to perform optimization of highly nonlinear problems such as those encountered in transient dynamic problems (e.g., crash). They typically involve many difficulties such as: •

The computational expense associated with repetitive costly finite element analysis (Kurtaran et al., 2002)

•

The difficulty in computing the sensitivity, due to the presence of numerical noise and the often encountered convergence uncertainties in explicit dynamic codes.

In view of the aforementioned difficulties, researchers (Gu.,2001, Yang et al.,2001, Kurtaran et al., 2002, Sobieszczanski-sobieski et al., 2000) have used metamodels and design of experiments (DOE) to optimize (or simply improve) their design, *

Graduate Student, Department of Mechanical and Aerospace Engineering, [email protected] ♣ Assistant Professor, Departement Mécanique des Systèmes [email protected] § Distinguished Professor, Department of Mechanical and Aerospace Engineering, [email protected]

especially for crashworthiness. The metamodel (or response surface) is used to replace the responses extracted from expensive simulations by a simple analytical model. The model, which is a function of the design variables, is built based on sampling points selected from the design space with a chosen design of experiments. In addition, response surfaces allow the removal of the numerical noise due to the simulations. An important aspect of most nonlinear problems is that the structural behavior is extremely sensitive to small variations in design or imperfections. As an implication, some responses may be discontinuous due to the presence of bifurcations and limit points (Braibant et al., 2002, Missoum et al., 2004). This hampers the usage of gradient-based methods or response surface approximations to perform optimization. Traditionally, safety factors are used to account for uncertainties in deterministic optimization approaches. The use of safety factors assists in pushing the optimal design from the boundary of the infeasible domain. Hence, safety factors are a measure of safety in deterministic approaches. Alternatively, probabilistic approaches such as reliability-based design optimization can be used. They incorporate information about the uncertainties that typically appear in material properties, loading conditions, geometry and simulation. Therefore, they provide more accurate measures of safety. Another

1 American Institute of Aeronautics and Astronautics Copyright © 2004 by S. Missoum, P. Ramu and R.T. Haftka. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

advantage of this approach is that a particular reliability level can be specified. Due to its advantages, reliability-based design optimization (RBDO) has triggered a strong interest in the optimization community [Kharmanda et al., 2004, Royset et al., 2003, Koch et al., 2001, Youn et al., 2004] In this research, we investigate an approach to identify regions in the design space where the dynamic response is continuous. This is achieved by identifying clusters in the discontinuous response space. The cluster corresponding to “unwanted” behaviour can be identified. Each cluster translates into specific region in the design space. The boundaries between the different regions can be distinctly represented using explicit separation functions in terms of the design variables. This allows defining the boundaries of the failure domain. Then, response surface techniques can be used in the region of interest to perform optimization. In addition, as the regions in the design space are associated with various dynamic behaviours, the identification of specific regions allows the designer to specify a particular behaviour (e.g., no global buckling).

hull in defining the boundaries of the failure domain explicitly, the proposed method also includes the possibility of using response surface and gradientbased optimization and the avoidance of repeated costly finite element simulations. Section II discusses the discontinuities and clusters in the response for a simple shallow truss with geometric nonlinearity. Section III discusses the cluster formation in the response space, regions in the design space, and the definition of the separation functions (lines and convex hull). Reliability-based design optimization and response surface approximations are discussed in Section IV. Section V demonstrates the proposed method on a column impacting a rigid wall. Discussions on the results and the possible enhancements on the proposed method are presented in Section VI and concluding remarks are provided in Section VII. II. Discontinuities and clusters: a simple example As a demonstrative example exhibiting discontinuities, the classical symmetric shallow twobar truss (Figure 1) is used (Crisfield M.A., pp 3-7). The member cross-sectional area is A and the applied force is F. The displacement under the point of application of the force is u.

The explicit functions in terms of design variables dividing the regions in the design space, serve as constraint for an optimization problem or as a limit state for reliability-based design. Missoum et al., (2004) used straight lines and ellipse to define the boundary of the regions in the design space. However, it was shown that straight lines and ellipses might be too conservative for some problems. Indeed, among the set of points in a design space, a specific region of failure points is identified. The boundaries of this region when defined using straight lines or ellipse encloses all the failure points but it might also encompass many acceptable points. In order to overcome this difficulty, this work proposes the use of a convex hull to define the boundary of the region corresponding to failure points. A convex hull of failure points is the smallest convex set containing all these points. In the design space, the separation functions are the explicit boundaries that define the failure domain. In the context of reliability-based design optimization (RBDO), once these functions are defined, the failure domain is known explicitly and hence failure probability can be computed easily based on the randomness of the design variables. In addition to the advantages provided by the convex

F u

Figure 1. Two-bar truss The truss exhibits the typical snap-through behavior as depicted in the force / deflection (F(u)) diagram on Figure 2. Two curves are plotted corresponding to trusses with cross sectional areas A and A + ∆ A. It can be seen that for the same load F, the system

F

A A+∆A

u1

u2

Figure 2. Force – deflection diagram. Two-bar truss might exhibit snap-through (displacement u2) or not (displacement u1). This can happen for

2 American Institute of Aeronautics and Astronautics

infinitesimally small variation ∆A, hence the displacement u is a discontinuous function of A. This might be a serious limitation if optimal design is considered.

optimization and/or a limit state function if uncertainties are considered.

II. 2 Design of experiments and response study In order to study the behavior of the response, the design space needs to be sampled. The Latin Hypercube Sampling (LHS) technique (Wang., 2003, Butler., 2001) is used to achieve this. The ranges of the two variables are: A ∈ [70 , 150 mm] F ∈ [-700, -400 N ]

The design of experiments is constituted of 404 points. The couples (F, A) used are represented on Figure 3 with the axes scaled based on the maximum values Fmax and Amax. LHS is a space filling technique that generates points within the domain. However, there might be a lack of information on the boundaries of the sampling space. In order to enhance the quality of the sampling, the four vertices of the design space are also included in the design of experiments.

Figure 4. Displacement u with respect to force and area for the two bar truss problem. Two clusters corresponding to stable and buckling behaviors

Figure 5. Division of the two bar truss design space into two regions. The R region corresponds to a stable behavior

Figure 3. LHS design of experiments for the two bar truss problem. 404 points The corresponding values of the displacements, u, are plotted on Figure 4. As the displacement is discontinuous with respect to the force and the area, two clusters are generated. The cluster that contains the circled dots is associated with a snap-through behavior. Every point in each cluster corresponds to a point in the design space (Figure 3). The set of points corresponding to the two clusters is represented in Figure 5. The clusters result in two regions in the design space that can be separated, in this case, by a straight line. If we want to optimize the two-bar truss enforcing no snap-through, then the design space is limited to the R region. The function separating the two regions can be used as a constraint for

III. Identification of clusters and definition of separating functions III.1 Cluster identification In order to identify the regions of interest in the design space, the clusters created by the discontinuities in the response space need to be identified. Statistical methods are available to find clusters within a cloud of points. This work uses one of the most widely used methods, the K-means algorithm (Hartigan J.A., 1979). The basic idea of the method is to minimize the sum of the Euclidean distances of the points of a cluster to its centroid. The number of clusters to be identified is an input parameter. In the two-bar truss example, the cluster separation as depicted on Figure 4 is obvious. There are cases

3 American Institute of Aeronautics and Astronautics

where the separation is not as clear as will be seen in the other example treated in this paper. III.2 Separating functions Once the clusters are identified in the response space, the corresponding points in the design space can be isolated to form regions. In order to achieve this, simple geometric entities such as lines and convex hulls can be used that confine the points that belong to a cluster in a simply shaped domain. The construction of the separating functions is explained in the sequel with the following notation: Ns: Total number of sampling points S : Set of points with “unwanted” bifurcation points D: Maximum Euclidean distance in S

{

D = max Pi − Pj 2 ;(Pi ,Pj ) ∈ S2

}

A, B: two most distant points in S

( {

(A,B) = Arg max P i −Pj ;(Pi ,Pj ) ∈ S2 2

})

III.2.1 Linear separation function In a two dimensional space, the algorithm used for the linear separating functions is: Step 1: Find the maximum distance D and the corresponding points A and B. Step 2: Define the equation g(x, y) = 0 of the line (L) going through A and B. Step 3: Create two lines (L1) of equation gr1(x, y) = 0 and (L2) (gr2(x, y) = 0) with the same slope as (L) so that: ∀Pi ∈ S, g r1 (Pi ) ≤ 0 and g r2 (Pi ) ≥ 0 (1) An example of the definition of linear separation function in the design space (x, y) of an arbitrary problem is provided in Figure 6.

It is noteworthy that when the two linear functions are defined as constraints, they cannot be used at the same time since the feasible space they define is disjoint. Two optimization problems have to be created with gr1 and gr2 as constraint respectively. The solutions of these two problems must be compared to find the optimum. Missoum et al (2004) used simple geometric entities such as lines and ellipses to define the separation function. It is shown that the simplest separation function, the straight line, is suitable for fairly simple problems like the two bar truss but the ellipse is too conservative for problems where the clusters are not that distinct as in the two bar truss problem. This work extends the idea by using a convex hull to define the separation function. III.2.2 Convex hull separation function Several methods are available to construct the convex hull in the design space. This work uses the “convhulln” function available in Matlab to construct the convex hull. This function is based on Qhull which is a program to compute convex hulls in arbitrary dimensions in a provably stable way. Qhull implements the Quickhull algorithm (Barber et al., 1996) which starts with three vertices that are known to be in the final convex hull, and creates two faces with these three points (one facing in each direction). Then, recursively for each face, the farthest vertex from that face is determined. The three faces formed by this point and the three edges on the first face are created, and the algorithm repeats on all vertices outside of these new three faces. The key idea is that, given a triangle of three points of the original set, the points inside this triangle do not belong to the facet of the convex hull. Hence, they can be discarded. An example of the definition of convex hull separation function in the design space (x, y) of an arbitrary problem is provided in Figure 7.

gr1

gr2

Figure 6. Example of linear separation functions in a design space (x, y)

Figure 7: Example of convex hull separation function in a design space (x, y)

4 American Institute of Aeronautics and Astronautics

IV. Reliability-based design optimization (RBDO) Generally, RBDO problems are formulated as follows: Min : f(x) x

Such that: gd (x) ≤ 0 Pf =P(gr(x) ≤ 0 ) ≤Pftarget

(2)

with x={xd, xr} where xd and xr are vectors of deterministic and random design variables respectively. gd is a vector of deterministic constraints and gr is a limit state function. P(gr(x)≤0) is the probability of failure and Pftarget is the maximum allowed or target probability of failure. The limit state function (gr(x)=0) divides the design space into a failure domain (gr(x)≤0) and a safe domain (gr(x)>0) and hence serves as a safety criterion. IV.1 Failure probability computation Reliability-based design involves the computation of a failure probability as shown in Eq (2). The widely used methods to compute the failure probability are Monte Carlo Simulations (MCS) or moment-based methods such as the First Order Reliability Method (FORM). Here, we use MCS to estimate the failure probability. MCS involves the generation of sample points depending on the statistical distribution of the variables. The sample points that violate the safety criterion are considered as failed. The failure probability is computed as: num( g r ( xˆ ) ≤ 0) Pf = (3) N Where, ˆx is the randomly chosen sample point, num( g r (ˆx) ≤ 0) is the number of samples for which g r (ˆx) ≤ 0 and N is the total number of samples. Another widely used measure of safety is the reliability index. It can be directly computed from the failure probability using the following inverse relationship:

β = - Φ-1 (Pf )

(4)

Where, Φ is the standard normal cumulative distribution function and β is the reliability index. The convex hull presented in Figure 7 represents graphically the boundaries of the failure region. Each facet of the hull can be replaced by explicit equations in terms of the design variables. This allows representing the failure domain with explicit boundaries. Each equation is a limit state equation and a point is considered to be failed if it violates all

the limit state functions simultaneously. The notion of violation is defined by assigning an inequality to each facet of the convex hull. When a point belongs to the failure region, all the inequalities have a definite sign. In order to check the sign of each inequality for a point within the hull, the computation of each limit state function at a definite point in the convex hull is sufficient. As a particular point, the centroid of the convex hull can be chosen. For every realization during MCS, the sign of each inequality is compared to the sign of the same inequality evaluated at the centroid. If all the signs are identical, then the point belongs to the failure domain. This procedure is summarized in Figure 8. Convex hull p: number of facets

For i = 1 to p Find the straight line equations (gi) based on the vertices of the convex h ll

xc= Centroid of the convex hull For i =1 to p Find the sign of gi at = sgn (gi (xc))

For any point x x belongs to the convex hull if For j = 1 to p sgn (gi(x))=sgn(gi(x)) Figure 8. Procedure to find points belonging to the failure region defined by convex hull In the case where the failure domain is approximated by straight lines, as defined in Section III.2.1, the sign of the inequalities can be easily found based on the points A or B. IV.2 Response surface approximations MCS often generates numerical noise due to limited sample size. Numerical noise in failure probability or reliability index eventually leads the optimization to converge to a spurious optimum. The accuracy of the MCS is dependent on the number of samples. When the target failure probability is very low and only limited sample size is used to perform MCS, it approximates the failure probability to zero. This

5 American Institute of Aeronautics and Astronautics

hinders the sensitivity computations in gradient-based optimization. These difficulties motivate the use of response surface approximations which employ lower order polynomials to approximate the failure probability or safety index in terms of the design variables. These response surfaces are termed as design response surface and are widely used by researchers performing RBDO (Sues et al., 1996, Qu et al., 2004). In the design space, failure probability changes by several orders of magnitude. This steep variation in the failure probability introduces additional challenges in requiring to fit response surfaces of higher order. This increases the computational expense. Alternatively, response surface can be fit to the reliability index as it does not suffer abrupt changes in magnitude. This work uses response surfaces to approximate the reliability index. The widely used metrics to measure the accuracy of the response surface are the R square and the Relative Maximum Absolute Error (RMAE). R square is used to check the global accuracy while RMAE gives a measure of the maximum local error. The expressions for these measures are given as: Ns

R 2 = 1−

∑(y

i

− yˆ i ) 2

∑ (y

2

i =1 Ns

the model has also been parameterized. However, a fixed rectangular section of height 50 mm and width 40 mm has been used. The analysis is performed with the explicit software ANSYS/LS-DYNA and the

Figure 9. Tube impacting a rigid wall. simulation time is 40 ms. Four masses of 15 Kg are located at the four corners at the rear of the tube. The tube is meshed with 3600 reduced integration Belytschko-Tsai shell elements. The tube can deform in various ways after impact onto the rigid wall. Here, the dynamic behavior is divided into two main categories: crushing (deformation along its axis) and global buckling. Two examples of these behaviors are given in Figures 10 and 11.

(5) i

− y)

i =1

RMAE =

max( y1 − yˆ 1 , y 2 − yˆ 2 ,... y n − yˆ n ) STD

(6)

where yi is the actual response value and yˆ i is the corresponding predicted value. y is the mean of the actual response and STD is the standard deviation of the actual response values. Figure 10. Crushing of the tube V. Transient dynamic example The proposed approach to define the boundaries of the failure domain is applied to a transient dynamic problem. V.1 Problem description The problem considered is a tube impacting a rigid wall with a velocity of 15 m/s (Figure 9). The objective of this work is to optimally design the tube so that no global buckling appears. That is, a constraint on the dynamic behavior has to be defined to enforce a crushing of the tube following its axis. The thickness t and the length L of the tube are chosen as design variables. Note that the section of

Figure 11. Global buckling of the tube

6 American Institute of Aeronautics and Astronautics

V.2 Design of experiments LHS is used to sample the design space. The ranges of the two variables are: L ∈ [0.3m , 1.0 m] t ∈ [1.0 mm, 5.0 mm]

The design of experiments, constituted of 100 points, is depicted in Figure 12. As presented in section II.2, the four vertices of the domain are added to the design of experiments. Therefore, the total number of sampling points is 104.

Figure 12. LHS design of experiments for the transient dynamic problem. 104 points V.3 Deletion of invalid analysis The reduced integration Belytschko-Tsai shell elements exhibit spurious modes. The work done byartificial forces to overcome these modes is termed as hourglass energy. In order to study the response, FE simulations are performed for the 104 points from the DOE. For each experiment performed, the ratio of hourglass energy over the total energy is stored. If the ratio is higher than 10%, the analysis is considered failed. Out of the 104 experiments, 7 failed this criterion and were not considered. V.4 Response study and clusters To detect the design for which global buckling occurs, the absolute values of the maximum transverse displacements |Uxmax| and |Uymax| were stored.

use of the cluster identification technique (K-means) as described in Section III is a less arbitrary way of selecting the points. When using K-means with two clusters, the clusters identified are represented in Figure 13. The circled dots correspond to points with potential global buckling. The response and the clusters are also plotted in a 3D diagram with respect to the values of length and thickness (Figure 14). The clusters in the response space translate into corresponding sets of failure and acceptable points in the design space as represented in Figure 15.

Figure 13. Response projected on the (response, length) subspace. Transient dynamic problem

Figure 14. Response plot with respect to the length and thickness for the transient dynamic problem

The sum |Uxmax| + |Uymax| is used as a response that encompasses the buckling in the x and y directions. This quantity should clearly exhibit a discontinuity compared to the case when there is crushing of the tube following its axis. For a clearer visualization of the response behavior, it is projected on the (response, length) subspace (Figure 13). The points with the highest response value (i.e., sum of displacements) correspond to designs with global buckling. At this stage, one could impose an arbitrary limit on the response to select the points that “seem” without buckling. However, the

Figure 15. Distribution of failure and acceptable points in the design space for the transient dynamic problem

7 American Institute of Aeronautics and Astronautics

V.5 Separation function Separating functions can be constructed in the design space to define the boundaries of the failure domain. When using straight lines, the resulting separating functions are depicted in Figure 16. The two corresponding line equations are:  L   t  g r1 (t , L) = −  − 1.016  + 1.474 = 0 Lmax    tmax  (7)  L   t  g r2 (t , L) = −  − 1.016  + 0.949 = 0  Lmax   tmax  (8)

Min L,t

:V

Such that: E/ET ≥ 0.99 Pf =P(L,t ɛ Ωf ) ≤Pftarget

(9)

Where V is the volume, E is the internal (absorbed) energy and ET is the total energy. Ωf is the domain limited by the two lines as defined by equations (7) and (8) or the region enclosed by the convex hull as shown in Figure 17. L and t are design variables. While L is deterministic, t follows a normal distribution with a mean defined as the current iterate of the optimization process and a standard deviation of σ=0.06*tmax mm. The target failure probability is 1×10-3. The RBDO problem presented in (9) involves the computation of a failure probability to evaluate the probabilistic constraint. Since the boundary of the failure domain is available as straight lines or as facets of the convex hull, the failure probability can be computed. This section discusses the RBDO performed by using straight lines or a convex hull approximating the failure domain.

Figure 16. Linear separation functions used to define the failure domain for the transient dynamic problem Based on the convex hull approach, the boundaries of the failure domain are defined as presented in Figure 17.

Figure 17. Convex hull separation function used to define the failure domain for the transient dynamic problem V.6 RBDO problem The RBDO problem considered consist of finding the length L and the thickness t of the tube so that:

V.6.1 Failure region bounded by straight lines Here, we consider the failure domain delimited by two lines as depicted in Figure 16. There are two methods by which the failure probability can be computed for the line separation function. One method is to use MCS. The other method employs equation (10). When MCS is used to estimate the failure probability, the two lines are to be represented as inequalities. It can be observed that the feasible region is disjoint for each limit state. Any sample that violates both the inequalities is considered to be failed. Based on MCS, the failure probability is estimated using the relation in (3). To include the energy ratio as a constraint in the RBDO problem, it is approximated with a response surface that not only removes the numerical noise but also prevents the repetitive calls to costly transient dynamic simulations. Another response surface could be used to approximate the probability of failure which is also known to be very noisy. However, due to acute variations of the probability of failure, it is usually recommended to fit the reliability index instead, as mentioned in Section IV.2. Both response surfaces are fitted in the (L/Lmax, t/tmax) space with second order polynomials. The response surface for the energy ratio is presented

8 American Institute of Aeronautics and Astronautics

in Figure 18. The surface in the figure shows the quadratic surface and the blue dots represent the values of exact energy ratio. The fitted response may not capture the actual energy ratio exactly, but it follows the trend and represents the region of interest (E/ET ≥ 0.99) well. Similarly, for the reliability index, a second order polynomial is used and provides a rather good approximation of the exact reliability index.

Figure 19. Failure domain (thick line) when t is the random variable absorption requirement is met and the probabilistic constraint is active. The results are presented in Table 1.

Figure 18. Response surface of approximated energy ratio with respect to thickness and length for the transient dynamic problem. The approximated reliability index and the energy ratio are used in the optimization process. The results are presented in Table 1. The accuracy of the energy ratio and the reliability index response surfaces are given in Table 2 based on the error measures expressed in equations (5) and (6). Alternatively, failure probability can be computed using (10). Since L is deterministic, the probabilistic problem is one-dimensional, and the failure domain is a segment. This situation is represented in Figure 19 for a given value L0 of L and the scaled thicknesses t1 and t2 that limit the failure domain, the failure probability can be expressed as:

V.6.2 Failure region bounded by a convex hull In the case where the boundaries of the failure domain are approximated by the facets of the convex hull, MCS is used to estimate the failure probability. As discussed in section IV.1, the facets of the convex hull are replaced by inequalities and the sample points that violate all the inequalities are considered to be failed. Based on this, the failure probability is computed from (3). Both, the energy ratio and failure probability suffer from numerical noise. Hence, response surface approximations are used to perform optimization. The response surface for the energy ratio is the same as the one fitted for the line separation function case. For the reliability index, a second order polynomial in terms of L/Lmax and t/tmax are fit and is presented in Figure 20.

 t2 − t  t −t  − Φ 1  Pf= Φ (10)   σ   σ  where, Φ is the standard normal cumulative function.

Table 1. Optimal design for the transient dynamic problem – Line separation function Method for Pf computation MCS* Eqn (15) 100,000 samples

t/tmax 0.96 0.98

Optimum L/Lmax V(mm3) 0.677 521640.3 0.665 522024.7

Failure Probability 0.001 0.001

Based on the allowable failure probability Pftarget =1×10-3, the result obtained is such that the energy

Figure 20. Response surface for the reliability index The blue dots are the values of the actual reliability index. The approximated surface for the reliability

9 American Institute of Aeronautics and Astronautics

index seems to represents the actual value well. A more precise description of the accuracy is given by the error statistics in Table2. The approximated reliability index and the energy ratio can now be used to solve the optimization problem defined in (9). The results are presented in Table 3. Table 2: Error metrics for the response surfaces Metrics R2 RMAE

Separation function Straight lines Convex hull Straight lines Convex hull

Reliability index 0.9343 0.9104 0.7552 0.9920

Energy ratio 0.8992 0.7639

It can be clearly observed that the optimal design obtained using a convex hull presented in Table 3 is lighter than the design obtained using straight lines presented in Table 1. This is due to the fact that probabilistic constraint pushes the design away from the separation function (line) that is already somewhat conservative compared to the convex hull.

constraints are active. It can be noted that there are two local minima, and the results reported in Table 1 corresponds to the better of the two. VI. DISCUSSION In the transient dynamic example discussed in this paper, the convex hull approach seems to define the boundaries of the failure region more accurately than the two line separation functions. However, this is valid only around the optimum found. In fact, parts of the actual failure region might not be represented by the convex hull. The most stringent example of this phenomenon is given by the two bar truss. The convex hull defining the boundaries of the failure region for the two bar truss is presented in Figure 22. 130

11

Table 3: Optimal design for the transient dynamic problem – Convex hull separation function Failure Optimum probability* t(mm) L(mm) V(mm3) 4.74 563.8 430418.6 0.001 * Computed using MCS – 100,000 samples

In the course of optimization, it was observed that the optimization algorithm converged to different optimal design for different starting points. In an attempt to understand this, the optimization problem

Figure 21. Graphical representation of the transient dynamic optimization problem stated in (9) was treated graphically and is presented in Figure 21. The constraints are represented by dashed lines and the contours of the objective function is presented in the solid lines. The feasible domain is shaded. The optimal solution is such that the objective function is minimum and both

Figure 22. Two bar truss. Failure region approximated by a convex hull. In this case, the failure region is actually a half plane bounded by a single line as represented in Figure 5. However, it can be observed in Figure 22 that the convex hull leaves the region above the facet formed by endpoints 130 and 11 as safe whereas failure occurs in these regions. Moreover, the facets of the convex hull which are common to the boundaries of the sampling space are artificial. For the transient dynamic case, investigation of the region above the facet formed by endpoints (16, 30) of Figure 17 shows that the same problem occurs. The reason for this is that the convex hull is constructed based on the sampling points generated by the design of experiments only. As the set of sampling points is fixed, no other information is provided. Hence, the definition of the failure domain by the convex hull is very dependent on the initial design of experiments. In addition, the bounds on the design variables might dictate the geometry of the convex hull as no samples is generated outside those bounds. Hence, the area of the failure region may be highly underestimated. For example, inspection of Figure 22 shows that a

10 American Institute of Aeronautics and Astronautics

reduction of the area ratio to 0.4 would lead to a safe design. Clearly this is not possible as the design would buckle for this area. However, this problem does not occur for the transient dynamic problem as the length is deterministic and the failure region in the vicinity of the optimum is not limited by the bounds on the design variables. In order to minimize the impact of these difficulties and to improve the use of the convex hull approach, the following points can be considered: • A solution to reduce the error in the definition of the failure region would be to remove some of the inequalities defining the convex hull that are artefacts of the fixed design of experiments. Despite the fact that for two dimensional cases, the physics of the problem and a graphical inspection might help to decide which inequality should be dropped, the following procedure, extendable to higher dimensions, can be used to identify the artificial facets of the convex hull:

• For the facets of the convex hull that follow the boundaries of the sampling space (e.g., facet (3, 30) of the transient dynamic problem), the same procedure as described previously can be used. However, in this case, the test point will be generated outside the domain. Note that this type of facet would not be limiting the failure region artificially if the random variable distributions were forced to be within the boundaries of the sampling space. Hence, a solution to this difficulty would be to perform the optimization in a design space subset of the sampling space. While the design space would represent the means of the design variables, the sampling space would encompass the entire finite random variable distributions. An example of the design space being a subset of the sampling space is depicted in Figure 24. In the work presented in this paper, design and sampling spaces were identical.

1) For every point in the design space, evaluate all the inequalities of the convex hull 2) If there exist a facet such that no point satisfies it while violating all the other facets, then this particular facet is suspicious as it might be artificial. In the transient dynamic problem, it is the case for the facets with end points (16, 30) and (27, 25). Note that facets (3, 9), (9, 15) and (15, 27) also fall into this category. 3) A sample satisfying this condition can then be generated. For example, for a two-dimensional problem, it could be generated at the centre of the candidate facet at an eccentricity away from the convex hull (figure 23). Test point

Figure 23. Example of test point generated to check whether a suspicious inequality is to be removed. Transient dynamic problem. 4) This sample is tested for its response. If it fails, the facet is an artificial facet and the corresponding inequality can be removed from the approximation of the failure region. If the response corresponding to the generated sample exhibits acceptable behaviour, then the facet is a natural facet of the convex hull.

Figure 24. Example of design space subspace of a sampling space. The design space represents the values of the means of design variables. • Since the convex hull is based on the samples generated by the DOE, it is recommended to use a space filling technique in which the samples are spaced as uniformly as possible. For instance, the Optimal Latin Hypercube Sampling (OLHS) where the minimum distance between points is maximized can be used (Butler, 2001). However, this improves the distributions of sample points inside the domain and special care should be taken for the boundaries. In this paper, the vertices of the domain were added, but other points on the boundaries could be added to the design of experiments. VII. Conclusion A convex hull approach to handle discontinuous response in nonlinear structural problems is proposed. The method consists of identifying regions of the design space that encompasses points for which the response is discontinuous. The boundaries

11 American Institute of Aeronautics and Astronautics

of the failure region are defined explicitly based on a convex hull.

(DAC), Pittsburgh, PA, DETC2001/DAC-21063.

The explicit knowledge of the failure domain can be used for reliability-based optimization. The facets of the convex hull are transformed into explicit limit state inequalities that must be violated simultaneously for failure to occur. Based on Monte Carlo Simulations, this allows an efficient calculation of the failure probability. The proposed method is applied to the optimal design of a tube impacting a rigid wall.

Hartigan, J.A. and M. A. Wong (1979) Algorithm AS 136: A K-means clustering algorithm. Applied Statistics 28, 100-108.

The next stages of this research involve the improvement of the techniques used for the definition of the failure region for two and three-dimensional problems. Also, the methodology should be applied to other types of uncertainties such as material properties, loads and presence of defects.

Koch, P. N. and L. Gu (2001, 18-19 June). Addressing uncertainty using the Isight probabilistic design environment. In First Annual Probabilistic Methods Conference, Newport Beach, CA.

Acknowledgments The authors would like to acknowledge Dr. Noulis Pavlopoulos and Mr. Max de Grandi for their financial support in this collaborative research effort between the University of Florida and Pôle Universitaire Léonard de Vinci.

ASME,

Paper

No.

Kharmanda, G., A. Mohamed, and M. Lemaire (2002). Efficient reliability-based design optimization using a hybrid space with application to finite element analysis. Structural and Multidisciplinary Optimization 24, 233-245.

Kutaran, H., A. Eskandarian, D. Marzougui, and N. E. Bedewi (2002). Crashworthiness design optimization using successive response surface approximations. Computational Mechanics 29, 409421. Madsen, H. O., S. Krenk, and N. C. Lind (1986). Methods of structural safety, Prentice-Hall, Englewood Cliffs, NJ.

References Barber, C. B., D. P. Dobkin, and H. T. Huhdanpaa (1996). The quickhull algorithm for convex hulls. ACM Transactions on Mathematical Software 22(4), 469-483. Braibant, V., M. Bulik, M. Liefvendahl, S. Molinier, R. Stocki, and C. Wauquiez (2002, 25-27 November). Stochastic simulation of highly nonlinear dynamic systems using the M-Xplore extension of the RADIOSS software. In Proceedings of the Workshop on Optimal Design, Laboratoire De Mècanique Des Solides, Ecole Polytechnique 91128 Palaiseau. Butler, A. N. (2001). Optimal and orthogonal latin hypercube designs for computer experiments. Biometrika 88( 3), 847-857. Crisfield, M. A. (1994). Non-linear finite element analysis of solids and structures. Volume 1. John Wiley and Sons. Gu, L. (2001, 9-12 September). A comparison of polynomial based regression models in vehicle safety analysis. In ASME Design Engineering Technical Conferences – Design Automation Conference

Missoum, S., S. Benchaabane, and B. Sudret (2004, 19-22 April). Handling bifurcations in the optimal design of transient dynamic problems. In 45th AIAA/ASME/ASC/AHS/ASC Structures, Structural Dynamics & Materials Conference, Palm Springs, CA. Qu, X. and R. T. Haftka (2004). Reliability-based design optimization using probabilistic sufficiency factor. Paper accepted for publication by Journal of Structural and Multidisciplinary Optimization Royset, J.O., A. D. Kiureghian, and E. Polak (2003). Successive approximations for the solution of optimal design problems with probabilistic objective and constraints. In Der Kiureghian, Madanat and Pestana (Eds.), Applications of Statistics and Probability in Civil Engineering, pp.1049-1056, Millpress, Rotterdam, ISBN 90 5966 004 8. Sobieszczanski-Sobieski, J., S. Kodiyalam, and R. J. Yang (2001). Optimization of car body under constraints of noise, vibration, and harshness (NVH), and Crash. Structural and Multidisciplinary Optimization 22(4), 295-306. Sudret, B. and A. D. Kiureghian (2000). Stochastic finite element methods and reliability. A state-of-the-

12 American Institute of Aeronautics and Astronautics

art report. Report No. UCB/SEMM-2000/08. Department of Civil and Environmental Engineering. University of California, Berkeley. Sues, R. H., D. R. Oakley and G. S. Rhodes (1996). Portable parallel computing for multidisciplinary stochastic optimization of aeropropulsion components. Final Report, NASA Contract NAS327288 Wang, G. G. (2003). Adaptive response surface method using inherited latin hypercube design points. Transactions of the ASME, Journal of Mechanical Design 125, 210-220. Yang, R.J., N. Wang, C. H. Tho, and J. P. Bobineau (2001, 9-12 September). Metamodeling development for vehicle frontal impact simulation. In ASME Design Engineering Technical Conferences – Design Automation Conference (DAC), Pittsburgh, PA, ASME, Paper No. DETC2001/DAC-21012. Youn, B. D., K. K. Choi, R. J. Yang, and L. Gu (2004). Reliability-based design optimization for crashworthiness of vehicle side impact. Structural and Multidisciplinary Optimization 26, 272-283.

13 American Institute of Aeronautics and Astronautics