n ..,. .. . vi

 

at»?

:1.

f
. ids?
u»... . 3....1.

4!.
in
, .

t1?
.4... EU 5 F
L...” )4 . .fik 5

V II I. 71 ‘1}.
“wavfi‘flqwu .
t

‘ I
.3 f . v E: .
fixtkrO ! .1} ‘ an}

. .5 .
duff it tilt:

..
.ZJI azufih

if. flit

.. :7.

. .4.
2r
:5an
1.1
. .
in...
Va)

 

in}. [‘1.
hr}: EVIL.“ 3‘01.

 

(to (NM

 

LIBRARY
Michigan State
University

 

 

 

This is to certify that the

dissertation entitled

METAMODEL-BASED OPTIMIZATION APPLIED TO CRASH
WORTHINESS AND SHEET METAL HYDROFORMING

presented by

Aravindhan Ravisekar

has been accepted towards fulfillment
of the requirements for

MS . Mechanical Engr.
degree in

 

 

 

'Mafi professor

 

Date 28 HAM) ZODB

MS U i: an Affirmative Action/Equal Opportunity Institution 0-12771

PLACE IN RETURN Box to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 cJCIRCJDeteDuepes-pJS

 

METALMODEL-BASED OPTIMIZATION APPLIED TO CRASHWORTHINESS
AND SHEET METAL HYDROFORMING

By

Aravindhan Ravisekar

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
MASTER OF SCIENCE
Department of Mechanical Engineering

2003

ABSTRACT

METALMODEL-BASED OPTIMIZATION APPLIED TO CRASHWORTHINESS
AND SHEET METAL HYDROFORMING

By

Aravindhan Ravisekar

A metamodel-based optimization approach is presented that can be used efficiently to
replace complex, computationally expensive computer analysis and simulation codes in
optimization. Kriging is used as the statistical interpolation technique. Optimization is
performed using a gradient-free algorithm based on a Bayesian technique. The proposed
approach is applied to problems of crashworthiness optimization and sheet metal
hydroforming. The objective in the problem of crashworthiness optimization is to
minimize the sum of maximum section force and maximum displacement in the rails of a
bumper-rail assembly of an automobile. The objective in sheet metal hydroforming is to
arrive at an optimal fluid pressure-punch stroke path that would reduce wrinkling failures
and allow for greater draw depths before rupturing. The proposed approach is efficient in

terms of computational time involved in optimization.

To My Family

iii

ACKNOWLEDGMENTS

I would like to thank, first and foremost, my advisor, Dr. Alejandro Diaz, for his
guidance and support throughout my graduate career and during the completion of this
thesis. His leadership, support, attention to detail and hard work have set an example, that
I hope to match in my career. I am also grateful to the members of my committee, Dr.

Andre Benard and Dr. F arhang Pourboghrat , for their interest in this research.
I also owe a great debt of gratitude to Ms. Cori I gnatovitch for her continued support

throughout this thesis work. Thanks also to my family for supporting me in my

educational pursuits and all my friends for their encouragement.

iv

TABLE OF CONTENTS

LIST OF TABLES .................................................................................. vii
LIST OF FIGURES ................................................................................ viii
CHAPTER 1
INTRODUCTION ................................................................................... 1
1.1 Motivation ............................................................................... 1
1.2 Background ............................................................................. 2
1.3 Problem Statement ..................................................................... 3
1.4 Selection of Optimization Approach ................................................ 4
1.5 Thesis Outline .......................................................................... 5
CHAPTER 2
LITERATURE REVIEW ............................................................................. 6
2.1 Design of Experiment ................................................................. 6
2.1.1 Factorial Designs .......................................................... 7
2.1.2 Latin Hypercubes .......................................................... 8
2.2 Metamodel Selection .................................................................. 9
2.2.1 Least Squares ............................................................ 10
2.2.2 Interpolating and Smoothing Splines .................................. 11
2.2.3 Fuzzy Logic and Neural Networks .................................... 11
2.2.4 Radial Basis Functions and Wavelets ................................. 12
2.2.5 Response Surface Methodology ....................................... 13
2.2.6 Kriging ..................................................................... 15
CHAPTER 3
KRIGING AND EFFICIENT GLOBAL OPTIMIZATION .................................. 16
3.1 Kriging ................................................................................. 16
3.1.1 History and Terminology ................................................ 16
3.1.2 General Kriging Approach .............................................. 17
3.1.3 Derivation of the Prediction Formula ................................. 20
3.1.4 Evaluation of Kriging ................................................... 25
3.2 Efficient Global Optimization ...................................................... 30
3.2.1 Expected Improvement Criterion ....................................... 32
3.2.2 Evaluation of Efficient Global Optimization ......................... 35

CHAPTER 4

CRASHWORTHINESS OPTIMIZATION OF A BUMPER-RAIL ASSEMBLY ........ 41
4.1 Background on Crashworthiness Analysis ........................................ 41
4.2 Finite Element Model ................................................................ 43
4.3 Problem Statement .................................................................... 45
4.4 Optimization ........................................................................... 46
4.4.1 Case 1 ..................................................................... 46
4.4.1.1 Discussion of Results .......................................... 47
4.4.2 Case2 ...................................................................... 52
4.4.2.1 Discussion of Results .......................................... 54
4.5 Conclusion .............................................................................. 58
Chapter 5
FINDING OPTIMAL FLUID PRESSURE-PUNCH STROKE PATH IN
STAMP HYDROFORMING ...................................................................... 60
5.1 Background on Stamp Hydroforming ............................................. 60
5.2 Summary of Previous Works ........................................................ 62
5.3 The Finite Element Model .......................................................... 65
5.4 Problem Formulation for Surrogate-Based Optimization ....................... 67
5.5 Problem Statement ...................... . .............................................. 71
5.6 Optimization ........................................................................... 72
5.6.1 Step 1 ...................................................................... 72
5.6.1.1 Discussion of Results ...... ' ..................................... 74
5.6.2 Step 2 ...................................................................... 76
5.6.2.1 Discussion of Results ........................................... 78
5.6.3 Step 3 ...................................................................... 80
5.6.4 Step 4 ...................................................................... 81
5.6.4.1 Discussion of Results ............................................ 82
BIBLIOGRAPHY ................................................................................... 84

vi

LIST OF TABLES

Table 1. Initial designs for case 1 .................................................................. 47
Table 2. Search points generated by EGO for case 1 ............................................ 50
Table 3. Initial designs for case 2 .................................................................. 54
Table 4. Search points generated by EGO for case 2 ............................................. 57
Table 5. Initial Designs for step 1 .................................................................. 73
Table 6. Initial Designs for step 2 .................................................................. 77
Table 7. Search points generated by EGO for step 2 ............................................ 79
Table 8. Initial designs for step 3 ................................................................... 80
Table 9. Initial Designs for step 4 .................................................................. 81
Table 10. Search points generated by EGO for step 4 ........................................... 82

vii

LIST OF FIGURES

“ Images in this thesis are presented in color”

Figure l. Full factorial designs ....................................................................... 7
Figure 2. 31 Fractional factorial design .............................................................. 8
Figure 3. Latin hypercube sampling (LHS) ......................................................... 9
Figure 4. Latin Hypercube sampling of 21 points ................................................. 26
Figure 5 True and the kriging model ............................................................... 27
Figure 6 True Vs Predicted responses ............................................................. 28
Figure 7 Cross validation predictions from the kriging model .................................. 30
Figure 8 Surface plot of Branin function .......................................................... 37
Figure 9 Contour plot of Branin fimction ......................................................... 38
Figure 10 Initial LHS points for Branin function .................................................. 39
Figure 11 Initial points and points generated by EGO ........................................... 40
Figure 12 Top view of the bumper-rail assembly before crash ................................. 43
Figure 13 Symmetric model of the bumper-rail assembly ....................................... 44
Figure 14 Latin Hypercube Sampling for case 1 ................................................. 46
Figure 15 Scatter Plot of true responses before Optimization for casel ....................... 48
Figure 16 Scatter Plot of true responses after Optimization for easel ......................... 50
Figure 17 Resultant sectional force Vs. time at optimal solution for case 1 .................. 51
Figure 18 X-displacement Vs. time at optimal solution for case 1 ............................. 52
Figure 19 Latin hypercube sampling for case 2 ................................................... 53
Figure 20 Scatter Plot of true responses before Optimization for case2 ....................... 54

viii

Figure 21 Scatter Plot of true responses after Optimization for case2 ......................... 55

Figure 22 Resultant sectional force Vs. time at optimal solution for case 2 .................. 57
Figure 23 X-displacement Vs. time at optimal solution for case 2 ............................. 58
Figure 24 Optimum fluid pressure—punch stroke path from previous work .................. 63
Figure 25 FLD Failure and distribution of strains with the initial pressure profile ........... 64
Figure 26 Full model created for the stamp hydroforrning process with a hemispherical

punch, using square blank .............................................................. 65
Figure 27 Quarter-model for the hydroforrning process using a hemispherical punch with

a round blank ............................................................................. 66
Figure 28 F LD Failure and distribution of strains for high values of pressure ............... 68
Figure 29 Fluid pressure—punch stroke path: Upper limit ....................................... 69
Figure 30 F LD Failure and distribution of strains: Upper Limit ................................ 70
Figure 31 Latin hypercube sampling for step 1 .................. . ................................. 73
Figure 32 Scatter Plot of true responses before Optimization for step 1 ....................... 74
Figure 33 Scatter Plot of true responses after Optimization for step 1 ......................... 75
Figure 34 F LD Failure and distribution of strains after step] ................................... 76
Figure 35 Latin hypercube sampling for step 2 ................................................... 77
Figure 36 Scatter Plot of true responses before Optimization for step 2 ...................... 78
Figure 37 FLD Failure and distribution of strains after step 2 ................................. 79

Figure 38 FLD Failure and distribution of strains for the optimized pressure profile. . . . . ...82

ix

CHAPTER 1

INTRODUCTION

1.1 Motivation

Computer-based simulation and analysis is used extensively in engineering for a variety
of tasks. Practical engineering analysis oflen requires running complex, computationally
expensive c omputer analysis and simulation c odes s uch as finite e lement analysis and
computational fluid dynamics models. Despite the steady increase in computing power,
these analyses seem to pose computational problems in terms of time consumed. As a
result many optimization problems in industries become expensive when solved by

conventional optimization techniques.

In engineering, performance of a product is governed by its design. Computational
expense due to time could negate many good features of the product design. Hence in
recent years, industries have started focusing more on finding alternate techniques to
solve the computationally expensive optimization problems. A widely accepted practice
is to build an inexpensive approximation to the expensive problem and solve or optimize
the inexpensive approximation instead of the original problem. These approximations are

normally referred to as coarse models and the original problem as a fine model.

There are again several ways to build a coarse model. In terms of finite element analysis

a coarse model would be one with a coarser mesh than the fine model. This readily tells

that there is considerable loss of accuracy in a coarse model. Hence, there is a necessity

to develop a coarse model, which closely approximates the true problem.

1.2 Background

Initially, function approximations used Taylor series expansion, which required the
gradient of the function. Most of the complex problems, which need approximate models,
are normally multi-modal and highly nonlinear. This makes it difficult to get the
derivative information out of these functions. As a result, optimization methods that
require derivative information are not suitable for these problems. Lately, researchers
([28], [29]) in the field have come up with the result that gradient-based approaches are
not compatible for the complex problems of modern .day. Hence, derivative-free

optimization methods are gaining popularity among researchers.

Curve-fitting techniques such as least squares and response surfaces were initially used to
build the approximations. Another solution to the problem of fitting highly nonlinear data
involves the 11 se 0 f splines, w here the p olynomials are d efined in a piecewise m anner
rather than a single expression. Common to all the modeling techniques mentioned above
is the assumption that the underlying data are p olynomial in nature. A method, which
provided a suitable alternate to above techniques, is kriging [9]. The origin of kriging is
in mining and geostatistics. A Detailed description of these methods is given in chapter

two.

1.3 Problem Statement

The objective of design optimization is to select the values of design variables

xi(i =1,....,n) subjected to various constraints in such a way that an objective fimction

f = f (x) attains an extreme value. This can be expressed in the abbreviated form

min {f(x):h(x)=0,g(x)_<_0} (1.1)
x93"

with ‘R”the set of real numbers, f an objective function,xe SR"a vector of 11 design

variables, g(x) a vector of in inequality constraints, h(x) a vector of p equality

constraints, and X = {x 6 SR" :h(x) = 0, g(x) S 0} the feasible domain.

The field of interest in this thesis is simulation-based optimization. Consider, for
example, a crashworthiness test performed on a vehicle model. There may be hundreds of
design parameters and millions of degrees of freedom. The huge computational time
involved could hinder finding an optimum solution to the test problem. Replacing such
expensive test by computer simulation may reduce the cost of the experiment. But these
computer simulations may prove expensive in terms of computational time involved.
Hence, there is a need to find an inexpensive solution. The present work focuses on
building statistical approximations by taking data from such expensive computer
simulations. The statistical models are referred to as “metamodels”[35], meaning model
of a model. The so-called metarnodels are considerably less expensive because they

require less numbers of evaluations of complex finite element models.

The functions approximated by metamodels are typically (a) noisy or discontinuous (i.e.,
non-smooth) and/or (b) requiring a long time to compute. The terms approximations,
surrogate models, and metamodels will be used synonymously throughout this thesis. The
major disadvantage of replacing the original computer-based simulation model by a
surrogate model is that the solution of the metamodel optimization problem may not be
the solution to the true problem (i.e., we are actually solving a different problem).
Optimizing a surrogate model by a gradient-free optimization method requires more
function evaluations compared to gradient-based techniques. This thesis focuses on

finding an optimization approach that can address the problems of this genre.

1.4 Selection of Optimization Approach

The basic requirement for the optimization approach is the selection of the design points
for b uilding the s urrogate. T he (I esign o f e xperiments p lays a m ajor p art in s urrogate-
based optimization. The performance of the surrogate depends on where the points are
selected in the design domain. The most widely accepted method called Latin hypercube
sampling is used for designing the experiments in this thesis. The next step is to build the
surrogate based on true function evaluations at the selected design sites. The metamodel
in this thesis is built based on a statistical interpolation technique called kriging. A
derivative-free optimization technique called Efficient Global Optimization (EGO) is
used to optimize the surrogate. The compatibility of EGO and kriging will be discussed

later in this thesis.

1.5 Thesis outline

The rest of the chapters are arranged as follows. In chapter two, various available
methods for designing the experiments and modeling the surrogates are described, and
their relative advantages a re discussed. C hapter three discusses the k riging m etamodel
technique in detail and a test function is used to evaluate it. Also chapter three describes
the Efficient Global Optimization technique. In chapter four an example problem on
design for crashworthiness is solved using the suggested optimization approach. In
chapter five an example problem on Sheet metal hydroforrning is solved using the

suggested optimization approach.

CHAPTER 2

LITERATURE REVIEW

As discussed in introduction one can use approximations to enhance design optimization.
The surrogate based optimization approach is normally divided in to three stages. The
first stage is the selection of initial design sites or design of experiments. Some common
methods for designing the experiment are discussed in this chapter. The design of
experiment stage is followed by the metamodeling stage where a surrogate function is
built. This thesis uses kriging for metamodeling. Other metamodeling techniques, which

are widely used, are discussed.

2.1 Design of Experiment (DOE)

Building an approximation for computer simulations involves (a) choosing a strategy to
sample the region of interest and (b) constructing an approximation model to the sampled
data. The region of interest is referred to as “design space”. It is bounded by the upper
and lower bounds of each design (input) variable being studied. Design of experiment
techniques allows a d esigner to s elect the p oints intelligently and in such a w ay as to
produce an accurate and statistically meaningfiil approximation. This chapter compares
several experimental design types, which are widely used in engineering design

community.

2.1.1 Factorial Design

There are basically two types of factorial designs, full factorial designs and fractional
factorial designs. Both these types are characterized by the terms factor and level. A
factor is the dimension of the design space. A level is a sub-division of a factor. Each

factor (the dimension of the design space) is fixed at certain number of levels. Some

examples of full factorial designs are 2" and 3" designs. The power k in the designs

stand for the number of factors or dimensions and the numbers (2 and 3) represent the
number of levels. Consider a 22 design i.e., a two-dimensional space with two levels for

each dimension, the points are at the comer of a square. In a 32 design the points are at

the center of each factor apart from the ends. These two designs are shown in figure 1.

f 9

t + .
‘ O

22 Full factorial design 32 Full factorial design

 

   

 

 

 

 

Figure l. Full factorial designs

The size of a full factorial experiment increases exponentially with the number of factors;
this leads to an unmanageable number of experiments. Fractional factorial designs are

used when experiments are costly and many factors are required. A fractional factorial

design is a fraction of a full factorial design. 3‘1”” is a fraction of 3" design. If we take

the case of 32 design, one fractional design of it is 3'. A 31 fractional factorial design is

shown in figure 2.

 

 

 

 

 

 

Figure 2. 31 Fractional factorial design

2.1.2 Latin Hypercubes

Latin hypercubes were the first type of design proposed specifically for computer
experiments [36]. A Latin hypercube is a matrix of 11 rows and k columns where n is the
number of levels being examined and k is the number of design variables (dimensions). A
level refers to the total number of points to be sampled. Each of the k columns contains
the levels 1,2...n, randomly permuted, each level being equally likely. Then these k
columns are matched at random to form the Latin hypercube. Latin hypercubes offer
flexible sample sizes while ensuring stratified sampling. A 10 x 2 Latin hypercube

sampling is shown in figure 3.

 

 

 

 

 

Figure 3. Latin hypercube sampling (LHS)
As it can be seen from figure 3 there is no definite pattern for the selection of points. This
totally eliminates any bias in the selection of design sites. As said previously LHS is well
suited for computer experiments because selection of random numbers and random match
of columns is easily accomplished with computers. Due to the above discussed features

Latin hypercube sampling is used for designing experiments in this thesis.

2.2 Metamodel Selection

In this section, a brief overview of different surrogate modeling techniques is given. The
purpose is not to explain each in detail. Rather, the aim is to illustrate the wide variety of
approximations available. It should be made clear that the most important technique for
this thesis is kriging (Section 2.2.6). While kriging is presented in this overview, the

details of the technique are lefi for the next chapter.

2.2.1 Least Squares

This was the first kind of metamodel brought in to use. Least squares use polynomials.

For a given data, one can determine the regression coefficients ,6,- in, say, a quadratic

model of the form ([19])
we = 130 + Ax + flzxz, (2.1)

where x is the input variable and y(x) is the predicted value of the input. The regression
coefficients are found by minimizing the mean of the sum of squared errors (MSSE) of

the predicted output values at x,
1 n A 2
MSSE = ; 2U: - y(xr)) . (2 -2)
‘-1
Equation 2.1 can be extended to any k th- order polynomial as

j:(x) = ,30 + Ax + ,62x2 + ..... + ,kak , (2.3)
Equation 2.3 in matrix form is written as
Xfl = y (2-4)
Solving for the vector ,6 gives the least square estimate
[3 = (XTX)'1 X7 y. (2.5)
with the optimal values of the coefficients chosen, a simple polynomial equation relating

the inputs and outputs is built. This type of metamodel can be extended to any higher

order polynomial by changing the number of regression coefficients.

10

2.2.2 Interpolating and Smoothing Splines

Splines, where the polynomials are defined in a piecewise manner rather than as a single
expression for the entire data set, provide an alternate solution to the problem of fitting
highly nonlinear data. Simply put, several low order polynomials are fit to the data, each
over a separate range defined by the knots or break points. Then boundary conditions are

placed on the polynomials to ensure that the pieces match up with a prescribed order of
continuity. Most often, the cubic polynomial pieces with C 2 continuity are used. By using

C 2 continuity, the pieces have identical value, slope, and curvature at the knot locations.
In the case of interpolating splines, the knots are taken at the data points so that the spline
passes exactly through each data point. Therefore, while interpolating splines may be
more true to the sampled data at the data points than least squares polynomials (i.e., a

single polynomial), they can exhibit more waviness in between the data points [28].
There has been considerable work done on another regression model known as the
smoothing spline. By adjusting a weight factor known as the smoothing parameter, a

compromise can be reached between the goodness of fit seen in interpolating splines and

the smoothness of the approximating functions seen in least squares models.

2.2.3 Fuzzy Logic and Neural Networks

Even though splines do not require the fit and analyze methods that least squares models

use, they still make the assumption that the data can be locally fit by a cubic polynomial.

11

There are many other methods that avoid these assumptions by letting the data be fit in a
more free form. Two popular methods in use are known as fuzzy logic and neural
networks. In the former, variables are classified as belonging to sets that define an input
as say, high, medium or low. Rules are then defined to describe how the dependent
variable responds to the various input sets. By describing how closely an input belongs

to one set or another, an approximate response can be defined.

For neural networks, only the set of input and corresponding output values are required.
As the name suggests, the technique is analogous to the way the human brain learns. The
input values are passed to a layer of nodes, which alter the data using, transforms such as

the sigmoid or tansig functions

The outputs of the hidden layer nodes are then passed to the output layer where they are
recombined and scaled. Nonlinear optimization techniques are used to “train” the model
parameters of the nodes so that the network can predict the output or response for a given
input accurately. Having multiple nodes and layers allows for highly nonlinear models.
Because no prior knowledge of the behavior of the response is necessary, these methods

have gained in popularity throughout various fields.

2.2.4 Radial Basis Functions and Wavelets

Two other metamodels that have found wide use in image filtering and reconstruction are

radial basis functions and wavelets. The radial basis function metamodel was originally

12

developed by Hardy in 1971 as an interpolating scheme. Fifteen years later, the work of
Dyn made radial basis functions more practical by enabling them to smooth data. As the
name suggests, the form of these metamodels is a basis function dependent on the
Euclidean distance between the sampled data point and the point to be predicted. Radial
basis functions cover a wide range of metamodels, and one particular choice of basis

function leads to the well-known thin plate splines.

Wavelet metamodels are also relatively recent in their development. They are similar to
Fourier transforms, but wavelet transforms have the ability to maintain space or time
information. They are also much faster and require fewer coefficients to represent a

regular signal.
2.2.5 Response Surface Methodology

A response surface model may be written in general as follows

;=f(x1,x2,...,xk)+£ (2.6)
where 8 represents the total error (the difference between the actual values and the
predicted values) and k is the number of predictor variables in the model. The function f
is normally chosen to be a low order polynomial ([42]), typically second order. Such a
second order model with k variables may be represented mathematically as follows

A k A k k A

f=xio+ Zing-+2211] xixj (2.7)
i=1 i=1j=l

l3

for k variables the second order polynomial of equation (2.7) will have

p = (k + l)(k + 2)/ 2 parameters. The intercept or constant term 20 is always included in

the parameters of a model. The model may be expanded to include higher order terms.

Once a model is chosen as a response model, the next step is to obtain the estimates of
the unknown parameters. This step will be accomplished by using a multiple linear

regression analysis. In general the error terms are unknown and the regression method is

A
used to get the estimates of the parameters. A.- is such that the sum L of the square of the

error ICI'HIS

L = Z 8,-2 (2-8)

is minimized . The parameters 2.- can be obtained by solving the following equation

3 =(XTX)71XTy (2.9)
The polynomial model as determined above will usually include redundant parameters.
These redundant parameters will reduce the accuracy of the model. Several methods exist
for selecting the best subset of parameters to be used in the final response model. This
final model is referred to as reduced model. Stepwise regression is a procedure employed

to get the final response model. Interested readers can refer to ([23] and [42]).

14

2.2.6 Kriging
Another type of metamodel that is currently gaining popularity is kriging. In all of the

above metamodels, a fundamental assumption is that y(x) = f (x) +8 , where the
residuals a are independent, identically distributed normal random variables. The main

idea b ehind k riging is that the e rrors in the p redicted v alues, 8,. , are 11 ot independent.

Rather, the error is taken to be a systematic function of x. The kriging metamodel,

y(x) = f (x)+Z (x), is comprised of two parts: a polynomial f (x) and a functional
departure from that polynomial Z (x) . While the idea behind kriging is simply put here,

the details are left for the next chapter.

15

CHAPTER 3

KRIGING AND EFFICIENT GLOBAL OPTIMIZATION

3.1 Kriging

This section describes kriging in detail. The model fitting procedures are explained and
the prediction equation of the method is derived. The final part of the chapter evaluates

kriging.

3.1.1 History and Terminology

There has been much interest of late in a form of metamodeling known as kriging. The
name refers to a South African geologist, D. G. Krige who first used the statistics-based
technique in analyzing mining data in the 1950s and 1960s. His work was furthered in the
early 1970s by authors such as Matheron, and formed the foundation for an entire field of

study now known as geostatistics [6].

At the end of the 1980s, kriging was taken in a new direction as a group of statisticians
applied the techniques of kriging to deterministic computer experiments [25], [46]. This
new research avenue was quickly explored by a variety of authors as the potential for
more general engineering work became apparent. Their process of experimental design
and using kriging models was named Design and Analysis of Computer Experiments

(DACE), much as RSM refers to a set of tools for modeling and optimizing a system.

16

3.1.2 General Kriging Approach

This section gives a general overview of kriging. The majority of the metamodeling

A
techniques rely on the decomposition y(x) = y(x)+£,. where the residuals a are
independent, identically distributed normal random variables. The main idea behind

kriging is that the errors in the predicted values, a, , are not independent. Rather, they

take the view that the errors are a systematic function of x. To understand this statement
better consider the following illustration where a quadratic equation is fit via least
squares to a given data set. For a narrow region around two sample points, let the value of

true function be y(x) and its associated regression model be f (x). At the sample point
xi, we know the true function value is y(xi)andi therefore the error is
g(xi) = y(xi) — f (x,) . At a location xi + 6 , for some small 6 , we only know the value
of regression fimction f (xi +5). Kriging theory posits that if the error at x,- is large,
then there is good reason to believe that the error at xi + 6 is also large. This is because

there is a systematic error associated with the functional form f (x) . The closer a point x

is to x,- the more likely f (x) is to deviate from y(x)by the same amount as g(xi) .

A kriging metamodel in multidimensions, which takes the form

y(x) = f(x)+Z(X) (3-1)

17

is comprised of two parts: a polynomial f (x) and a functional departure from that

polynomial Z (x) . This can be written as [25]

A d
y(x) = Z fljfj (x) + Z(x) (3-2)
j=l

where f1 (x) are the basis functions (i.e., the polynomial terms), fl j are the

corresponding coefficients, and Z (x) is a Gaussian process. More precisely,

Z (x) represents uncertainty in the mean of y(x) with E (Z (x)) = O. The f (x) term is
similar to the polynomial model in a response surface, providing a “global” model of the
design space. In many cases f (x) is simply taken to be a constant term ,6 [25]. While
f (x) “globally” approximates the design space, Z (x) creates “localized” deviations so
that the kriging model interpolates the sampled data points. The covariance matrix of
Z (x) is given by:

C0v(Z(w),Z(x)) = 02R(w, x) (3.3)

Here 02 is a scale factor called process variance that can be tuned to the data and

R( w, x) is known as spatial correlation function (SCF).

The choice of SCF determines how the metamodel fits the data. There are many choices

of R(w, x) , to quantify how smoothly the function moves from point x to point w. one of

the most common [21] SCF’s (shown here for one-dimensional points w and x) used in

DACE is

R(w — x) = e(‘glw‘xlp ), (3.4)

18

where 6 > 0 and O s p S 2. The choice ofthis SCF is motivated by the fact that the
sample path produced is extremely smooth and also the function is infinitely

differentiable. As with all choices of SCF, the function goes to zero as lw— x| increases.

This shows that the influence of a sampled data point on the point to be predicted
becomes weaker as their separation distance increases. The magnitude of 6 dictates how
quickly that influence deteriorates. For large values of 9 , only the data points very near
each other are well correlated. For small values of 6 , points further away still influence
the p oint to b e p redicted b ecause they are still w ell c orrelated. The p arameter p is a
measure of smoothness. The response becomes increasingly smooth with increasing
value of p. One should also note the interpolating behavior of the SCF. As R(x, x) = 1 ,

the predictor will go exactly through any measured data point.

To extend the SCF to several dimensions, common practice in DACE is to multiply the
correlation functions. The so-called product correlation rule applied to Equation (3.4)

results in

d
-6d|wd -xd’p

R(w,x): [DIe (3.5)
d=l

where the subscript d=1 to D refers to the dimension. The metamodels are robust enough
to allow a different choice of SCF for each dimension. In equation (3.5) the values of (9
and p are different in each dimension. But some authors [25] suggest using the same
correlation parameters for all dimensions. In the work at hand (9 and p are fit for each

dimension.

l9

3.1.3 Derivation of prediction formula

Several quantities and formulae have to be defined before deriving the prediction
formula. The polynomial term of the model in equation (3.2) comprises of the d x] vector

of regression functions

f(X) = If] (X), f2 (x)9°°’fd (x)]t

and the associated vector of constants

fl = [filuBZa-"anlt-

We next define the n x d expanded design matrix

q

FftOCI)
ft(x2)

 

 

[ft (xn )_

Notice that each row is one of f (x) vector and there are 11 rows for the n sampled points.

If we notate the stochastic process as
2 =[Z(x1),Z(x2), ..... ,Z(x,, )]’
then for the output of the sampled data, y = [y1, y2,...,y,,]' , one can rewrite equation
(3.2) as
y = F ,6 + z (3.6)
Let R be the n x n correlation matrix with i ,j defined by R(x,,xj )as in equation (3.5) and

let

20

rx = [R(x1,x),R(x2,x), ..... ,R(x,,,x)]t
be the n x 1 vector of correlations between the point at which to predict ,x, and all the

previously sampled points xi.

The goal of this derivation is to determine what is called the Best Linear Unbiased
Predictor (BLUP). To do so, we will need a measure of the prediction error to quantify
the best covariance model and a constraint to ensure unbiasedness. This can be written as
a minimization problem using Lagrange parameters for the constraints. The term

unbiasedness refers to the fact that the expected value of the linear predictor must equal

the expected value of Equation (3.6). A linear predictor is one of the form cx’ y. That is,

it uses a linear combination of the sampled data points, y, with the linear weighting

vector cx. Thus, an unbiased linear predictor is one for which E (ext y): of F )6 , since

E (Z (x)) = 0 by definition. This also means that by similar reasoning, E(;(x)) = f x’ ,B.
Equating these two for all ,6 , we obtain the set of d unbiasedness constraints:

F’cx =fx (3.7)
we now show that for any linear predictor c,’ y of Y(x), the mean square error of
prediction is

E(c,’ y — Y(x))2 = E cx‘yy‘c, + Y2(x) — 2cx‘yY(x)J
= Ec.‘(F/3 + 2)(F/3 + z)’c. + (m + Z(x))2 — 2c.‘(Ffl + z)(fx’fl + Z(x)) J

The second equation is obtained by substituting for y and Y (x) in the first equation.

Expanding the above equation we get:

21

E cxtFfl(Ffl)tcx +Cxtzztcx +2cxtFflzcx +fxtflfltfx
Z(x)2 + 2 fx’ ,62(x) — 2cx’Fm’p — 2cx’zfx’ ,8 - ZcxtFflZ(x) — 2c,’zZ(x)

Substituting the unbiasedness constraint from equation (3.7) the above equation

simplifies to:
= Elcx’FmFm’ c. + cx’zz’c. + fx’flfl’fx + Z(x)2 — 2c.’1~m‘fl — 2c.’22<x)J
Z (x)2 = 022 , the process variance and 22 = 022R. The above equation simplifies to
= E[(c,‘F,3 — fx’fl)2 + cx’zz’cx + Z(x)2 —- 2c,’zZ(x)J
By again using the unbiasedness constraint, the first term in the above equation is zero.
= Elcx’zz’cx + Z(x)2 - 2c,’zZ(x)J
By using the identity E(x + y) = E (x) + E ( y) :
= cx’ch zz’ J+ ElZ(x)2 J— 2cx’E[zZ(x)]

E12002]: 0'22 , [5122’]: 022R and E[ZZ(x)] = azzrx
E(cx’y- Y(x))2 =cx’azzkc, + 0,2 —2c,‘a,2r, (3.8)

Introducing the Lagrange m ultipliers ll,- , i = 1 to k for the c onstraints, the b est linear

unbiased predictor is the one that minimizes the error above. The formal statement is

min lcx’azchx +022 —2cx’azzrxJ ‘41:th —fo (3.9)

ex

22

subject to: F’cx —fx=0
Taking the partial derivative of the objective function with respect to c]r yields
0122ch -— 0'2er — F). = 0
Notice that the factor of two was removed from the first two terms because defining a

new choice of Lagrange multipliers A. = 2,1 is permissible. From optimization theory,

we now have two sets of equations to solve in the unknown vectors ex and 1 :
Ftcx = f x
azchx — F11 = 0’2er

Rewriting this system of equations in matrix form yields:

0 Ft ‘1 fx
[F 0' 2Rj[ 6. j = [02271:] (3'10)

Solving this system we get the BLUP as

—1
y(x)=cx’y=[fx’ rx’] [0 Ft] [0] (3.11)

F R y

This can also be written as ([24])

A t A t -l A
y(X)=fx fl+r xR (y-Ffl) (3-12)
where
A 1 1 1
,B=(F’R’ F)“ F’R" y (3.13)
is the generalized least square estimator of ,6. Equation (3.12) is the Best Linear

Unbiased Predictor (BLUP) i.e., the prediction equation for kriging. Also it is shown in

literature ([24]) that the mean square error (MSE) for this prediction can be calculated as:

23

A2 _1
a (x)=0221-[fx’ rx’] [1: 1:] [f] (3.14)

This estimate gives a measure of how reliable the prediction is at a given location. The
MSE is defined as the square of the expected value of the difference between the

predictor and the true value.

With equations (3.12) and (3.14) one can predict the response and the MSE at any untried

location. However the question still remains of what values to use for the parameters

022 in equation (3.14) and t9 and p in SCF of equation (3.4). The DACE approach

applies maximum likelihood estimation (MLE). While it will not be proven here, the

MLE of ,B is the generalized least squares estimator shown in equation (3.13).

The MLE of 0'22 is:

 

A 2 ( F3)‘R"< F2)
or, = y y (3.15)
n
The log-likelihood of the stochastic process is:
-1 A 2
—2— nln(0'z )+ ln(det(R)) (3.16)

Equation (3.16) is the function of the unknown SCF parameters 6 and p only. While any
values for the parameters create an interpolative model, the “best” kriging model is found
by solving the d-dimensional unconstrained non-linear maximization problem given by

equation (3.16).

24

The value of p is often kept fixed at 2 [35] for all dimensions. The authors have
suggested that by just changing the values of 6 one can get a good prediction from
BLUP. Also it has been suggested by authors [35] that the polynomial approximation can

be reduced to a simple constant term ,6 without loss in model fidelity. This

simplification results in the following prediction formula

3m = fl+ r‘xR‘Ry—ffi) (3.17)

where f is an n x 1 vector of one’s and

[Ai=(f’R‘1f)‘1f’R’1y (3.18)

This thesis uses the above prediction formula given by equation (3.17).
3.1.4 Evaluation of kriging

This thesis is aimed at building approximations to functions that are expensive to
compute, for e xample f mite-element c odes. It is c onvenient, h owever to take e xample
functions from the optimization literature ([28]) to validate the kriging model. They

demonstrate many qualitative features of real fimctions. They are often highly non-linear.

This thesis uses one such example from literature:
f = 2+O.Ol(x2 —x12)2 + (1 —x1)2 + 2(2—x2)2 + 7sin(0.5x1)sin(0.7x1x2) (3.19)
the function is defined over the range x,- 6 [0,5]. Latin hypercube sampling was

employed for the initial data sampling. The number of initial points sampled was fixed at
21. This follows the rule of thumb given in the literature [31], ten times the number of

active design variables. The initial sample is shown in figure 4

25

 

 

 

 

X1

Figure 4. Latin Hypercube sampling of 21 points
Using this sampled data and the response of the complex function (equation 3.19)

corresponding to this data, approximation was built using kriging.

 

 

(a) True function

26

 

 

(b) kriging model

Figure 5 True and the kriging model

The three—dimensional surface plots of the true and the kriging models are shown in
figure 5. Figure 5a gives the surface plot of the true function and figure 5b gives the
surface plot of the kriging function. Graphically it can be seen that kriging does a decent
job in approximating this complex function with considerably less number of sample
points. The surface plot of the true function helps us visualize its multi—modality and non-
linearity. Though the above picture gives a qualitative evaluation of kriging, a model is
not completely valid without some quantitative results. The most straightforward method
of evaluating the above kriging model quantitatively is to use additional data points. 100
additional points were sampled using Latin hypercube. A simple measure of root mean

square error is given by the following formula

27

 

 

 

 

(3.20)

where y:- is the predicted value from the kriging model, yiis the true response and n, is
the number of extra points used to assess the error in the kriging model. The value for
the RMSE for the kriging model was 4.9944. Another method to visualize error is to
make a plot between the actual function values and the predicted values. A linear trend
line is used to fit the data. If the model closely represents the actual fiinction, then the

slope of the trend line would be close to 45°.

 

 

 

 

 

4O
35 1 ’
30 +
0
~ 0
'3 25 o
‘6 o
- 20
E 15 O .0 9. . O o.
n. o
to. °
or”. O '
10 O Q ..
5 O O O . O Q
9 o
0 ‘ . . ‘ . . . T .
O 5 10 15 20 25 30 35 40
true

Figure 6 True Vs Predicted responses

28

It can be seen from figure 6 that there are some outliers in the plot, but in an overall sense
the model seems to be a good approximation. Another technique called cross-validation
can be used to evaluate the kriging model. It is gaining popularity recently as it does not
require any additional sample points. The initial sampled data points are used to evaluate
the model. Cross—validation is used to assess the model accuracy by removing each
sample point used to fit the model one at a time from the data set, rebuilding the model
without the sample point and computing the difference between the predicted value at

that sample point and the actual value [20]. The cvnnse is then taken as:

 

 

 

 

cvrmse = (3.21)

where n3 is the number of sample points, y,- is the point currently not in the data set, and

A
yi is the predicted value yi. The value for the CVRMSE for the kriging model was

5.1645.

29

40
35 .
30
25 -
20 J
15 w

 

>

 

 

 

 

true res

Figure 7 Cross validation predictions from the kriging model

The line is extremely well fit proving that the cross-validated responses are close to the

actual responses.

3.2 Efficient Global Optimization

Once the metamodel is constructed, then an optimization algorithm can be employed to
optimize the surrogate. This section describes the Efficient Global Optimization (EGO)

technique, an optimization algorithm, used in this thesis.

Efficient Global Optimization uses a surrogate of the true function to generate the search
points that leads the algorithm to the solution of the problem. Since it uses a stochastic

model (kriging) for the surrogate, it falls in to the general class of optimization

3O

algorithms referred as Bayesian analysis algorithms. These global searching algorithms
use statistical models from an initial data sample to determine where to evaluate the

functions.

The simplest way to use surrogates for optimization is to create the metamodel and then
find the minimum (local search) of the metamodel. But the problem with this method is
that, it totally depends on the accuracy of the surrogate and can easily lead to a local
optimum. By simply finding the minimtun of the surrogate, the method does not
acknowledge the uncertainty associated with the metamodel. To eliminate this problem,
the method must put some emphasis on sampling (global search) where the uncertainty
associated with the metamodel is high. The uncertainty of the metamodel is measured by
the root mean squared error of the predictor. EGO uses Expected Improvement criterion

to strike a balance between the global and local search.

The EGO algorithm proceeds as follows:
1. Choose a small initial experimental design (set of points, say It) spread over the
entire input space and get the true function responses at these points.
2. Model (kriging) the true function using all the previous function evaluations.
3. Find the maximum of the “expected improvement” criterion. The location of
maximum is the new design point.
4. Evaluate the true function at the new design point.

5. Go to step two and continue till the algorithm converges to a solution.

31

The most important part of the algorithm is step three. The next design point found by the
algorithm is one with largest “expected improvement” in the current minimum. In simple
terms the next point for optimization [30] is one which has a response less than the

current minimum or is the point where the uncertainty is maximum.

3.2.1 Expected Improvement Criterion

Let us now explain the expected improvement criterion in detail. Let f {gm be the current

minimum function value, after the initial sample (n) evaluations. If the function is
sampled at x to determine y(x) then the improvement I over the current minimum is

defined as

I = {friiin — y(x) If y(x) < frflim (3.22)

0 otherwise

The expected improvement is then defined as [32]
fn’iin n
E(1)= I (fmin - y(X))¢(y(X))dy (323)

where ¢(y(x)) is the probability density function representing uncertainty about y. To

A
predict y at an untried location the kriging predictor y(x) (3.17) is used.

in) = 3+ r’xR“(y — 1,3) (3.24)

the mean squared error of the predictor is given by (3.14)

32

52(x)=0’22[1—[1 rx’] [(1) 11:] [H] (3.25)

More details on the kriging predictor and its derivation are available in section 3.1 of this

thesis. It is assumed [30] that y(x) has a Gaussian distribution with mean given by the
predictor, y (omitting the dependence on x for notational simplicity) and variance given

by 32(x). As a result ((y(x)-y)/s) is standard normal. By applying some tedious

integration by parts on the right hand side of equation (3.23), one can get an expression

for the expected improvement in closed form as shown below [11],

k n A n A if S>0
(min—m ffll +s¢ fmh’y

E(I) = 4 s (3.26)

 

[0 1'f s = 0
In the above equation (D( . ) is the cumulative density function. The first term in the
equation (3.26) is the difference between the current minimum and the predicted value

multiplied by the probability that the predicted value is smaller than the current

minimum. The second term is the product of the root mean squared error and the

probability that y is equal to f 13m . The second term is large when the predicted value is

closer to friiin and when s is large, i.e., when there is much uncertainty about the
predicted value. Thus the expected improvement will tend to be large at a point with
predicted value smaller than f n’iin and/or where there is much uncertainty associated

with the prediction.

33

The expected improvement criterion is then generalized [32] by introducing a parameter

gas,

0 otherwise

18 ={(f,:m —Y(x))g y(x)<fn’;in (3.27)

The greater the value of g, the more global the algorithm will tend to search. The

generalized form of equation (3.23) becomes,

fn'iin
E(Ig ) = I (fgm, - y(x))g ¢(y(X))dy

It has been shown that [32] the generalized expected improvement is given by,

 

 

A g-k
E(18)=s8 f (—1)"[ g! ] frgi" —y TK (3.28)
k=O k!(k — g)! s
where
A A k-l
Tk =—¢ mi: _y f'zi‘; _y +(k—1)Tk_2 (3.29)

The factor 3" in equation (3.28) gives the root mean squared error of the prediction, 8,
greater weight as g increases. Hence increasing the value of g shifts emphasis towards
global se arch. Although this statement is widely accepted in the literature, there is n 0
standard proof available for this. Too high of a value could prevent EGO from
converging on a good solution in a reasonable amount of time. Too low of a value could
lead EGO to overlook areas of high uncertainty as it searches too locally. It is proposed

[27] to start with a large g value, reducing the value towards zero.

34

The expected improvement criterion has to be maximized to find the next design point in

the search for the optimum solution. This thesis uses multiple starts of the simplex

algorithm for maximizing E(Ig ). At any existing design point the root mean squared

error is zero (the predictor is deterministic), and hence so is E (I g ). Thus the starting

points for simplex algorithm are taken from locations halfway between the existing
design points [32] in each dimension. Since the original design is space filling it is
ensured that the local starting points are spread over the entire design space. This does
not guarantee to find the global maximum, but since we are trying to find only the

location of next design point there is no need to find the global maximum.

Once the expected improvement is maximized and hence the location of next design
point determined, the kriging metamodel is updated with the information from the new
design point. Again the same procedure is repeated to find the next design point using
the information from the updated kriging model. The EGO algorithm is said to have

converged when the design point found is the same as one of the previous design point.

3.2.2 Evaluation of Efficient Global Optimization

This section evaluates the Efficient Global Optimization technique by using a two-

dimensional function called Branin function [4]. The function being moderately multi-

35

modal, would test the Approximating capability of the surrogate and the two—dimensional

form of the function would enable viewing the effects.

The analytical form of the Branin function is given below

y(x1,x7) : (x2 — 1x12 + 3x1 — 6)2 + 10(1 — i)Cos(x1)+10 (3.30)
" 47,2 n 87:

The graphical form of the function is given below

 

 

0 5 10 15 10

a) True function

 

 

O
0‘!

10 15

X
2

b) Kriging function

Figure 8 Surface plot of Branin function

The x ranges are —5 s x] s 10 and 0 3 x2 S 15. It has been found that ([5]), that there
exists three global minima for the Branin function and they are at points: (-rt, 12.25),
(9.4248225) and (rt, 2.25). This can be seen from figure 9. The metamodel is built using
the statistical technique, kriging. Initially the function is sampled at 21 points, generated
using Latin hypercube. The location of these 21 points can be seen in figure 10. The
choice of 21 is motivated by rule of thumb [30], “10 times the number of active
variables”. Considering 21 points instead of 20 conveniently spaces points at 5% of the

range.

37

 

 

 

 

 

Figure 9 Contour plot of Branin function

The true and the kriging metamodel are shown in figure 8. It can be seen that the kriging
approximates the function very closely. EGO proceeds by sequentially selecting one
point at a time. The 22“‘1 point is selected using the information from the initial kriging
metamodel. Once this point is selected, then the metamodel is updated and the updated
metamodel is used to select the next point in the optimization routine. The optimization

continues until the convergence.

38

 

 

 

 

 

Figure 10 Initial LHS points for Branin function

EGO algorithm took 12 additional points (figure 1 1), from the initial sampling to arrive at
a global minimum. The algorithm found the global minimum at (9.426,2.47) and this is
closest to the point (9.4248225). The value of this minimum is 0.397910. The algorithm
also found the other two global minima at (3.14375,2.26799) and (-3.14558,12.29152).

The true function values at these points are 0.398012 and 0.397938.

39

 

 

 

 

 

Figure 11 Initial points and points generated by EGO.

40

CHAPTER 4

CRASHWORTHINESS OPTIMIZATION OF A BUMPER-RAIL ASSEMBLY

In this chapter, the suggested optimization approach is used to solve an optimization
problem in crashworthiness analysis. The suggested optimization approach uses Latin
Hypercube for sampling the design points, kriging for building the metamodel and

Efficient Global Optimization for Optimizing the metamodel.

4.1 Background On Crashworthiness Analysis

As was discussed in chapter one, there are many complex optimization problems in
industries that require optimization through approximation. Optimizing a vehicle or a part
of it for crashworthiness is one such issue. The automotive industry always faces the
issue of maintaining vehicle safety standards. To meet these standards, and for studying
the behavior of a vehicle in crash, it is required to perform crash test on real vehicles.
This is an expensive problem in terms of time and cost involved. Lumped mass models
have been used since early 1970s [13] for analysis and design of automotive structures
for safety during crash. In these models, the major nonstructural components are
represented by lumped masses and the major deformable structural components are
modeled as nonlinear spring elements, typically represented with the force-deformation
data obtained from experiments. But this technique required an improvement in terms of
the energy absorption of the structure. With the advances in computing speed, finite

element models emerged as an alternate to lumped mass models. Finite element models

41

had a better performance in terms of representing the energy absorption of the actual

structure.

Yang et a1 [44] presented a feasibility study of using numerical optimization methods to
design structural components for crash. Both single and multiple objective formulations
were used for numerical Optimization, which resulted in better designs. It was found that
crash optimization was feasible but costly and that finite element mesh quality was

essential for successful crash analysis and optimization.

The commercial analysis packages require a detailed model of the entire vehicle. Such
models may take many hours of computing time. Due to the greater details needed in
building such models, optimization cannot be performed until the later stages in the
design process. To overcome this problem it was proposed [8] that the optimization
should take place much earlier in the design cycle. A lattice model was used to
approximate the b ehavior o f the continuum s tructure. It was shown b y S oto and D iaz
[39], that a lattice model could represent at least the relevant aspects of the complex

behavior of a structure, during an impact event.

Another approach to solve the complex crashworthiness problems is to build an
approximation of the complex finite element model and optimize the approximate model.
This thesis is developed based on this idea. Chapter two discusses some of the commonly

used approximation techniques.

42

4.2 Finite Element Model

A detailed multi-purpose finite element model of a 1994 Chevrolet C-1500 pick-up truck
was develOped at the National Crash Analysis Center (NCAC). The model is the first of
its kind developed specifically to address vehicle safety issues, including front and side
performance. It was decided to use the bumper-rail assembly from this model for this
thesis. The rail mounting which connected the bumper and the rails was replaced with

beam connectors.

 

1....

Figure 12 Top view of the bumper-rail assembly before crash

43

As the bumper-rail assembly is symmetric about the x-axis, the model was replaced by
the symmetric part. The rails were then divided in to two parts, the first part comprising
the portion of the rails in front of the rail connector and the second part comprising the
portion behind the rail connector. The rail connector is the green part shown in figure 12

that connects the two rails together. Figure 13 shows the model after the modifications.

 

Figure 13 Symmetric model of the bumper—rail assembly
The part in red is bumper, the part in blue is front—rail, the green part is the back—rail, the

orange part is the rail-mounting (beams) and the golden part is the rail connectors. A

mass of 100 pounds was added at the end of the rails. The bumper-rail assembly was
given an initial velocity 0 f20 mph. A planar rigid wall is defined on the nose ofthe
assembly so that the impact occurs right after the start of the analysis. The analysis is
given a termination time of 0.153 in a commercial crash analysis software. The total

computing time was 14 minutes on a Sun Ultra 80 machine.
4.3 Problem statement

The optimization problem that was proposed and solved in this thesis is defined below
Objective: To minimize the sum of maximum section force and maximum displacement

in the rails of a bumper-rail assembly.

find {151 1x2}
min F(x)

subject to: l 5 x1 5 70mpa
2.0 5 x2 5 5.0mm

where F (x)= (Maxdisp/200.0 + Maxforce/30000.0). The objective function is the sum

of maximum displacement in the rails and maximum section force in the rails. The
maximum force and the maximum displacement for each set of input variables are
obtained by using a commercial crash analysis software.

The design variables {x1, x2} are defined below:
1. x1= yield strength of the rail-mounting (beams)

2. x2 = sheet thickness of fiont-rails

45

 

4.4 Optimization

The next stage is to apply the optimization approach. The first step is to sample design
points in the design domain. Two different cases each with different number of initial
points were studied. Once the initial points were sampled, then the true fiinction was

evaluated at those points to get the actual response values.

4.4.1 Case]

The Latin hypercube sampling was used to sample 10 initial points and the model was

evaluated at those points to get corresponding response values.

5 t r r I I I

 

4.5 - —

x2 3.5 - ' -

2.5

I
O
l

2 41 r L L 1 L

0 10 20 30 40 50 60 70
x1

 

 

 

Figure 14 Latin Hypercube Sampling for case 1

46

 

 

 

 

 

 

 

 

 

 

 

 

 

9 x1 9‘2 yr

1 50.24 3.56 2.0349
2 63.18 3.31 2.0587
3 48.57 4.03 2.0128
4 9.27 2.20 2.2321
5 37.56 2.74 2.1544
6 26.26 2.47 2.2058
7 2.96 4.64 3.3142
8 18.04 4.72 2.0234
9 29.05 4.28 2.1701
10 63.02 2.92 2.1125

 

 

 

 

 

 

Table]. Initial designs for case 1

Kriging model was built with the sample points and the true responses at those sample

points.

4.4.1.1 Discussion of Results

Figure l 5 gives a s catter p lot 0 f t he response v alues against the d esign p oint n umber.
Though it may be statistically incorrect, it helps us to visualize the variation of the section
force. The minimum value of the response from the initial sample is 2.01275. The

corresponding values of the design variables are

47

 

Yield Strength of the rail-mounting beams = x1: 48.57

Sheet thickness of front-rails =x2 = 4.03

 

 

 

3.4 ' 7 "
O

3.2 - -

3 - -
8
9.

9 2.6 - -
3

*— 24 - -

4
2.2 - ' 5 3 9 -
2 ' . “b
l - 3 a
2 . o ..
1 8 l I I l 1 l l l l

 

'0 1 2 3 4 5 6 7 8 9 10
Design point number

Figure 15 Scatter Plot of true responses before Optimization for case]

This point is statistically important as it governs the initial search direction for the
optimization. Three possible cases are discussed depending on the location and
magnitude of this point.

Case]: If the current minimum lies in the immediate neighborhood of the global
minimum, then the o ptirnal v alues o f the d esign v ariables will b e closer to the design
values corresponding to the current minimum. Since, Efficient Global Optimization starts

with a point that is closer to the global minimum, it will take less number of searches to

48

 

arrive at the solution. Also the chances of optimization algorithm being struck in a local
minimum are avoided.

CaseZ: If the current minimum is closer to a local minimum and/ or is a local minimum,
then the optimization starts the search near a local minimum. Most of the optimization
routines have a fair chance of being mislead and ending up in a local minimum. But this
is taken care of in Efficient Global Optimization, because the next point in the search is
selected by not just focusing on improving the current minimum, but also on exploring
the unknown regions of the firnction.

Case3: The current minimum is actually the global minimmn. This is an ideal case. One
of the initial sample points has the optimal values for the design variables. In this case,
Efficient Global Optimization would try to search in the “little known” regions of the

function and would then converge to the global minimum.

The first search point in optimization and the subsequent points are selected based on the
Expected Improvement algorithm. As described earlier the points selected should have a
response value less than the current minimum and / or should be in the unexplored

regions of the design domain.

The EGO algorithm took 7 additional points to arrive at the “best possible” solution.

 

 

 

 

‘1 x1 3‘2 y.
11 70.00 5.00 2.0116
12 56.87 5.00 2.0185

 

 

 

 

 

 

49

 

 

 

 

 

 

 

 

 

 

13 16.75 2.00 2.4665
14 39.64 5.00 2.0872
15 70.00 2.59 2.1698
16 45.67 2.00 2.5818
17 70.00 5.00 2.0116

 

 

Table 2. Search points generated by EGO for case 1

 

3.4 r

True response
N N N 9°
.b O) (I) (p l\)

N
w
I

M
l
C
.m

 

F

1.8

.33-

5%

.V

L

 

 

6

8

Design point number

Figure 16 Scatter Plot of true responses after Optimization for case]

The EGO algorithm converged to a point with response value equal to 2.00116.The

optimum values for this case is presented as:

x1: 70.0mpa

2:2: 5.0011111]

50

 

35

 

o.)
a 0
d5—

 

 

 

N
O

 

\
l

...A
LII

 

 

 

Resultant Force(E+3)
8
“1

i
l \Viw \WWWMwfl

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14
Time

U)
1.x

 

 

 

 

 

 

 

 

 

 

 

 

Figure 17 Resultant sectional force Vs. time at optimal solution

The optimal value for yield strength of rail-mountings is 70mpa and the optimal sheet

thickness for the front-rails is equal to 5.0mm.

51

200

 

 

 

 

 

 

 

 

 

 

 

 

 

 

RE
RH

150 RE
.3 /
(u
E /
g 100
'5
.9
C.‘
x

50

0
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14
Time
Figure 18 X-displacement Vs. time at optimal solution
4.4.2 Case 2

In case 2 the number of initial sample points was increased to 20. This is in accordance
with the rule of thumb [31], ten times the number of active design variables. Latin

hypercube sampling of the 20 points and their corresponding response values are shown

in table 3.

52

2.5

 

 

 

 

r I 1 I ' I I
O
O
.. . _
O
O
- -i
O
O
.. . .4
O
O
' o
.. . d
O
l I l l l l.
0 10 20 30 40 50 60
X

1

Figure 19 Latin hypercube sampling for case 2

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1‘ x1 3‘2 yr

1 18.92 4.15 2.0834
2 7.57 3.58 2.4606
3 51.26 4.29 2.0099
4 10.08 3.74 2.3447
5 67.36 2.91 2.1152
6 47.74 2.39 2.2117
7 42.16 4.96 2.0649
8 33.21 4.69 2.1727
9 55.01 4.51 2.0135
10 61.00 2.11 2.4233
11 3.16 2.51 2.3075
12 24.11 3.30 2.1577

 

 

 

 

 

 

53

70

4.4.2.1 Discussion of Results

True response

 

 

 

 

 

 

 

 

 

 

 

 

13 29.97 3.47 2.1092
14 64.53 3.96 2.0164
15 37.76 3.19 2.0863
16 17.69 2.74 2.2207
17 43.68 4.79 2.0579
18 12.82 2.79 2.2518
19 58.25 3.93 2.0178
20 27.10 2.23 2.2577

 

 

Table 3. Initial designs for case 2

 

 

 

 

2.5 1 . .
3
10
o
2.4 - ~
4
o
1
2.3 - 3 -
1B 17'
o
5 W
2.2 ‘ B ‘
O 12
o
5 13

. . ' 15 _

2.1 3 7 . 17
' o
14 19
2 a l 9. I . l .
0 5 10 15 20

Design point number

Figure 20 Scatter Plot of true responses before Optimization for case2.

54

Figure 2 0 gives a good picture 0 f t he v ariation of the response v alues. The minimum
value of the response from the initial sample is 2.0099. The corresponding values of the
design variables are

Yield Strength of the rail-mountings =x1= 61.11mpa

Sheet thickness of front-rails =x2 =3.24mm

The EGO algorithm took 37 additional points to arrive at the optimal solution.

 

 

 

 

3 r l I l I
e
2.8- -
325» O -
8 o
0.
52m ' -
3 ' 0
e a
l‘22~ 0 '. ' ~
. ' o. o
e . O. O o
2~ e 0 O O 0% ~M~fi -
18 1 k r r I
0 10 20 30 40 50 60

Design point number
Figure 21 Scatter Plot of true responses after Optimization for case2
The optimal values of the design variables are
x1 =21 .21mpa

XZ =5.00mm

 

1‘ x1 3‘2 Yr

 

 

 

 

 

 

 

55

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

21 8.49 5.00 2.8684
22 70.00 5.00 2.0116
23 70.00 3.44 2.0475
24 57.18 5.00 2.0211
25 70.00 4.55 2.0076
26 37.14 4.09 2.0816
27 10.47 2.00 2.2478
28 47.26 3.59 2.0336
29 70.00 4.16 2.0116
30 50.57 5.00 2.0088
31 57.10 3.19 2.073
32 65.13 4.30 2.0080
33 54.89 4.21 2.0096
34 67.90 4.40 2.0067
35 67.23 4.74 2.0119
36 68.84 4.42 2.0093
37 61.39 4.30 2.0082
38 66.79 4.31 2.0075
39 63.24 4.31 2.0080
40 55.25 4.36 2.0105
41 59.45 4.29 2.0080
42 52.90 4.22 2.0103
43 67.60 4.37 2.0065
44 48.85 5.00 2.0026
45 48.08 5.00 2.0023
46 48.42 5.00 2.0022
47 24.76 3.99 2.1543
48 63.46 5.00 2.0157
49 52.44 3.62 2.0301

 

 

 

 

56

 

 

 

 

 

 

 

 

 

 

50 61.87 4.58 2.0110
51 48.60 4.95 2.0033
52 37.84 2.00 2.6138
53 21.58 5.00 1.9939
54 20.24 5.00 1.9987
55 21.21 5.00 1.9888
56 21.19 5.00 2.0066
57 21.21 5.00 1.9888

 

 

 

 

 

Table 4. Search points generated by EGO for case 2

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

35
3. i.
n 2-
m
+
Hi
01 2
8
o
u.
E 15
B
a 10
a)
“ W]
5 P1 A A A Awfi w v _
V\, U V \4 U V‘\/
0
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14
Time

Figure 22 Resultant sectional force Vs. time at optimal solution for case 2

57

 

The EGO algorithm converged to a solution which is less than the Optimum found in case
1. From this it can be concluded that the optimum found in case] is not the actual
solution. This helps in inferring that the accuracy of kriging model depends on the
number of initial sample points used. A plot Of the resultant sectional force against the

time, for the optimal values of the design variables is shown in figure 22.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

200
/Rma
150 ha.
/
E
(D
E /
8 100'
to
a
g
C.‘
X
50
o
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14

Time

Figure 23 X-displacement Vs. time at Optimal solution for case 2

4.5 Conclusion

After analyzing the two cases it can be concluded that the minimum value Of Objective

function is] .9888. Also it can be inferred that the model built with 20 initial design points

58

is more accurate than the model with 10 initial sample points. This is because kriging is
an interpolative model. Hence the larger number Of points the more accurate is the model.
The work in this thesis has given an important lead to the application of surrogate -based
optimization in crashworthiness analyses. It has been proved by the results from this
chapter, that, statistical models can be used as surrogates while modeling complex
problems in the design for crash worthiness. It can be observed that the optimum from
both the cases are at the boundaries of the domain. This could be avoided by using

different weights in the Objective function.

59

CHAPTER 5
FINDING OPTIMAL FLUID PRESSURE-PUNCH STROKE PATH IN SHEET
METAL HYDROFORMING
In this chapter the suggested optimization approach is applied to a sheet metal
hydroforming process to arrive at an optimal fluid pressure-punch stroke path. This

would reduce wrinkling in the finished part.

5.1 Background on sheet metal hydroforming

In sheet metal stamping process a punch and a blank holding die are used to plastically
deform a blank sheet of metal in to desired shape. This technique is being widely used to
fabricate commercially important structures, 6. g. automotive, aerospace etc... Depending
upon the geometry, volume and material type, different methods are used to form sheet
metals. The sheet metal forming encounters several types of failures, such as rupturing,

necking, wrinkling and springback that are undesirable.

Sheet metal stamp hydroforming is currently considered as a viable alternative to
conventional sheet metal stamping. In the fabrication Of complex structures, stamp
hydroforming exhibits many advantages compared to conventional stamping. These
advantages include improved forrnability due to the applied fluid pressure, low wear rate
of tools, better distribution of plastic deformation, and finally significant economic

savings due to less tooling [45]. In stamp hydroforming, one or more of the sheet metals

60

are supported with a pressurized viscous fluid to assist with the stamping of the part and a

female die is not required. The advantage of using pressurized fluid is discussed below:

(1) It supports the sheet from the start to the end Of the forming process, thus
yielding a better finished part.

(2) It delays the onset of material failure.

(3) It h as the p otential to reduce wrinkle formation when applied to o ne 0 r

both sides Of the sheet metal.

It had been demonstrated [18] that the use of hydrostatic pressure prevents the initiation
and spreading Of microcracks within the metallic material. Based on the success in using
a hydrostatic pressure to delay the onset of fracture, the idea of stamp hydroforming was
investigated as an alternative for conventional stamping. The sheet metal hydroforming
investigated in this thesis is a process in which a part is formed by a simple hemispherical
punch. The basic setup consists of a clamping mechanism and a fluid chamber. The fluid
chamber forms the upper half of the die. But in this case, the die half doesn’t take the
shape of the female part. The fluid chamber is lowered and the work piece is clamped
securely between the two die halves, thereby creating a seal for the upper fluid chamber.
The fluid is injected in to the upper chamber and given an initial pressure. As the punch
moves up, the sheet metal begins to deform and takes the shape of the punch, in this case
a hemispherical shape. The pressure is controlled via a pressure transducer and is used as

a means of forcing the sheet to conform to the shape of the punch. Once the punch

61

reaches the prescribed depth, the fluid is drained and the chamber is raised to remove the

formed part fi'om the die.

5.2 Summary of Previous works

The challenge in using sheet hydroforming process lies in finding an appropriate fluid
pressure-punch stroke path, which will avoid rupturing or wrinkling in the formed part.
Tirosh and Yossifon (1997-1998) [41] have performed many simple analyses of the
hydroforming deep drawing process as applied to the formation of cups from metallic
materials. The goal of their analyses was to establish an optimmn fluid pressure path
relative to the punch stroke that would prevent part failure due to rupture or wrinkling.
The findings from their studies indicated that e xcessive and insufficient fluid p ressure
lead to premature failure of hydroformed parts. The purpose was to arrive at a
predetermined path for the fluid pressure relative to the punch stroke that can be followed

to produce parts that are defects free.

Through their work they were able to show that rupture instabilities occur when the fluid
pressure was too high. The fluid pressure constrained the motion of the part and forced
the punch through the material. LO, Hsu and Wilson (1993) [17] expanded upon the
earlier work of Tirosh by applying the deep drawing hydroforming theory to the analysis
of hemispherical punch hydroforming process. These were followed by much of research

through numerical simulations and experimental techniques to arrive at the Optimal path.

62

The work in this thesis is aimed at finding Optimal path for the fluid pressure, taking lead
from the work done by Abedrabbo [1]. A rigorous finite element model was developed in
his work for better understanding of the deformation of the sheet metal during the
forming process. This model served as a tool to verify the experimental results. An
Optimal fluid pressure- punch stroke path was developed using the finite element model

and is shown in figure 24.

 

 

 

 

Optimal Path
0.012
o
. o 0
A 0.01 . o o
N o
a 0.008 « .
o
g 0.006 .
o
g 0.004 4 .
tL
0.002 - I
.0
0 4M . . . .
0 5 10 15 20 25 30 35
time (ms)

Figure 24 Optimum fluid pressure—punch stroke path from previous work

The optimal pressure path was categorized as one that eliminated wrinkling as seen from
the finite element model. The forming limit diagram for this model is shown in figure 25.
In figure 25 the tearing limit is given by red line in the F LD and the wrinkling limit is the

green line.

63

 

 

 

 

 

 

 

 

 

 

 

Major Engineering Strain(%)

 

 

 

 

 

 

 

 

Minor Engineering Strain(%)

Figure 25 FLD Failure and distribution of strains with the initial pressure profile

Each point in the FLD represents the strain state of a particular element in the finite
element model. It can be observed that a large number of elements fall below the
wrinkling limit, With the model being very safe regarding the failure due to tearing. First
let us discuss some of the drawbacks of the procedure described above. The method is too
expensive in terms of the computational time involved. Also the procedure failed to
identify the pattern that governs the behavior of the sheet metal, with the changes in the
fluid pressure. The total number of elements in this finite element model, that fall below
the wrinkling limit is 163. The proposed optimization approach is aimed at improving

this number.

It was then decided to use the optimum pressure profile developed as the reference for the
present work and use the surrogate—based Optimization to arrive at a comparatively better

pressure value for each time step, in the context of reducing wrinkling.

5.3 Finite element model

The finite element model of the punch and die is used as the complex model for the

surrogate based optimization.

 

Figure 26 Full model created for the stamp hydrofomiing process with a hemispherical
punch, using square blank

65

Two types of finite element models were available, a full model and a quarter model and

are shown in figure 26 and figure 27.

 

Figure 27 Quarter-model for the hydroforming process using a hemispherical punch with
a round blank

The full size finite element model used approximately 10500 four- noded shell elements,
while the quarter model used 4600 elements. The punch, die and the blank holder were

created using rigid material.

66

5.4 Problem formulation for surrogate-based optimization

This section describes the technique adopted in applying the surrogate-based Optimization
method for Optimizing the fluid pressure-punch stroke path in stamp hydroforming. The
objective for surrogate-based Optimization was selected as minimizing the number Of
elements of the finite element model that stay below the wrinkling limit subjected to the
constraint that, strains in all elements should stay below the tearing limit. The fluid
pressure is selected as the design variable for the Optimization. A more detailed

explanation is given in the next paragraph.

The termination time for the stamping process was observed to be equal to 27ms i.e., for
the given initial condition and the required punch depth, the process takes 27ms to
complete the forming. It was then proposed to divide the total time into ten stages,
starting from Oms and ending at 27ms with each stage covering a time span of 3ms. The
aim is to find an Optimum pressure value for each time step. It was decided to input
pressure values for two subsequent stages at once, and analyze the results from the finite
element model. This makes the Optimization problem two dimensional, since two design
variables are u sed at input. It s hould b e n oted that, increasing the n umber o f p ressure
stages at input increases the dimensionality of the Optimization problem. The total

number of stages devised initially restricts the maximum dimension of the problem.

67

A detailed study of the hydroforming process was conducted to obtain more insight in to
the process. High values of pressures were chosen randomly using the supplied pressure
profile as the datum and the analysis was preformed. The Forming Limit Diagram (FLD)
from this analysis- at time equal to 27ms is shown in figure 28. In figure 28 the tearing

limit is given by red line in the FLD and the wrinkling limit is marked by green line.

120

 

 

 

 

“—4

 

 

 

 

 

 

 

Major Engineering Strain(%)

i

 

 

 

 

 

 

 

0 , ’ ’

 

-i o 0
Minor Engineering Strain(%) .

Figure 28 FLD Failure and distribution of strains for high values of pressure

From figure 28 it is observed that even when large number of elements have crossed the

tearing limit. there are still sizeable number of elements that fall below the wrinkling

68

limit. From this I concluded that it is impossible to achieve a totally wrinkle free forming

for this material using the hemispherical punch.

The constraint in the problem is satisfied by suitable selection of upper bounds on the
design variables. After studying the effects from different pressure profiles the upper
bound on the pressure profile was obtained. The limiting fluid pressure-punch stroke path

Obtained is shown in figure 29.

Pressure Profile: Upper limit

0.012
001: ° 0 o o

0.008 1 ,

 

Pressure (GPa)
.0 9
§ §
0

g:
§
O

 

 

it

 

1 0 15 20 25 30 35

time (me)

0
01

Figure 29 Fluid pressure—punch stroke path: Upper limit

Each point in figure 29 represents the upper bound on pressure at each time stage. The
Forming Limit Diagram (FLD) associated with this pressure profile is shown in figure 30.
The analysis of the supplied finite element model revealed that the contact between the
punch and the blank takes place at 3ms after the start of the forming process. Hence it
was decided to start the Optimization of the pressure profile from time equal to 3ms and

end at the final time of 27ms.

69

 

 

 

 

 

 

 

 

Major Engineering Strain(%‘)

 

 

 

 

 

 

 

 

 

-20 -10 0 10
Minor Engineering StrainCO/o)

Figure 30 FLD Failure and distribution of strains: Upper Limit

The value of pressure at time equal to Oms or at the beginning of the process was fixed at
0.0psi. The value of pressure at the time value of 3ms was fixed at 16.7psi. This is
selected based on the initial pressure profile supplied and the elaborate investigations

made to understand the process.

Out of 10 time stages, the value of pressure at first 2 stages had been fixed. The optimal

pressure values for the remaining 8 stages had to be found. Since the problem is two

dimensional, i.e.. two subsequent pressure values are given as input at once, the

70

optimization must be performed in 4 steps. Hence at each step the optimum values for

two subsequent pressures would be determined.

5.5 Problem Statement

Objective: To minimize the total number of elements in the finite element model of the
sheet metal that fall below wrinkling limit in the FLD at the end of time stage, subjected

to the constraint that strains in all the elements should remain within the tearing limit.

Optimization is performed in four steps and the design variable in each step is the
pressure values for two subsequent time stages. Hence in each time stage the design
variables are
x = {171,192}

where p] and p2 are the pressures.
Reframing the problem in standard form,

min y(x)

x
Subject to

[151:1 SUI
12$XZ $112

where y(x) represents the total number of elements that fall below the wrinkling limit,
x1 and x2 are the design variables and 11, 12 ,u] and u2 are the simple bounds on the

design variables. The constraint that all elements should remain within the tearing limit of

the FLD is taken care by suitable selection of the upper bounds on the design variables.

71

In other words the strain state in all elements, when the design variables are kept at their

upper bounds would stay below the tearing limit in the FLD.

5.6 Optimization

This section describes the process of Optimization preformed in four steps and the results

Obtained in each step.

5.6.] Step 1

Step 1 involves the optimization of pressures for times 6ms and 9ms. The simple bounds

on the design variables for this step are given below:

0.001313 x1 S 0.002

0.00328 3 x2 S 0.004
The lower bounds are selected from the reference pressure path supplied and the upper

bounds are selected such as to satisfy the constraint in the problem. First, Latin hypercube
sampling was employed to sample 8 points in the design domain. The design domain is
two dimensional, each dimension being the pressure for each time stage. The Latin
hypercube sampling and the corresponding values of the Objective function are shown in

table 5.

 

N 11 9‘2 Yr

 

 

1 0.001325 0.003365 490

 

 

 

 

 

 

72

 

 

 

 

 

 

 

 

 

 

 

2 0.001789 0.003982 498
3 0.001516 0.003626 495
4 0.001845 0.003857 498
5 0.001920 0.003798 497
6 0.001591 0.003448 494
7 0.001467 0.003732 494
8 0.001706 0.003489 497

 

Table 5. Initial Designs for step 1

Figure 31 shows the location of the initial sampling graphically. The points are well

distributed in the design domain.

LHS

 

0.00398 ~
0.00388 ~
0.00378 1
0.00368 .
0.00358 «
0.00348 -
0.00338 4

 

0.00328

 

0.00131

Figure 31 Latin hypercube sampling for step 1

0.00151

x1

After sampling the points the metamodel was built.

73

0.00171

0.00191

 

 

5.6.1.1 Discussion of Results

Figure 32 gives a scatter plot of the response values against the data point number for the

initial sampling.

 

 

 

 

499 - ~
2 4
I O
5 a
497 - O Q
(D
U)
S 3
£1495 - O -
E 5 7
o o o
5 493 - .
491 - -
1
b
4891 2 3 4 5 6 7 8

Design point number

Figure 32 Scatter Plot of true responses before Optimization for step 1

The first search point in optimization and the subsequent points are selected based on the
Expected Improvement algorithm. As described earlier the points selected should have a
response value smaller then the current minimum and / or should be in the unexplored
regions of the design domain. The EGO algorithm took 5 additional points to arrive at the

solution.

74

 

 

 

 

500 - -
4
498 ~ 3 o 8 ~
3 o
a) 496 - .
g 3
8 ' 5 7 12
5494 - . . . _
(D
E
1— 492 - -
1
4900 11 -
3 W O l?
488 - ~
2 4 6 8 10 12 14

Design point number

Figure 33 Scatter Plot of true responses after Optimization for step 1

Efficient Global Optimization converged to a solution value of 488. The pressure profile

was then updated with the optimum values from this step.

The FLD of the sheet for termination time up to 9ms with the optimized pressure profile

at the input is shown in figure 34. The strain in many elements is close to and below the

wrinkling limit. This is because of low value of punch force at the start of the process.

75

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

25._, m \_--.-__M____im
é: \ /W i
C
"(:3 2 . -._ ._...,._1L_.-..--._-_...._-_.z_..t_.. If .
35' /
g»
‘C 15~r~*:"‘lr; ‘m
0)
(1)
E
LL] 1 .. ““311: ---------------- -w e --
.6
8
E 5 __._.-...- __- .. - 2 Mg.--

0 tr ,

-1 5 -10 -5 0 5 10 1 5

Minor Engineering Strain(%)

Figure 34 FLD Failure and distribution of strains after step]

5.6.2 Step 2

Step 2 involves the optimization of pressures for times 12ms and lSms. The simple

bounds on the design variables for this step is given below:

0.00496 3 x1 3 0.0055
0.00646 s x2 3 0.0076

Latin hypercube sampling was used to sample 5 points in the design domain. . The Latin
hypercube sampling and the corresponding values of the objective function are shown in

table 6.

76

 

 

 

 

 

 

 

 

 

 

n 1‘1 x2 yi
1 0.005446 0.006624 264
2 0.005145 0.007232 265
3 0.005006 0.007114 264
4 0.005317 0.007567 264
5 0.005196 0.006827 265

 

 

 

Table 6. Initial Designs for step 2
LHS

 

0.00746 ~
0.00726 .
Q: 0.00706 1
0.00686 1

0.00666 1

 

0.00646

 

 

0.00496

0.00516 0.00536

x1

0.00556

Figure 35 Latin hypercube sampling for step 2

Figure 35 shows the location of the initial sampling graphically.

77

5.6.2.1 Discussion of Results

Figure 36 gives a scatter plot of the response values against the data point number for

the initial sampling.

 

 

 

 

266 . . .
2
265 - o i 50
9’6
8
g 2640 . -
3
F.
263 - .
2621 2 3 4 5

Design point number
Figure 36 Scatter Plot of true responses before Optimization for step 2
The variation in response in the present step being very minimal the metamodel tuned out
very accurate. The lack of large variations in the fimctional response is due to two
reasons. First, as the forming process continued from step one, the contact between the
punch and the sheet increased and this caused the strains to move towards the tearing
limit in the FLD. There was no addition to the number of elements below the wrinkling

limit but the number started decreasing due to the increased contact.

78

Efficient Global Optimization algorithm took just two additional points to converge to the

solution. Table 7 shows the two points generated by EGO and their corresponding

response values.

 

N 3‘1 X 2 .Vi

 

 

6 0.005439 0.0076 263
7 0.005439 0.0076 263

 

 

 

 

 

 

 

Table 7. Search points generated by EGO for step 2

U)
C)

 

 

L...~——.~. w -— .. -. “~—

 

 

/
/

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

32‘
E
g 20g..- . -
(D
U)
E \.__
615 - --—
e)
.E
8‘
UJ 1 ..‘QK.
.5
5 m- At...
0 ' T _
~15 -10 -5 0 5 10 15

Minor Engineering Strain(%)

Figure 37 FLD Failure and distribution of strains after step 2

79

After the completion of this step, the optimal pressure profile was updated for termination
time up to 15ms. The FLD of the sheet for termination time up to lSms with the

optimized pressure profile at the input is shown in figure 37.

5.6.3 Step 3

Step 3 involves the optimization of pressures for times 18ms and 21ms. The simple
bounds on the design variables for this step is given below:

0.0088 5 x1 5 0.0105
0.0103 5 x2 3 0.0110

The Latin hypercube sampling and the corresponding values of the objective function are

shown in table 8.

 

 

 

 

 

 

N 3‘1 3‘2 Yr

1 0.009480 0.01000 189
2 0.009108 0.01009 189
3 0.008879 0.010286 189
4 0.009153 0.010194 189

 

 

 

 

 

 

Table 8. Initial designs for step 3

The response of the objective function remained the same with different pressure inputs.
As discussed previously the increased punch force caused the strains in all the elements

to be pushed towards the tearing limit. Also it can be inferred that the pressure values at

80

this step have no control over the number of elements falling below wrinkling limit. This
negates the possibility of building a metamodel for this step and hence obtaining
optimum pressure values for times 18ms and 21ms. As a result, one of the sampled points '

(sample point four) was selected randomly and the optimum pressure curve was updated.

5.6.4 Step 4

This step involves the optimization of the pressure profile for times 24ms and 27ms. The
simple bounds are shown below:

0.01040 S x1 5 0.01070
0.01090 S x2 S 0.01110

The Latin hypercube sampling and the corresponding values of the objective fimction are

shown in table 9.

 

 

 

 

 

 

 

1‘ x1 3‘2 Yr

1 0.01096 0.01109 148
2 0.01053 0.01103 145
3 0.01041 0.01095 142
4 0.01055 0.01097 146

 

 

 

 

 

Table 9. Initial Designs for step 4

81

5.6.4.1 Discussion of Results

Efficient Global Optimization algorithm took 2 additional points to converge to the

solution. The two additional points generated by EGO are shown in table 10.

 

 

 

 

N 11 3‘2 Yr
5 0.01044 0.01098 140
6 0.01043 0.01098 139

 

 

 

 

 

 

Table 10. Search points generated by EGO for step 4

 

 

 

 

  

 

 

 

 

Major Engineering Strain(%)
VI

 

 

 

 

 

 

 

10.2-..”
5... .._..
1
~20 ~10 0 10

Minor Engineering Strain

Figure 38 FLD Failure and distribution of strains for the optimized pressure profile

82

The algorithm converged to the solution with just two additional points. This could be
due to smaller variations in the response among the sampled points. With the completion
of this step, the optimal fluid pressure—punch stroke path was obtained for the total
process termination time of 27ms. The FLD showing the distribution of strains in all the
elements at the completion of the forming process (time: 27ms) is shown in figure 38.

The optimal solution obtained after four steps is equal to 139 i.e., at the completion of the
forming process the total number of elements that stayed below wrinkling limit is 139.
When compared to the result from the work of Abedrabbo [1] wherein this number was
equal to 163. This shows that the optimization approach suggested could be used to find

an optimal fluid pressure—punch stroke path in sheet metal hydroforming.

83

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

BIBLIOGRAPHY

N.E. Abedrabbo, “Experimental and Numerical Investigations of Stamp
Hydroforming and Ironing of Wrinkling in Sheet Metal Forming”, A Master’s
Thesis, Michigan State University, 2002.

Z.A. Bedewi, N.E., Kan, CD, and Marzougui, D., “Development and Evaluation of
a C-1 500 Pick-Up Truck Model for Roadside Hardware Impact Simulation,”
Presented in the FHWA Simulation Conference, Langley, VA, (0 F HWA and GWU
July 1996.

Z.A. Bedewi, A., Kan, A., and Marzougui, D., “Validation of a Non-linear Finite
Element Vehicle Model Using Multiple Impact Data,” Presented and Published in
the International ASME Winter Annual Congress and Exposition, Atlanta GA,
November 17-22, 1996.

D.G. Brooks and William A. Verdini, “ Computational Experience with
Generalized Simulated Annealing Over Continuous Variables,” American Journal
of Mathematical and Management Sciences, Vol.8, pp.425-449, 1988.

DD. Cox and Susan John, “ SDO: A Statistical Method for Global Optimization,”
In M.N. Alexandrov and M.Y. Hussaini (Eds), Multidisciplinary Design
Optimization: State of The Art, pp.315-29. SIAM, Philadelphia, 1997.

N. Cressie, “ Spatial Prediction and Ordinary Kriging,” Mathematical Geology, 20
(4), pp.405-421, 199

A. Farhang-Mehr, S.Azarm, “A Sequential Information-Theoretic Approach to
Design of Computer Experiments,” 9th AIAA/ISSMO Symposium on
Multidisciplinary Analysis and Optimization, Atlanta, Georgia, 2002

CL. Ignatovich, A.R. Diaz, C.A. Soto. Strategies In Design For Enhanced
Crashworthiness. Proceedings of DETCOO, 2000 ASME Design Engineering
Technical Conference, September 10—13, 2000 — Baltimore, Maryland.

T. Jansson, L. Nilsson and M. Redhe, “Using Surrogate Models and Response
Surfaces in Structural Optimization,” Nordic Seminar of Computation Mechanics,
Lund, October 18-20, 2001.

R. Jin, Wei, Chen and Timothy W. Simpson, “Comparative Studies of
Metamodeling Techniques Under Multiple Modeling Criteria,” Journal of Structural
Optimization, Vol. 23, No.1, pp.1-l3, 2001.

DR. Jones, Schonlau, M., Welch, W., “Efficient Global Optimization of Expensive
Black-Box Functions”, Journal of Global Optimization, Vol.13, pp.455-492, 1998.

84

[12] DR. Jones, “A Taxonomy of Global Optimization Methods Based on Response
Surfaces,” Journal of Global Optimization, Vol.21, No.4, pp.345-383, 2001.

[13] M.M.Kamal and J.A.Wolf , “Modern Automotive Structural Analysis,” Van
nostrand Reinhold company, 1982.

[14] SM. Kang, Lebbelici, and Yusuf, “CMOS Digital Integrated Circuits: Analysis
and Design,” second edition, McGraw-Hill Higher Education, 1998.

[15] CH. Kim, A.R. Mijar and J .S. Arora, “Development of Simplified Models for
Design and Optimization of Automotive Structures for Crashworthiness,” Structural
and Multidisciplinary Optimization, Volume 22, Issue 4, pp.307-321, 2001.

[16] C.H.Kim , “ Development of Simplified Models for Automotive Crashworthiness
Simulation and Design Using Optimization,” A Ph.D. dissertation, university of
Iowa, 2001.

[17] SW. Lo, T.C.Hsu and W.R.D. Wilson, “An Analysis of the Hemispherical Punch
Hydroforming Process,” Journal of Materials Processing Technology, Vol.37,
pp.225-239,1993.

[18] EA. McClintock, “A Criterion for Ductile Fracture by the Growth of Holes,”
Journal of Applied Mechanics, Vol. 35, pp.363-371,1968 -

[19] T.J. Mitchell and Max D. Morris, “The Spatial Correlation Function Approach to
Response Surface Estimation,” Proceedings of the 24th conference on Winter
Simulation, p.565-571, December 13-16, 1992, Arlington, Virginia, United States.

[20] T.J. Mitchell and Max D. Morris, “Bayesian Design and Analysis of Computer
Experiments,” Statistica Sinica , Vol.2, pp.359-379.

[21] DC. Montgomery and EA. Peck, “Introduction to Linear Regression Analysis,”
second edition, Willey series ,1992

[22] RH. Myers, Andre I. Khuri and Walter H. Carter, Jr., “Response Surface
Methodology: 1966-1988,” Technometrics, Vol. 31 (2), pp. 137-157.

[23] W. J. Roux, Nielen Stander and RT. Hafika, “Response Surface Approximations
for Structural Optimization,” International Journal for Numerical Methods in
Engineering, Vol.42, pp.517-534, 1998.

[24] J. Sacks, William J. Welch, Toby J. Mitchell and Henry P. Wynn, “Design and

Analysis of Computer Experiments,” Statistical Science, Vol.4, Issue.4, pp.409-
423, 1989.

85

[25] RR. Salagame and Russel R. Barton, “Factorial Hypercube Designs for Spatial
Correlation Regression,” Journal of Applied Statistics, Vol. 24, No.4, 1997, 453-
473.

[26] M.J. Sasena, “Optimization of Computer Simulations Via Smoothing Splines and
Kriging Metamodels” A Master’s Thesis, University of Michigan, 1998.

[27] M.J. Sasena, Panos Papalambros and Pierre Goovaerts, “Exploration of
Metamodeling Sampling Criteria for Constrained Global Optimization,” Journal of
Engineering Optimization, Vol.34, No.2, pp. 263-278, 2002.

[28] M.J. Sasena, “Flexibility and Efficiency Enhancements for Constrained Global
Design Optimization with Kriging Approximations,” A Ph.D. Dissertation,
University of Michigan, 2002.

[29] K. Scheinberg, “Derivative Free Optimization Method,” IBM Watson Research
Center.

[30] M. Schonlau, "Computer Experiments and Global Optimization" A Ph.D.
Dissertation, University of Waterloo, 1997.

[31] M. Schonlau, Welch, W., Jones, D., “A Data-Analytic Approach to Bayesian
Global Optimization”, American Statistical Association Proceedings, Section of
Physical Engineering Sciences, 186-191, 1997.

[32] M. Schonlau, Welch, W., Jones, D., "Global Versus Local Search in Constrained
Optimization of Computer Models," in New Developments and Applications in
Experimental Design, N. F loumoy, W.F. Rosenberger, and W.K. Wong (editors),
Institute of Mathematical Statistics, Hayward, California, Vol. 34, 11-25, 1998.

[33] T.W. Simpson, Jesse Peplinski, Patrick N. Koch and Janet k. Allen, “ On the Use of
Statistics in Design and Implications for Deterministic Computer Experiments,”
Proceedings of DETC’97, 1997 ASME Design Engineering Technical Conference.

[34] T.W. Simpson, Timothy M. Maurey, John J. Korte and Farrokh Mistree,
“Comparison of Response Surface and Kriging Models for Multidisciplinary
Optimization,” 7th Symposium on Multidisciplinary Analysis and Optimization, St.
Louis, MO, AMA, Vol.1, pp.381-391, 1998.

[35] T.W. Simpson, Janet k. Allen and Farrokh Mistree, “Spatial Correlation
Metamodels for Global Approximation in Structural Design Optimization,” 1998
ASME Design Engineering Technical Conferences - Design Automation
Conference, Atlanta, GA, ASME, Paper No. 98DETC/DAC-5615, 1998.

[36] T.W. Simpson, “A Concept Exploration Method for Product Family Design” A
Ph.D. Dissertation, Georgia Institute of Technology, 1998.

86

[37] T.W. Simpson, Dennis K.J. Lin and Wei Chen, “Sampling Strategies for Computer
Experiments: Design and Analysis,” International Journal of Reliability and
Applications, Vol.2, No.3, pp.209-240, 2001.

[38] T.W. Simpson, Timothy M. Mauery, John J. Korte and Farrokh Mistree, “Kriging
Models for Global Approximation in Simulation-Based Multidisciplinary Design
Optimization” AIAA Journal, Vol. 39, No. 12, December 2001.

[39] C.A.Soto and A.R.Diaz , “ Basic Models for Topology Design Optimization in
Crashworthiness Problems,” Proceedings of DETC 99, Las Vegas.

[40] A. Srivastava, Hacker, K., Lewis, K., and Simpson, T., “Development of a Kriging
Based Surrogate Approximation Method for Large Scale Systems,” International
Journal for Product and Process Improvement, 1999.

[41] J. Tirosh, S.Yossifon, R. Eshel and A.Betzer, “Hydroforming Process of Uniform
Wall Thickness Products,” ASME Journal of Engineering for Industry,
Vol.99,pp.685-691,1997.

[42] G. Venter, Raphael T. Haftka and James H. Starnes, Jr., “Construction of Response
Surfaces for Design Optimization Applications,” Proceedings of the 6th
AIAA/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization,
Bellevue, Washington, September 4 - 6, 1996, Paper No. 96-4040, Part 1, pp. 548-
564.

[43] W.J. Welch, R.J. Buck, J. Sacks, H.P. Wynn, W.J. Mitchell, and MD. Morris,
“Screening, Predicting, and Computer Experiments,” Technometrics, 34:15—25,
1992.

[44] R.J.Yang, L.Tseng, L.Nagy, and J.Cheng , “Feasibility Study of Crash
Optimization,” Advances in Design Optimization, DE-Vol.69-2, 1994.

[45] Y. Youssef, and J. Denault, “Thermoformed Glass Fiber Reinforced Polypropane:
Microstructure, Mechanical Properties and Residual Stresses,” Polymer
Composites, Vol. 19, No.3, pp.301-309,1998.

[46] AK. Zaouk, Nabih E. Bedewi, Cing-Dao Kan, and Hansjoerg Schinke, “Evaluation
of a Multi-purpose Pick-up Truck Model Using Full Scale Crash Data with
Application to Highway Barrier Impacts,” 29th International Symposium on
Automotive Technology & Automation, Florence, Italy, 1996.

87