Open access

Concurrent Process-Product Design Optimization Using Coupled Nonlinear Finite-Element Simulations

Written By

Ali Najafi and Masoud Rais-Rohani

Submitted: February 14th, 2012 Published: October 10th, 2012

DOI: 10.5772/47852

Chapter metrics overview

2,245 Chapter Downloads

View Full Metrics

1. Introduction

Since such FE simulations and the accompanying structural design optimization studies rely on computer-based geometric model of the structural component and the tabulated material properties, the initial state (e.g., internal stresses and strains) of the component/material is most often ignored, which results in exclusion of the manufacturing process effects on the product (e.g., plastic strain, residual stress, thinning, springback, etc.). The selected manufacturing process and the choice of process parameters can also change the material microstructure (e.g., dislocation density, texture), thereby affecting the macroscale (e.g., stress-strain) behavior of the material and the structural component (Najafi et al., 2012; Oliveira et al., 2006).

A practical way to alleviate this shortcoming is to perform coupled process-performance simulations in a sequential manner whereby both changes in the material and component can be properly modeled and tracked from one stage to the next for a more accurate prediction of the structural performance measures (Noels et al., 2004). Consequently, process parameters can be evaluated based on both the process objectives and performance criteria. Coupling of the material, process, and performance models is an important step in modeling the actual physical behavior of the material and structure while facilitating the integrated material-process-product design (Olson, 1997; McDowell et al., 2007; Acar et al., 2009).

Coupled quasi-static analyses can be performed in some implicit finite element analysis (FEA) codes such as Abaqus/Standard where the boundary conditions can be changed and the material state can be tracked and passed from one solution step to the next by using the restart option, for example. The same capability also exists in some explicit FEA codes such as LS-DYNA and Abaqus/Explicit for coupled transient dynamic simulations. The increasing demand for coupled simulations has resulted in many commercial software codes to have an option to map some solutions as the initial state. However, when some combinations of dynamic and quasi-static analyses (i.e., explicit and implicit solvers) are required, a separate data flow management (DFM) program and strategy are required. The DFM procedure becomes more complex in the context of design optimization when the coupled analyses need to be performed numerous times for design space exploration in search of the optimum design. A particular case being considered in this chapter is the concurrent design of a coupled process-product system where both manufacturing process and product performance attributes are to be optimized by finding the optimum values of the manufacturing process on product design variables.

Coupled process-product (performance) simulation studies include those that have considered the manufacturing effects associated with forming and springback on crush/crash performance of tubular components (Oliveira et al., 2006; Kellicut et al., 1999; Kaufman et al., 1998; Grantab, 2006; Krusper, 2003). For example, Kellicut et al. (1999) considered springback, thinning, and other parameters such as plastic strain and residual stresses in bending-crush simulations of hydroformed tubes and showed that the plastic strain has the greatest effect on the crush behavior. Mayer (2004) and Williams et al. (2005) also performed integrated hydroforming-crush simulations, whereas Ryou et al. (2005) used ideal forming solution to extract the stress and strain responses from forming process and a hybrid membrane/shell method to pass the information to impact simulation. They improved the computation time while preserving the accuracy of the FE simulations. Simunovic and Aramayo (2002) studied the crash response of energy absorbing components of the ultra light steel auto body vehicle models and showed that by including the history effects the energy absorption properties can change even though the difference in the overall response was relatively modest.

Oliveira et al. (2006) and Williams et al. (2005) performed experimental and computational study of s-rail tubes and discovered that both the maximum and mean crush force values will change as a result of the manufacturing process effects. Bottcher and Frik (2003) did a similar study and showed that metal forming data is required in crash simulation of front rail panel of a vehicle model, especially in high strength dual phase steel due to its rapidly hardening characteristic. Krusper (2003) and Dagson (2001) performed analysis on a simple bar while considering the springback response of the material. Most of the studies cited only considered a material model with isotropic hardening while a few included the effect of kinematic hardening on the crush response. Recently, Williams et al. (2010) studied the effect of combined isotropic/kinematic hardening and strain rate sensitivity along with an anisotropic yield surface to study the crush behavior of hydroformed tubes. They showed that the combination of these parameters increases the capability of models to correctly predict the energy absorption performance of the crush tube. However, they did not investigate the capabilities of the model in coupled process-performance simulations.

This chapter presents the different steps in performing sequential coupled nonlinear FE simulations and their application in multi-objective process-product design optimization of thin-walled structures. The sheet stamping simulation includes both the deep drawing and the springback stages of the manufacturing process. Performance simulation considers the energy absorption response characteristics of the component under an impact (crash) loading condition. The coupled simulations involve both nonlinear explicit and implicit FEA. The developed computational framework is used for the analysis and design of a double-hat tube. A design sensitivity analysis is performed to investigate the effect of manufacturing process parameters and geometric attributes on the process and performance responses that are affected by the manufacturing process and geometric design parameters. The ensuing nonlinear design optimization problem is cast in a multi-objective formulation and solved for Pareto optimum design points using a multi-objective genetic algorithm.


2. Modeling of manufacturing effects

Constitutive models describe the stress-strain relationship for a given material and the influence of various factors such as temperature and strain rate. Plastic deformation requires variables that define the history of stress and temperature in the material. The history can be defined through functional analysis and mathematical theories known as the theory of material with memory (Lubliner, 2006). Manufacturing effects in most continuum-level material models are considered implicitly through the state variables defined in the model. For example, in the classical plasticity model, the total strain is written as an additive decomposition of elastic strain and plastic strain. Considering a piecewise linear isotropic hardening law derived from the stress-strain data, the equivalent plastic strain would represent, in a limited sense, the history or the manufacturing effect. Inclusion of manufacturing effects in the continuum plasticity models emerges by specifying the initial state for the numerical integration of the evolution equation of the state variable.

The constitutive relation of the rate dependent classical plasticity model with isotropic hardening can be represented through the following equations:

  • Elastic stress-strain relationships in three-dimensional space

  • Yield surface (closure of elastic domain in the stress space)

  • Flow rule and hardening law

ϓ̥vp=ϙ̥ fϡ,Ϙϡ, Ϙ̥=ϙ̥E3
ϙ̥ϡ,Ϙ=1Cϡϡy1P-1 if fϡ,Ϙ0 0 if fϡ,Ϙ<0 E4

where the bold symbols ϡand ϓrepresent the second-rank stress and stain tensors, respectively, and Cis the forth-rank tensor of isotropic elasticity. The superscript vp designates the viscoplastic component. The ̥ terms indicate time derivatives. The yield stress is denoted by ϡy with parameters H and ҝ representing isotropic hardening modulus and plastic multiplier, respectively. Moreover, Cand Pin the hardening law equation are known as material constants associated with rate sensitivity, which are found based on experimental data on a specific material. As mentioned previously, the manufacturing effects emerge as the consequence of numerical integration of the model described by Eqs. (1) through (4).

The numerical integration scheme in the deviatoric space is defined based on updating the current state of stresses and state variables (with subscript n) by calculating their increment. The deviatoric part of total strain tensor eis updated as


where Iis the identity tensor and tris the trace of the second rank tensor. The trial stress tensor sn+1t is calculated from the Hook’s law given by


where the enp is the plastic strain at the nth increment and Ϛis the shear modulus. Then, the calculated stress magnitude value at n+1 increment sn+1t is evaluated with the yield surface fn+1t at the same increment. If fn+1t=sn+1t-23ϡy+0 is met, the material is in the elastic region and the trial values are the admissible stress values and, therefore, the plastic multiplier ϙ̥ or βϙn+1 is zero with


And if fn+1t>0, the trial values should be corrected due to an additional plastic strain and the yield surface is approximated through the first-order Taylor series expansion as


where βϙn+1 is computed based on Newton-Raphson scheme (Wang and Budiansky, 1978; Simo and Hughes, 2000; Souza et al., 2008). Therefore, the next state of material (with subscript n+1) is derived by adding the increments βto the previous state of stress as

Ϙn+1=Ϙn+23 βϙn+1E9
ϓn+1vp=ϓnvp+βϙn+1 nn+1E10
βϡn+1=ϙ¯ trβϓn+1I+βsn+1t-2Ϛ βϙn+1 nn+1E11

where trrepresents the trace of the tensor, trial values are denoted with t superscript, and the normal to the yield surface is found asnn+1= sn+1tsn+1t. By havingϘn=ϓnvp, the manufacturing effect can be considered through a non-zero state of ϡn and ϓnvp at the beginning of the analysis as ϡ1 andϓ1vp.

Traditionally, the initial values for the first increment are considered to be zero assuming that the material deformation starts from a conceptual pristine state. However, incorporating proper initial values emanating from forming and springback simulations in the crush analysis couples the process and performance from the material standpoint.


3. Manufacturing process and product performance simulations

Coupled process and performance simulations are conducted sequentially by using Abaqus/Explicit for the deep drawing (loading) simulation, Abaqus/Standard for the springback (unloading) simulation under isothermal condition, followed by Abaqus/Explicit for the crush simulation. The illustrative example considered in this chapter is a double-hat, thin-walled tube that is modeled by joining two identical single hat sections.

3.1. Deep-drawing simulation

For forming or deep-drawing simulation, two sets of blank/holder/die/punch geometries are defined in the FE model with the model for one set shown in Fig. 1. A single die set is mirrored with respect to blank plane to simulate the forming of both pieces simultaneously. There are two basic steps in using the explicit FEA for the forming simulation of each single hat section depending on the specified boundary conditions.

The first step defines the gripping of the blank between the die set and the holders. While keeping the punch and dies fixed in their respective positions, the holding forces are increased linearly from zero to the specified value matching the corresponding manufacturing process control parameter. In this stage, the kinematic contact formulation is used because of the computational efficiency of the formulation (Rebel et al., 2002; Oden and Kikuchi, 1982; Oden and Pires, 1983; Bayram and Nied, 2000; Wriggers, 2006). Due to the simplicity of the geometry, die and holder surfaces are defined through standard analytical rigid surfaces. Contact surfaces are defined on both surfaces of the blank by considering the surface offset due to the blank thickness. Penalty formulation is used in tangential contact and a friction coefficient is defined as a manufacturing process parameter for an equivalent representation of both surface roughness and draw beads.

Figure 1.

Sheet-forming FE model of a single hat section

Figure 2.

Formed hat section with the boundary conditions used for the subsequent springback simulation

The second step performs the deep drawing simulation by applying a constant velocity to the punch that is configured in the model to match the final geometry of the product (excluding the springback effect). Hence, the results of the previous step are directly transferred to this step where the boundary conditions on the fixed punch in the direction normal to the blank surface are removed and a constant velocity is applied to the punch to form the single hat section as shown in Fig. 2. The punch velocity is assumed to be constant for a linear displacement. This setup does not impose a uniform strain rate in all the elements. Thus, rate sensitivity of the material will not have a uniform effect on the structure. Dies remain clamped and the holders are fixed in all degrees of freedom except the direction parallel to the punch movement. In this direction, the constant holding force is applied to preserve the constant gripping force throughout the drawing process.

Rupture and thinning are two responses that can be calculated from the results of the deep-drawing simulation. Rupture is calculated by extracting the principal major and minor plastic strains in each element and taking the difference between the major strain and that extracted from the forming limit diagram (FLD) using the following equation

R=i=1IRi2=i=1I(ϓ1i-ϳ(ϓ2i))2 ϓ1i>ϳ(ϓ2i) 0 ϓ1iϳ(ϓ2i) E13

where ϳϓ2is the equation representing FLD curve and ϓ1i and ϓ2i are the principal major and minor strains at each integration point through the thickness, respectively, which are calculated in the deep-drawing simulation when the termination time reaches the limit. The parameter I represents the total number of integration points in the mesh. The FLD (Lee et al., 2008) that is used in this study is that for AZ31 magnesium alloy sheet and it is assumed to behave linearly in both compressive and tensile plastic strains as shown in Fig. 3.

Figure 3.

FLD for AZ31 at two different strain rates (Lee et al., 2008)

Thinning is measured by taking the difference between the final thickness of each element and its corresponding initial value and is calculated using a single metric, T defined as


where to and ti are the initial and final shell thicknesses, respectively, and N is the total number of elements in the mesh. Since the shell thickness in the blank is assumed to be constant for all the elements, tois the shell thickness assigned to the elements.

3.2. Springback simulation

The state variables and geometric information from the deep drawing simulation are transferred and treated as the initial state in the unloading stage (i.e., removal of all the tooling parts from the workpiece) for the springback analysis. Springback process is modeled as a quasi-static problem considering the stress distribution captured from deep drawing, dynamic effects, and the contact conditions. Additionally, all the rigid surfaces including punch, dies, and holders are removed from the FE model, which makes the model more suitable for implicit FEA considering the quasi-static nature of the springback phenomenon and the absence of highly nonlinear factors in the model. Convergence of the nonlinear implicit FEA are guaranteed by defining the boundary conditions in such a way that the two edges of the hat section are held fixed perpendicular to the actual normal surface. As shown in Fig. 2, the equilibrium condition is achieved by constraining the model in all the transverse directions. The boundary conditions defined in this stage are designed such that the effect of the force required to assemble a non-fitted double-hat section is already considered. The residual stresses and geometric attributes are updated during quasi-static analysis while the other computational state variables such as plastic strain remain unchanged. Similar to the deep-drawing simulation, the springback analysis is performed separately and simultaneously on the two identical single hat sections. There is no interaction between the two hat sections, however, in both the deep drawing and springback simulations. The two hat sections are then assembled in the next stage to produce a double-hat crush tube.

Springback is calculated by comparing the nodal coordinates obtained in the last step of the deep-drawing simulation with those in the last step of the springback simulation. A single springback metric, S representing the deviation of the nodal coordinates is calculated as

S=MaxXi-Xo2+Yi-Yo2+Zi-Zo2 i=1,NNE15

where X,Y,Zare the Cartesian coordinates with subscripts i and o representing the end of springback and deep-drawing stages, respectively with NN as the total number of nodes in the mesh.

An in-house FORTRAN code is used to automate the procedure to extract the rupture and thinning results from the Abaqus binary file, calculate the principal strains, and incorporate the equations mentioned above (without using Abaqus CAE) in the deep-drawing and springback simulations.

3.3. Joining and trimming

The two hat sections are joined longitudinally using tie contact formulation and trimmed by removing the outer flange elements as shown in Fig. 4. The tie contact constrains the master and slave surfaces similar to the multiple constraint points when the clearance between two surfaces is below the tolerance specified as input. If the surfaces are out of the prescribed tolerance, the interaction becomes a contact formulation. A preliminary study showed that switching the master and slave surfaces would not affect the crushing behavior of the tube. Once the distance between the two surfaces becomes more than the clearance tolerance, constraints are removed and contact formulation is activated similar to the conventional contact definition. It is worth noting that this kind of joining represents a fictitious weld seam along each joint line since the thermo-mechanical process involved in an actual welding process is not modeled here.

Figure 4.

Joining and trimming of the two formed sections to generate the final double-hat tube geometry

3.4. Crush simulation

The explicit FEA for crush simulation is performed by holding one end of the tube fixed while applying an axial load through a moving rigid wall at the other end. The rigid wall moves with a constant speed to simulate constant loading rate defined with a prescribed displacement as shown in Fig. 5. The component geometry, residual stresses, and state variables at the start of crush stage are those accumulated through all the manufacturing stages discussed previously.

Figure 5.

General setup for crush simulation

Six contact interaction sets among the elements are defined in the crush simulation, including interactions between the lower hat section and the rigid wall, the upper hat section and the rigid wall, interaction between the upper and lower hat sections, tie contact between the assembly edge of the upper and lower surfaces, and separate self contact interaction for the upper and lower hat sections. For all of the aforementioned contact interaction sets, penalty function formulation in both normal and tangential directions is used. Despite the computational cost, penalty function provides a proper flexibility for the explicit FEA to find a stable time step that is affected by severity of the contacts. Moreover, the maximum ratio of thickness-to-element length is used to overcome the difficulty of the fine mesh density that results in relatively thick shell elements.

The contact force history of the rigid wall during the crush simulation is used to calculate the maximum crush force, Pmax while the mean crush force, Pm is found by dividing the area under the crush force-displacement curve by the effective crush distance. The mathematical equation for finding the mean crush force is expressed as

Pm=1ϒeff0tFtDtdt E16

where Ftis the instantaneous contact force normal to the rigid wall surface and D(t)is the instantaneous cross-head axial displacement of the rigid wall. The final cross-head displacement of the rigid wall represents ϒeff and it is assumed to be 125 mm, or 50% of the tube length. The maximum crush force is the largest value of Ftfound after applying the SAE filtering of 60 Hz to filter the noise in the force data. These results are extracted and filtered using a Python scripting application available in Abaqus and used as input to an in-house FORTRAN code to calculate the mean crush force values.

Figure 6.

Sequence of coupled nonlinear FE simulations

Figure 6 shows the general four-step sequence of coupled simulations from the initial sheet forming to crush. The rupture and thinning responses are the outputs of step 1, springback is the result of step 2 while the mean and max crush force values are the outputs of step 4. Other than the joining of the two hat sections and removal of the excess tabs, no changes are induced in the tube model in step 3.

3.5. Material model

The material model uses piecewise linear isotropic hardening. The constant for the linear kinematic hardening is calculated based on the slope of a line connecting two adjacent points on the stress-strain curve. The material model uses von Misses yield surface and a one-dimensional stress-strain input is considered as equivalent von Misses stress versus effective plastic strain. Coupling scheme is utilized by transferring residual stresses and the equivalent plastic strains as the material state variables. The yield surface expands due to the isotropic hardening assumption in the model; therefore, the instantaneous yield point varies during the loading process. The yield point at the end of the forming simulation is captured by finding the plastic strain.

The AZ31 magnesium alloy sheet data is used for all the simulations. For modeling the stress-strain response at various strain rates, the stress-strain curves for two extreme rates (i.e., 1 s-1 and 4300 s-1) are considered with those for the other rates found through interpolation. The elastic modulus, Poisson’s ratio, and density are 45 GPa, 0.33, and 1.738 kg/m3, respectively. The true stress-true strain curves for the two extreme strain rates are shown in Fig. 7. Adiabatic heating is not considered in any of the simulations.

Figure 7.

AZ31 magnesium alloy sheet stress-strain curves for two different strain rates

3.6. Effect of manufacturing process on product performance

To examine the role of the history effects on the axial crush response, two simulation cases are compared. In the first case, separate stand-alone performance simulation that does not include any history effects is performed, whereas in the second case, a sequential coupled process-performance simulation is performed that includes residual stresses, plastic strains, thinning, and springback information from process simulation (as manufacturing effects) together with a piecewise linear isotropic hardening material model in performance simulation. The 250-mm long tubes are modeled using the plane-stress shell element formulation. They are held fixed at one end and axially loaded with a flat rigid wall at the other end that is moving with a constant speed of 5 m/s. Both self-contact and surface-to-surface contact between rigid wall and tube are specified. A classical multi-linear kinematic hardening material model is used for this comparison.

Figure 8.

Crush behavior with and without consideration of manufacturing process effects

Figure 8 shows the two crush force-distance curves as well as the corresponding crush modes. The results clearly show that both the crush response and the collapse mode change due to the inclusion of the history effects. The maximum crush force increases by approximately 10% from 77 kN to 85 kN, whereas the mean crush force decreases from 20 kN to 18 kN when the history effects are considered.


4. Sensitivity analysis

A design sensitivity analysis is helpful in capturing the main effects of the individual process and product design variables on both the manufacturing and performance attributes (Najafi 2011). The product design variables are the tube cross-sectional dimensions (i.e., width, height of a single hat section, corner radius, and initial blank thickness) shown in Fig. 9, whereas the process design variables are the holding force, punch velocity and workpiece / die set friction coefficients.

Figure 9.

Description of geometric design variables for a 250-mm long tube

The friction coefficients for the holders, dies, and punch are assumed to be equal but can be treated as different design variables. The width design variable defines the punch width, the corner radius defines the die and holder corner radii, the thickness design variable is assigned directly to the shell elements that define the blank and the height design variable is controlled by the punch travel distance in the direction normal to the blank surface (this parameter also determines the simulation termination time as well as the prescribed punch velocity). Holding force defined as a manufacturing process parameter is the amount of maximum incremental force in the first step of deep drawing. The rate of holding force application is kept constant in all the simulations. Punch velocity is assumed to be constant in the direction perpendicular to the sheet metal; this parameter along with the height determines the deep-drawing simulation termination time. Friction coefficients are assigned to the contact tangential definition. Both kinematic and penalty tangential contact formulations produced the same response in the deep drawing simulation.

Width (mm)Height (mm)Corner Radius (mm)Thickness (mm)Holding Force (kN)Punch Velocity (m/s)Friction Coefficient
Radius ±15%5527.55.751.75306.00.225
Force ±15%5527.55.01.7534.56.00.225
Velocity ±15%5527.55.01.75306.90.225
Coefficient ±15%5527.55.01.75306.00.25875
Upper Bound70357.52.55010.00.35
Lower Bound40202.51.01002.00.1

Table 1.

The values assigned to design variables for sensitivity analysis

The sensitivity results are shown in Fig. 10. In each case, the sensitivity values are found by perturbing one design variable by +/-15% from its corresponding average value while holding the remaining design variables fixed at their respective average values shown in bold numbers in Table 1. Rupture is found to be most sensitive to the sheet thickness followed by the corner radius. In contrast, the friction coefficient, punch velocity, and holding force appear to have minimal effect. The rupture response was found to have a direct relationship with some parameters such as thickness and punch velocity and inverse relationship with others, corner radius being the most notable. The global measure of thinning is affected the most by changes in the corner radius, followed by blank thickness and height. In comparison, the manufacturing process parameters appeared to be less influential. Springback response is most sensitive to changes in blank thickness followed by the corner radius and friction coefficient. Both the maximum and mean values of the crush force increase as a result of increasing the blank or tube thickness. Generally, sensitivities to the geometric parameters seem to be greater than those of the manufacturing process parameters.

Figure 10.

General sensitivity of rupture (a), thinning (b), springback (c), maximum crush force (d), and mean crush force (e) to the process and product design variables


5. Multi-objective design optimization

In a traditional process optimization problem, the manufacturing process objectives are optimized by varying the process control parameters (Sun et al., 2010; Wei and Yuying, 2008; Hu et al., 2008; Konak et al., 2006) whereas in product optimization, the geometry (e.g., shape, size) is altered to enhance the product performance. However, in a coupled process-performance optimization problem as considered here, both manufacturing- and performance-level attributes are optimized. When faced with competing objectives, the resulting multi-objective optimization problem becomes one of finding not just one but a collection of non-dominated design points that form the Pareto frontier. By specifying a particular target value for each objective, the multi-objective optimization problem is expressed as

min x=x1,,x7min x=x1,,x7Rx-TR2, Tx-TT2, Sx-TS2, (Pmx-TPM)2,Pmaxx-TPX2, (Mx1,,x4-TM)2E17
s.t.40 x1 70 mm20 x2 35 mm2.x3 7.5 mm1.0 x4 2.5 mm10 x5 50 kNx6 10 m/s0.x7 0.35E18

where design variables x1 to x7 represent width, height, corner radius, thickness, holding force, punch velocity, and friction coefficient, respectively with each having both lower- and upper-bound side constraints. Also, R(x) is rupture, T(x) is thinning, S(x) is springback, Pm(x) is the mean crush force, Pmax(x) is the maximum crush force, and M(x) is mass with TR, TT, TS, TPM, TPX, and TM as the corresponding target values, defined later in this section. The blank length is always equal to the tube length of 250 mm, whereas the blank width is selected to be twice as long as a single hat section’s perimeter (developed width).

Considering the computational complexity and cost of coupled nonlinear FE simulations, reduced-order or surrogate models are used to approximate the responses defined in Eq. (17). Different metamodeling techniques have been developed and reported in the literature. Although they vary in terms of complexity and accuracy, they all rely on measured responses at the selected training points in the design space to find the unknown coefficients of a specific metamodel such that the approximation error is less than an acceptable threshold.

Radial basis functions (RBF) have been shown (Fang et al., 2005; Wang and Shan, 2007; Acar and Rais-Rohani, 2008; Parrish et al., 2012) to be suitable for approximating highly nonlinear responses using relatively small number of training points. In RBF formulation (Fang and Horstemeyer, 2006), the approximate response f̠y at a design point defined by the normalized design variable vector y is found as


where yi is the vector of normalized design variables at the ith training point, ||y yi|| = riis the Euclidian norm or distance from design point yto the ith training point, and M is the total number of functions included in the summation. The λi parameters are the unknown interpolation coefficients with φ representing the radially symmetric basis function that can take different forms. We considered both the multiquadric ϳr=r2+c2and Gaussian ϳr=exp(-cr2)basis functions, where c is a tuning parameter that can vary in the range of 0 < c ≤ 1 depending on the selected response.

In the coupled process-performance FE simulations, it is necessary to perform the deep-drawing, springback, and crush simulations in sequence. However, once a stand-alone surrogate model for each response is built, all the responses can be evaluated simultaneously, which provides considerable computational cost savings in the design optimization analysis. Using the Latin hypercube sampling (LHS) to produce a uniform distribution of design points, fifty training points were generated in the seven design-variable (7-dimensional) design space. Table 2 lists the training points and the corresponding values of the selected design variables. The maximum and minimum values for each design variable are shown in bold. Ten additional random design points are also generated as test points to measure the accuracy of the surrogate models.

Six responses are extracted for each set of simulations at each training point with the calculated response values listed in Table 3, where the value selected as the target for each response is shown in bold.

Two error metrics are considered. For the cross-validation normalized root-mean-square error (NRMSE) estimation (Lin et al., 1999), a metamodel is created using all except one training point and then the predicted response at the omitted point is compared to the corresponding true response value to measure the approximation error. This process is repeated for all the training points and the average is used as the overall error of the metamodel. NRMSE is calculated using

Pointx1 (mm)x2 (mm)x3 (mm)x4
x6 (m/s)x7Pointx1 (mm)x2 (mm)x3 (mm)x4
x6 (m/s)x7

Table 2.

Selected training points and design variable values in the process-product design space


where K is the number of training points, fiis the actual response obtained from the coupled FE simulations (expect for mass), and fi̠ is the response predicted by the model that excludes the contribution of the ith point with fmax and fmin as the maximum and minimum values of the response, respectively.

The second error metric is obtained by fitting a metamodel based on all the training points and taking the absolute value of the difference between the approximate and the true response values at each of the selected test points with the average error of all the test points as the final indicator of metamodel accuracy.

PointRTSPmax (kN)Pm

Table 3.

Results from sequential coupled process-performance simulations at the training points

The metamodels were tuned by selecting the parameter c and the RBF formulation that produced the least error for each response. Table 4 shows the RBF tuning parameter and formulation used for each response and the corresponding NRMSE and average error estimates. It is seen that the average errors for the test points also validate the surrogate models created for the optimization problem. In order to enhance the accuracy of the metamodels for thinning and spring back responses, the actual responses are transformed using a logarithmic function as presented in Table 4.

Error (%)
Rupture, RGNone1.0002.29.2
Thinning, TGln(T)0.0012.02.4
Springback, SMln(1x107 S)0.5001.83.5
Max Crush Force, PmaxMNone0.5004.52.0
Mean Crush Force, PmMNone0.1003.75.6
Mass, MGNone0.0104.61.8

Table 4.

Metamodel type, parameter, and approximation error for each response


6. Optimization results and discussion

The surrogate-based design optimization problem in Eq. (17) is solved using the multi-objective genetic algorithm (MOGA) toolbox in MATLAB (Fonseca and Fleming, 1993) with a randomly generated initial population size of 105 representing different combinations of design variable values within the specified bounds in Eq. (17). The subsequent generations are populated using the tournament selection algorithm with a crossover fraction of 80% using intermediate crossover function and the termination function tolerance of 1e-4. Stopping criterion is set at generation number 1400. The specific steps taken in the application of MOGA to this problem are as follows:

  1. Design variables expressed in real number are converted into bit strings.

  2. A random initial population is generated.

  3. Using a fitness function, members of the population are examined by

  4. assigning a rank to each solution based on non-dominated front (Sun et al., 2010).

  5. assigning a fitness value based on Pareto ranking.

  6. calculating the niche count of each solution.

  7. calculating the shared fitness value of each solution.

  8. normalizing the fitness values by using share fitness values.

  9. Using a stochastic method to select parents for the next generation.

  10. Performing crossover and mutation operations.

  11. Establishing a new population.

  12. Evaluating the population attributes.

  13. Continuing steps 3 to 7 to evaluate all the objectives.

  14. Selecting half of the individuals that have a higher rank than the rest.

  15. Continuing the solution process until a stopping criterion is satisfied based on the average change in the spread of Pareto solution being less than the tolerance specified.

The solution to the optimization problem in Eq. (17) for the Pareto optimum set converged after 139 GA generations and 14,701 function calls. The forty-five design points forming the Pareto frontier are listed in Table 5 in particular order. For ease of detection, the Pareto ID number and the corresponding best value for each objective are highlighted in the table. The results show that Pareto ID nos. 3, 4, 5, 6, 7, and 9 have the best values for rupture, springback, thinning, max crush force, mass, and mean crush force, respectively. A closer examination of the response values in Table 5 reveals that the response values at two or more design points may be very close to each other or nearly equal. This is because of the nonlinear distribution of the points on the Pareto frontier. For all the objectives, the preferred value is the smallest one except for the mean crush force. This is because the tube’s energy absorption capacity improves by increasing the mean crush force value. However, it is preferred to reduce the max crush force for such components.

The range of variation over the Pareto frontier is different from one objective to another. Specifically, the relative variations from the best to the worst values are approximately 757%, 20%, 1354%, 133%, -59%, and 145% for rupture, thinning, springback, max crush force, mean crush force, and mass, respectively. The most significant range of variation in

Design VariablesResponses
Pareto IDx1 (mm)x2 (mm)x3 (mm)x4
x7 RTSPmax (kN)Pm
4552.6234.494.552.2631.666.560.27 3337.918.850.69177.964.20.24

Table 5.

Design points on the Pareto optimum frontier

the springback reveals the high sensitivity of this response, relative to the rest, to changes in the design variable values. In comparison, thinning seems to be affected the least by the changes in design. Such information helps in identifying the critical responses for both process and product design considerations.

Given the six-dimensional space of the process-product criteria space, the process and performance objectives are plotted separately and shown in Fig. 11. In addition, a sample subset of the Pareto set is selected with the individual tube geometries and crush deformation modes shown in Fig. 12. Among the six points shown, design points 6, 11, and 33 are shown among the best choices to minimize the mass and maximum crush force.

The results indicate that the Pareto set consists mostly of a tube design that is larger in total height than width with approximately 72% and 67% having larger thickness and longer corner radius than the respective average values, respectively. The general trend appears to be toward a tube design model with dissimilar width and height dimensions, which can be traced to two contributing factors: (1) while only a portion of width is work hardened, the entire height section undergoes plastic deformation during the forming process; and (2) the flanges (short tabs) in the double-hat geometry (see Fig. 9) influence how the different sides of the tube deform and contribute to the crush energy absorption. In most cases, it appears to be preferable for the holding force to be less than its average value in approximately 60% of the Pareto set. A nearly equal percentage prefers a lower friction coefficient while a slightly higher portion (roughly 67%) prefers a higher punch velocity than the respective average values.

Figure 11.

Distribution of Pareto optimal set in performance (a) and process (b) criteria subspaces

To measure the approximation error in the optimum process and performance objective values, a complete verification simulation was performed on seven samples with the relative errors shown in Table 6. Most error values are fairly low, less than 5%, with the highest reaching 14.7%, which is reasonable for the types of simulation involved and the approximation techniques used in solving the design optimization problem.

Figure 12.

Selected design points from the Pareto set

Pareto IDR

Table 6.

Relative Error in Responses at Selected Pareto Points


7. Conclusion

A methodology for concurrent process-product design optimization using coupled nonlinear finite-element (FE) simulations was presented and applied to a sheet formed component made of AZ31 magnesium alloy. All FE simulations were performed using the Abaqus Explicit and Standard solvers. Surrogate models based on radial basis functions were developed for process and performance response approximations to facilitate the numerical multi-objective design optimization process.

The results of this investigation lead to the following conclusions:

  • Material and component geometry variations can be modeled using a sequence of coupled nonlinear FE simulations with careful transfer of state variables and other information from one simulation stage to the next.

  • Both the manufacturing and geometric design variables can have significant influence on the energy absorption behavior of the formed tube considered.

  • All process responses (i.e., rupture, thinning, and springback) were greatly influenced by the initial blank thickness value and corner radius.

  • The results of the multi-objective optimization problem highlighted different levels of conflict among the process and performance objectives considered. Moreover, the variation of objectives over the Pareto frontier indicated differing levels of sensitivity to changes in the design properties with springback being the most noteworthy.



This material is based on the work supported by the US Department of Energy under Award Number DE-EE0002323. Authors also acknowledge the collaboration provided through the SIMULIA Research & Development program under which licenses of Abaqus were provided.


  1. 1. AcarE.Rais-RohaniM.2008Ensemble of Metamodels with Optimized Weight Factors. Structural and Multidisciplinary Optimization, 373279294
  2. 2. AcarE.Rais-RohaniM.NajafiA.MarinE. B.BammannD.2009Concurrent Design of Product-Material Systems Using Multilevel Optimization. Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Palm Springs, CA, USA, May 4-7, 2009
  3. 3. BayramY. B.NiedH. F.2000Enriched Finite Element-Penalty Function Method for Modeling Interface Cracks With Contact. Engineering Fracture Mechanics, 65541557
  4. 4. BottcherC. S.FrikS.2003Consideration of Manufacturing Effects to Improve Crash Simulation Accuracy, 4th European LS-Dyna Users Conference, 2003
  5. 5. ChengJ. H.KikuchiN.1985An Analysis of Metal Forming Processes Using Large Deformation Elastic-Plastic Formulations. Computer Methods in Applied Mechanics and Engineering, 4971108
  6. 6. ChungW. J.ChoJ. W.BelytschkoT.1998On the Dynamic Effects of Explicit FEM in Sheet Metal Forming Analysis. Engineering Computations, 15750776
  7. 7. DagsonN.2001Influence of the Forming Process on the Crash Response of a Roof Rail Component, Master Thesis, Department of Solid Mechanics, Linkping University, March 2001.
  8. 8. FangH.HorstemeyerM. F.2006Global Response Approximation with Radial Basis Functions. Engineering Optimization, 384407424
  9. 9. FangH.Rais-RohaniM.LiuZ.HorstmeyerM. F.2005A Comparative Study of Metamodeling Methods for Multi-objective Crashworthiness Optimization. Computers and Structures, 8321212136
  10. 10. FonsecaM.FlemingP. J.1993Genetic Algorithms for Multiobjective Optimization: Formulation, Discussion And Generalization, Proceedings of the 5th International Conference on Genetic Algorithms, San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 416423
  11. 11. GrantabR.2006Interaction between Forming and Crashworthiness of Advanced High Strength Steel S-Rails, PhD dissertation, University of Waterloo, Canada, 2006
  12. 12. HuW.YaoL. G.HuaZ. Z.2008Optimization of Sheet Metal Forming Processes by Adaptive Response Surface Based on Intelligent Sampling Method. Journal of Materials Processing Technology, 1971-37788
  13. 13. KaufmanM.GainesD.KundrickK.LiuS. D.1998Integration of Chassis Frame Forming Analysis Into Performance Models to More Accurately Evaluate Crashworthiness. SAE Proceedings, paper 980551
  14. 14. KellicutA.CowellB.KavikondalaK.DuttonT.IregbuS.SturtR.1999Application of the Results of Forming Simulation in Crash Models. Proceedings of Numisheet 99 Conference, Besancon, France, 509514
  15. 15. KonakA.CoitD. W.SmithA. E.2006Multi-Objective Optimization Using Genetic Algorithms: A Tutorial. Reliability Engineering & System Safety, 919Special Issue- Genetic Algorithms and Reliability, 9921007
  16. 16. KrusperA.2003Influences of the Forming Process on the Crash Performance- Finite Element Analysis, Master Thesis, Department of Structural Mechanics, Chalmers University of Technology, May 2003
  17. 17. LeeS.KwonY. N.KangS. H.KimS. W.LeeJ. H.2008Forming Limit of AZ31 Alloy Sheet and Strain Rate on Warm Sheet Metal Forming, Journal of Materials Processing Technology, 2011-3431435
  18. 18. LiK. P.CardenW. P.WagonerR. H.2002Simulation of Springback, International Journal of Mechanical Sciences, 44103122
  19. 19. LinY.KrishnapurK.AllenJ. K.MistreeF.1999Robust Design: Goal Formulations and a Comparison of Metamodeling Method, Proceedings of DETC-99, 1999 ASME Design Engineering Technical Conferences, Las Vegas, ND, September 1215
  20. 20. LublinerJ.2006Plasticity Theory, Revised Edition
  21. 21. MayerR.2004Theoretical Effects of Hydroforming on Crashworthiness of Straight Sections. ASME Applied Mechanics Division, 255
  22. 22. Mc DowellD. L.ChoiH. J.PanchalJ. H.AustinR.AllenJ. K.MistreeF.2007Plasticity-Related Microstructure-Property Relations for Materials Design. Key Engineering Materials, 340-3412130
  23. 23. NajafiA.2011Coupled Sequential Process-Performance Simulation and Multi-Attribute Optimization of Structural Components Considering Manufacturing Effects, PhD Dissertation, Computational Engineering, Mississippi State University, 2011.
  24. 24. NajafiA.MarinE. B.Rais-RohaniM.2012Concurrent Multi-Scale Crush Simulations with a Crystal Plasticity Model. Thin-Walled Structures, 53176187
  25. 25. NoelsL.StainierL.PonthotJ. P.2004Combined Implicit/Explicit Time Integration Algorithms for the Numerical Simulation of Sheet Metal Forming. Journal of Computational and Applied Mathematics, 168331339
  26. 26. OdenJ. T.KikuchiN.SongY. J.1982Penalty-Finite Element Methods for the Analysis of Stokesian Flows. Computer Methods in Applied Mechanics and Engineering, 31297326
  27. 27. OdenJ. T.PiresE. B.1983Numerical Analysis of Certain Contact Problems in Elasticity With Non-Classical Friction Laws. Computers & Structures, 16481485
  28. 28. OliveiraD.WorswickM.GrantabR.WilliamsB.MayerR.2006Effect of Forming Process Variables on the Crashworthiness of Aluminium Alloy Tubes. International Journal of Impact Engineering, 325826846
  29. 29. OlsonG. B.1997Computational Design of Hierarchically Structured Materials. Science, 277533012371242
  30. 30. ParrishA.Rais-RohaniM.NajafiA.2012Crashworthiness Optimization of Vehicle Structures with Magnesium Alloy Parts, International Journal of Crashworthiness, 173259281
  31. 31. RebelG.ParkK. C.FelippaC. A.2002A Contact Formulation Based on Localized Lagrange Multipliers: Formulation and Application to Two Dimensional Problems. International Journal for Numerical Methods in Engineering, 54263297
  32. 32. RyouH.ChungK.YoonJ.HanC.2005Incorporation of Sheet-Forming Effects in Crash Simulations Using Ideal Forming Theory and Hybrid Membrane and Shell Method. Journal of Manufacturing Science and Engineering, 1271182192
  33. 33. SimoJ. C.HughesT. J. R.2000Computational Inelasticity, Springer, New York
  34. 34. SimunovicS.AramayoG.2002Steel Processing Properties and Their Effect on Impact Deformation of Lightweight Structures. Computer Science and Mathematics Division, Computational Materials Science Group, American Iron and Steel Institute
  35. 35. SouzaNeto. E. A.PericD.OwenD. R. J.2008Computational Methods for Plasticity: Theory and Applications, John Wiley & Sons
  36. 36. SunG.LiG.GongZ.CuiX.YangX.LiQ.2010Multiobjective Robust Optimization Method for Drawbead Design in Sheet Metal Forming. Materials & Design, 31419171929
  37. 37. van denBoogaard. A. H.MeindersT.HuétinkJ.2003Efficient Implicit Finite Element Analysis of Sheet Forming Processes, International Journal for Numerical Methods in Engineering, 56810831107
  38. 38. Wang-MN.BudianskyB.1978Analysis of Sheet Metal Stamping by Finite Element Method. Journal of Applied Mechanics, Transaction ASME, 457382
  39. 39. WangG.ShanS.2007Review of Metamodeling Techniques in Support of Engineering Design Optimization, Journal of Mechanical Design, 1294370380
  40. 40. WeiL.YuyingY.2008Multi-Objective Optimization of Sheet Metal Forming Process Using Pareto-Based Genetic Algorithm. Journal of Materials Processing Technology, 2081-3499506
  41. 41. WilliamsB. W.SimhaC. H. M.AbedrabboN.MayerR.WorswickM. J.2010Effect of Anisotropy, Kinematic Hardening, and Strain-Rate Sensitivity on the Predicted Axial Crush Response of Hydroformed Aluminum Alloy Tubes. International Journal of Impact Engineering, 376652661
  42. 42. WilliamsB. W.OliveiraD. A.WorswickM. J.MayerR.2005Crashworthiness of High and Low Pressure Hydroformed Straight Section Aluminum Tubes. Proceedings of SAE World Congress, paper 2005-01
  43. 43. WriggersP.2006Computational Contact Mechanics, Springer; 2nd edition

Written By

Ali Najafi and Masoud Rais-Rohani

Submitted: February 14th, 2012 Published: October 10th, 2012