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

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

$$\stackrel{\text{̥}}{ϡ}=\mathbb{C}:(\stackrel{\text{̥}}{ϓ}-{\stackrel{\text{̥}}{ϓ}}_{}^{\mathit{vp}})$$E1

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

$$f\left(ϡ,Ϙ\right)=\left|ϡ\right|-\left({ϡ}_y+\mathit{HϘ}\right)≣0$$E2

• Flow rule and hardening law

$${\stackrel{\text{̥}}{ϓ}}_{}^{\mathit{vp}}=\stackrel{\text{̥}}{ϙ}\mathrm{}\frac{\complement f\left(ϡ,Ϙ\right)}{\complement ϡ},\stackrel{\text{̥}}{Ϙ}=\stackrel{\text{̥}}{ϙ}$$E3

$$\stackrel{\text{̥}}{ϙ}\left(ϡ,Ϙ\right)=\left\{\begin{array}{c}\frac{1}{C}\left[{\left(\frac{\left|ϡ\right|}{{ϡ}_y}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$P$}\right.}-1\right]\mathrm{}\mathit{if}\mathrm{}f\left(ϡ,Ϙ\right)\le 0\\ \mathrm{}0\mathrm{}\mathit{if}\mathrm{}f\left(ϡ,Ϙ\right)<0\end{array}\mathrm{}\right.$$E4

where the bold symbols $ϡ$and $ϓ$represent the second-rank stress and stain tensors, respectively, and $\mathbb{C}$is the forth-rank tensor of isotropic elasticity. The superscript vp designates the viscoplastic component. The $\stackrel{\text{̥}}{∍}$ terms indicate time derivatives. The yield stress is denoted by ${ϡ}_y$ with parameters H and ҝ representing isotropic hardening modulus and plastic multiplier, respectively. Moreover, $C$and $P$in 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$▯\left(∍\right)$. The deviatoric part of total strain tensor $e$is updated as

$$e_{n+1}^{}={ϓ}_{n+1}^{}-\frac{1}{3}\mathit{tr}\left({ϓ}_{n+1}^{}\right)I$$E5

where $I$is the identity tensor and $\mathit{tr}\left(\mathrm{∍}\right)$is the trace of the second rank tensor. The trial stress tensor $s_{n+1}^t\mathrm{}$is calculated from the Hook’s law given by

$$s_{n+1}^t=2Ϛ\left(e_{n+1}^{}-e_n^p\right)$$E6

where the $e_n^p$ is the plastic strain at the nth increment and $Ϛ$is the shear modulus. Then, the calculated stress magnitude value at n+1 increment $‗s_{n+1}^t‗\mathrm{}$is evaluated with the yield surface $f_{n+1}^t$ at the same increment. If $f_{n+1}^t=‗s_{n+1}^t‗-\sqrt{\frac{2}{3}}\left[{ϡ}_y+\mathit{HϘ}\right]≣0$ is met, the material is in the elastic region and the trial values are the admissible stress values and, therefore, the plastic multiplier $\stackrel{\text{̥}}{ϙ}$ or ${\mathrm{\beta }ϙ}_{n+1}$ is zero with

$${\left(\mathrm{∍}\right)}_{n+1}^{}={\left(\mathrm{∍}\right)}_{n+1}^t$$E7

And if $f_{n+1}^t>0,\mathrm{}$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

$$f_{n+1}^t\ncong f_n^t+\frac{\complement f}{\complement \varnothing \mathrm{ϙ}}{\mathrm{\beta }ϙ}_{n+1}$$E8

where ${\mathrm{\beta }ϙ}_{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 $\mathrm{\beta }\left(∍\right)$to the previous state of stress as

$${{\mathrm{Ϙ}}_{n+1}={\mathrm{Ϙ}}_n+\sqrt{\frac{2}{3}}\mathrm{\beta }ϙ}_{n+1}$$E9

$${ϓ}_{n+1}^{\mathit{vp}}={ϓ}_n^{\mathit{vp}}+{\mathrm{\beta }ϙ}_{n+1}\mathrm{}{n}_{n+1}^{}$$E10

$$\mathrm{\beta }{ϡ}_{n+1}=\overline{ϙ}\mathrm{}\mathit{tr}\left({\mathrm{\beta }ϓ}_{n+1}^{}\right)I+{\mathrm{\beta }s}_{n+1}^t-2Ϛ{\mathrm{\beta }ϙ}_{n+1}\mathrm{}{n}_{n+1}^{}$$E11

$${ϡ}_{n+1}={ϡ}_n+\mathrm{\beta }{ϡ}_{n+1}$$E12

where $\mathit{tr}\left(∍\right)$represents the trace of the tensor, trial values are denoted with t superscript, and the normal to the yield surface is found as$n_{n+1}^{}=\mathrm{}\raisebox{1ex}{$s_{n+1}^t$}\!\left/ \!\raisebox{-1ex}{$‗s_{n+1}^t‗$}\right.$. By having${\mathrm{Ϙ}}_n={ϓ}_n^{\mathit{vp}}$, the manufacturing effect can be considered through a non-zero state of ${ϡ}_n$ and ${ϓ}_n^{\mathit{vp}}$ at the beginning of the analysis as ${ϡ}_1$ and${ϓ}_1^{\mathit{vp}}$.

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=\left\{\begin{array}{c}\coprod _i^{=}{R_i}^2=\coprod _i^{=}{({ϓ}_1^i-ϳ({ϓ}_2^i\left)\right)}^2\mathrm{}{ϓ}_1^i>ϳ\left({ϓ}_2^i\right)\mathrm{}\\ 0\mathrm{}{ϓ}_1^i≣ϳ\left({ϓ}_2^i\right)\mathrm{}\end{array}\right.$$E13

where $ϳ\left({ϓ}_2\right)$is the equation representing FLD curve and ${ϓ}_1^i$ and ${ϓ}_2^i$ 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

$$T=\coprod _i^{=}{T_i}^2=\coprod _i^{=}{\left(\frac{t_i-t_o}{t_o}\right)}^2$$E14

where $t_o$ and $t_i$ 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, $t_o$is 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=\mathit{Max}\left(\sqrt{{\left(X_i{-X}_o\right)}^2+{\left(Y_i{-Y}_o\right)}^2+{\left(Z_i{-Z}_o\right)}^2}\mathrm{}\right)\mathrm{}⇿i=1,\mathit{NN}$$E15

where $X,Y,Z$are 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, P_max while the mean crush force, P_m 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

$$P_m=\frac{1}{{\Upsilon }_{\mathit{eff}}}{\cup }_0^{t}F\left(t\right)D\left(t\right)\mathit{dt}\mathrm{}$$E16

where $F\left(t\right)$is the instantaneous contact force normal to the rigid wall surface and $D\left(t\right)$is the instantaneous cross-head axial displacement of the rigid wall. The final cross-head displacement of the rigid wall represents ${\Upsilon }_{\mathit{eff}}$ and it is assumed to be 125 mm, or 50% of the tube length. The maximum crush force is the largest value of $F\left(t\right)$found 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/m³, 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.

Design; Variable	Width (mm)	Height (mm)	Corner Radius (mm)	Thickness (mm)	Holding Force (kN)	Punch Velocity (m/s)	Friction Coefficient
Width	41.25	27.5	5.0	1.75	30	6.0	0.225
±15%	63.25	27.5	5.0	1.75	30	6.0	0.225
Height	55	20.625	5.0	1.75	30	6.0	0.225
±15%	55	31.625	5.0	1.75	30	6.0	0.225
Corner	55	27.5	3.75	1.75	30	6.0	0.225
Radius ±15%	55	27.5	5.75	1.75	30	6.0	0.225
Thickness	55	27.5	5.0	1.3125	30	6.0	0.225
±15%	55	27.5	5.0	2.0125	30	6.0	0.225
Holding	55	27.5	5.0	1.75	22.5	6.0	0.225
Force ±15%	55	27.5	5.0	1.75	34.5	6.0	0.225
Punch	55	27.5	5.0	1.75	30	4.5	0.225
Velocity ±15%	55	27.5	5.0	1.75	30	6.9	0.225
Friction	55	27.5	5.0	1.75	30	6.0	0.16875
Coefficient ±15%	55	27.5	5.0	1.75	30	6.0	0.25875
Upper Bound	70	35	7.5	2.5	50	10.0	0.35
Mean	55	27.5	5.0	1.75	30	6.0	0.225
Lower Bound	40	20	2.5	1.0	100	2.0	0.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

$$\begin{array}{c}\min \\ x=x_1,‧,x_7\end{array}\begin{array}{c}\min \\ x=x_1,‧,x_7\end{array}\left\{\begin{array}{c}{\left(R\left(x\right)-\mathit{TR}\right)}^2,{\left(\mathrm{}T\left(x\right)-\mathit{TT}\right)}^2,\mathrm{}{\left(S\left(x\right)-\mathit{TS}\right)}^2,\\ \mathrm{}{{(P}_m\left(x\right)-\mathit{TPM})}^2,{\left(P_{\mathit{max}}\left(x\right)-\mathit{TPX}\right)}^2,\mathrm{}{(M\left(x_1,‧,x_4\right)-\mathit{TM})}^2\end{array}\right\}$$E17

$$\begin{array}{l}\text{s}.\text{t}.\text{4}0\le x_{\text{1}}\le \text{ 7}0\text{ mm}\\ \text{2}0\le x_{\text{2}}\le \text{ 35 mm}\\ \text{2}.\text{5 }\le x_{\text{3}}\le \text{ 7}.\text{5 mm}\\ \text{1}.0\le x_{\text{4}}\le \text{ 2}.\text{5 mm}\\ \text{1}0\le x_{\text{5}}\le \text{ 5}0\text{ kN}\\ \text{2 }\le x_{\text{6}}\le \text{ 1}0\text{ m}/\text{s}\\ 0.\text{1 }\le x_{\text{7}}\le 0.\text{35}\end{array}$$E18

where design variables x₁ to x₇ 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, P_m(x) is the mean crush force, P_max(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 $\stackrel{\text{̠}}{f}\left(y\right)$ at a design point defined by the normalized design variable vector y is found as

$$\stackrel{\text{̠}}{f}\left(y\right)=\coprod _i^{=}{ϙ}_iϳ(‗y-y_i‗)$$E19

where y_i is the vector of normalized design variables at the ith training point, $\left|\right|y—yi\left|\right|=r_i$is the Euclidian norm or distance from design point $y$to 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 $ϳ\left(r\right)=\sqrt{r^2+c^2}$and Gaussian $ϳ\left(r\right)=\exp (-c{r}^2)$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

Point	x1 (mm)	x2 (mm)	x3 (mm)	x4; (mm)	x5; (kN)	x6 (m/s)	x7	Point	x1 (mm)	x2 (mm)	x3 (mm)	x4; (mm)	x5; (kN)	x6 (m/s)	x7
1	60.44	27.19	4.91	1.26	31.93	5.86	0.18	26	52.98	31.62	4.94	1.34	21.35	4.90	0.29
2	40.00	35.14	4.58	1.60	36.68	7.76	0.23	27	70.88	21.22	5.42	2.19	23.66	6.47	0.27
3	53.14	25.15	5.72	2.27	31.18	6.26	0.22	28	64.74	25.38	5.59	1.63	31.58	4.29	0.26
4	56.48	27.61	4.86	1.58	29.09	7.02	0.25	29	61.85	34.06	5.98	1.29	22.27	6.81	0.16
5	62.84	34.42	3.74	1.42	30.85	6.70	0.18	30	47.60	21.91	6.12	2.24	32.31	6.32	0.24
6	44.85	24.15	6.41	1.84	38.00	4.69	0.25	31	39.66	32.03	3.91	1.93	30.46	7.09	0.21
7	56.30	32.39	5.13	1.91	29.42	6.60	0.19	32	46.59	21.26	3.73	1.74	23.37	4.71	0.26
8	42.46	23.00	5.10	2.14	26.58	5.77	0.17	33	57.54	30.83	4.58	1.82	25.82	6.76	0.24
9	63.58	26.23	4.68	2.19	38.47	4.62	0.22	34	43.90	29.55	5.52	1.62	26.97	6.92	0.20
10	54.23	33.49	3.99	1.28	38.96	6.56	0.20	35	65.50	31.16	4.81	1.55	24.56	7.45	0.24
11	48.30	29.98	6.03	2.00	36.15	6.05	0.29	36	66.23	35.09	6.16	1.23	34.19	5.82	0.28
12	43.44	26.91	4.38	2.11	34.52	4.42	0.24	37	50.61	23.55	3.68	1.51	28.36	4.79	0.16
13	42.84	29.06	5.19	1.78	32.73	6.13	0.22	38	54.81	19.89	6.50	1.66	21.68	5.26	0.18
14	67.81	20.25	5.27	1.33	33.95	5.38	0.21	39	38.59	32.69	5.76	2.15	33.21	5.56	0.27
15	46.24	24.20	6.35	1.37	28.02	7.30	0.23	40	50.36	26.68	5.65	1.76	35.66	5.99	0.17
16	69.79	20.21	4.43	1.53	27.13	4.25	0.27	41	55.16	21.57	3.82	1.88	24.72	5.67	0.19
17	51.31	33.00	4.31	2.03	25.58	4.95	0.19	42	48.47	20.71	3.52	1.70	35.24	6.40	0.21
18	70.76	19.49	4.48	1.71	35.85	7.37	0.17	43	61.19	22.46	5.43	2.08	26.22	7.70	0.28
19	52.31	28.09	5.35	2.05	37.52	5.03	0.27	44	67.06	23.50	3.59	2.02	22.06	5.15	0.26
20	41.12	22.63	4.72	1.46	37.64	5.47	0.25	45	58.65	35.70	5.87	1.44	30.33	5.10	0.16
21	59.51	28.37	4.15	1.98	23.96	5.60	0.21	46	49.22	25.76	4.07	1.39	37.08	7.53	0.21
22	57.95	29.43	5.81	1.49	22.61	7.18	0.28	47	68.86	28.68	4.17	1.95	34.92	7.59	0.19
23	65.67	30.75	4.26	2.09	25.26	6.21	0.23	48	41.67	33.31	6.23	1.68	27.49	5.29	0.20
24	45.56	25.88	6.28	1.80	29.67	7.26	0.28	49	60.22	24.72	3.96	1.41	33.33	4.49	0.17
25	64.11	34.45	5.05	2.23	23.08	4.55	0.23	50	69.41	30.27	5.92	1.86	28.69	6.94	0.28

Table 2.

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

$$\mathit{NRMSE}=\raisebox{1ex}{$\sqrt{\frac{1}{K}\coprod _i^{=}{\left(f_i-\stackrel{\text{̠}}{f_i}\right)}^2}$}\!\left/ \!\raisebox{-1ex}{$(f_{\mathit{max}}-f_{\mathit{min}})$}\right.$$E20

where K is the number of training points, $f_i$is the actual response obtained from the coupled FE simulations (expect for mass), and $\stackrel{\text{̠}}{f_i}$ is the response predicted by the model that excludes the contribution of the ith point with $f_{\mathit{max}}$ and $f_{\mathit{min}}$ 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.

Point	R	T	S	Pmax (kN)	Pm; (kN)	M; (kg)	Point	R	T	S	Pmax; (kN)	Pm; (kN)	M; (kg)
1	538.3	36.6	0.21	95.5	30.8	0.13	26	1003.5	88.4	1.83	101.0	34.7	0.14
2	1911.8	173.2	0.29	116.3	43.3	0.16	27	1672.9	109.0	0.31	160.9	62.1	0.22
3	1938.2	119.0	0.68	160.2	67.4	0.21	28	881.1	62.5	1.72	123.9	41.6	0.17
4	1370.0	116.8	0.79	116.8	43.1	0.16	29	455.3	32.0	2.56	109.6	33.8	0.15
5	1477.6	141.8	3.60	123.2	35.7	0.16	30	1521.7	81.5	1.16	144.6	63.7	0.18
6	870.2	45.3	2.10	114.8	54.0	0.15	31	3328.9	298.8	1.18	136.8	57.0	0.18
7	1697.4	139.7	0.34	154.0	60.0	0.20	32	2123.5	165.4	0.40	102.7	36.2	0.14
8	1774.8	107.1	0.92	132.0	54.3	0.17	33	2073.2	184.0	0.44	142.8	51.7	0.19
9	2098.0	161.2	0.50	167.0	66.3	0.22	34	1112.4	85.7	1.78	112.4	54.9	0.15
10	1253.6	117.3	3.54	100.6	30.1	0.14	35	1282.8	111.4	1.48	132.0	42.5	0.17
11	1663.9	123.8	1.36	145.5	54.7	0.19	36	633.2	59.6	1.63	108.4	33.9	0.15
12	2470.5	186.3	0.50	140.6	53.6	0.18	37	1427.3	112.1	1.03	96.6	28.8	0.13
13	1651.3	124.6	1.62	121.3	50.0	0.16	38	476.2	15.9	1.48	103.1	42.2	0.14
14	538.4	32.1	0.62	92.3	30.0	0.13	39	2113.2	154.7	0.85	153.6	64.9	0.20
15	573.7	34.9	1.16	86.4	31.6	0.12	40	1003.8	66.3	1.81	122.2	46.3	0.16
16	1115.1	78.7	0.59	108.5	37.0	0.15	41	2158.6	168.0	0.83	121.7	45.1	0.16
17	2118.4	184.5	1.77	158.6	50.3	0.21	42	2313.8	180.7	0.50	100.3	34.9	0.13
18	1288.4	92.9	0.32	120.6	39.6	0.17	43	1779.5	122.4	0.42	145.2	60.6	0.20
19	1971.0	141.4	0.77	148.5	57.9	0.20	44	3037.0	251.1	0.70	149.3	49.0	0.20
20	1148.4	81.1	1.34	84.0	29.9	0.11	45	528.6	36.1	1.04	123.7	45.6	0.17
21	2188.3	184.8	0.94	151.8	51.8	0.20	46	1497.3	133.9	1.21	93.0	30.5	0.12
22	938.5	76.0	1.68	115.4	43.1	0.15	47	2242.6	195.4	0.61	162.3	59.7	0.22
23	2675.3	228.7	1.16	175.1	54.0	0.23	48	818.8	57.2	2.95	124.5	47.1	0.16
24	1213.2	82.4	1.66	119.2	45.0	0.16	49	1017.6	78.1	0.66	100.6	30.9	0.14
25	1940.7	171.3	1.422	195.1	68.1	0.26	50	1265.0	103.4	0.81	158.0	53.4	0.21

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.

Response	RBF; Type^a	Response; Transformation	c	NRMSE; (%)	Average; Error (%)
Rupture, R	G	None	1.000	2.2	9.2
Thinning, T	G	ln(T)	0.001	2.0	2.4
Springback, S	M	ln(1x10^7 S)	0.500	1.8	3.5
Max Crush Force, Pmax	M	None	0.500	4.5	2.0
Mean Crush Force, Pm	M	None	0.100	3.7	5.6
Mass, M	G	None	0.010	4.6	1.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 Variables									Responses
Pareto ID	x1 (mm)		x2 (mm)	x3 (mm)	x4; (mm)	x5; (kN)	x6; (m/s)	x7		R		T	S	Pmax (kN)	Pm; (kN)	M; (kg)
1	43.35		27.93	4.86	1.28	25.06	6.95	0.20		848.0		17.28	2.25	86.0	36.2	0.11
2	48.98		30.77	6.36	1.55	22.50	6.24	0.23		624.2		17.04	2.32	112.5	46.2	0.15
3	40.44		31.85	5.55	1.31	26.93	5.12	0.18		455.4		16.65	3.16	95.9	38.8	0.12
4	58.36		33.16	5.35	2.09	32.40	6.85	0.23		2217.0		18.38	0.24	171.7	64.0	0.23
5	65.41		22.66	6.28	1.28	26.72	6.08	0.22		684.0		15.93	1.52	93.7	32.3	0.13
6	51.72		24.55	6.09	1.30	27.16	5.76	0.21		474.8		16.26	1.63	84.0	30.3	0.12
7	54.92		22.86	5.55	1.28	25.34	6.73	0.21		594.6		16.54	1.09	85.1	29.1	0.11
8	65.80		35.05	6.32	2.26	32.40	6.68	0.20		1543.5		18.16	1.21	195.4	71.4	0.27
9	63.30		35.61	6.44	2.26	32.41	6.68	0.20		1533.6		18.19	1.14	193.7	71.6	0.27
10	58.99		34.75	6.13	2.24	31.87	6.57	0.23		1924.1		18.27	0.61	186.9	69.5	0.25
11	40.05		31.91	6.33	1.65	28.18	5.00	0.18		569.7		16.78	3.39	119.3	47.8	0.15
12	43.97		31.87	3.88	2.25	29.99	6.46	0.27		3902.1		18.99	1.21	162.3	62.0	0.21
13	41.42		31.79	5.98	1.85	28.20	5.22	0.19		1143.4		17.41	2.98	133.9	53.4	0.17
14	44.79		31.86	3.74	2.13	29.97	6.42	0.27		3751.6		19.06	1.23	153.8	57.7	0.20
15	39.44		31.91	6.46	1.35	28.18	5.00	0.18		762.9		16.08	3.49	100.3	41.5	0.12
16	44.16		31.91	3.75	1.88	28.41	6.52	0.23		3130.2		18.91	1.88	134.9	51.3	0.18
17	40.52		31.89	4.45	1.99	28.33	5.08	0.28		2700.7		18.75	1.57	141.0	54.5	0.18
18	42.16		32.76	4.24	2.25	30.28	6.22	0.22		3606.3		18.83	1.76	164.0	63.7	0.21
19	55.93		33.29	5.12	2.25	31.21	6.57	0.24		2764.3		18.55	0.44	181.7	67.0	0.24
20	57.84		30.86	5.35	2.08	28.88	6.61	0.22		2046.6		18.21	0.51	166.7	63.6	0.22
21	40.67		31.85	5.88	1.63	27.44	5.15	0.19		811.9		17.13	3.09	117.3	46.9	0.15
22	41.36		31.88	6.13	1.68	28.18	5.16	0.18		724.3		17.01	3.20	121.7	48.4	0.16
23	48.48		31.92	3.91	2.16	30.45	6.46	0.27		3685.4		19.00	1.03	160.3	58.6	0.21
24	63.50		34.16	5.71	2.26	32.02	6.65	0.20		1983.8		18.29	0.78	192.6	70.3	0.26
25	40.06		31.96	6.34	1.51	29.18	5.57	0.22		503.7		16.94	2.88	108.1	44.2	0.14
26	57.74		33.56	5.22	2.17	32.15	6.61	0.23		2484.5		18.48	0.34	178.5	65.7	0.24
27	49.24		32.26	5.00	2.25	30.16	6.60	0.26		3062.1		18.60	0.80	170.3	64.7	0.23
28	40.60		32.47	3.78	2.19	28.72	6.34	0.26		3871.6		19.04	1.52	156.8	61.3	0.20
29	54.19		32.99	5.00	2.12	31.53	6.66	0.25		2706.8		18.57	0.42	168.1	62.1	0.22
30	42.42		31.88	5.24	2.02	29.14	6.09	0.19		2195.0		18.19	2.24	146.3	59.3	0.19
31	39.63		31.85	6.15	1.81	28.19	5.01	0.18		951.1		17.18	3.27	130.6	52.5	0.17
32	42.26		33.15	4.45	2.24	31.63	6.62	0.19		3423.0		18.75	1.93	164.0	64.4	0.21
33	41.73		27.49	3.92	1.81	28.97	5.97	0.28		2648.1		18.71	1.45	119.0	44.8	0.16
34	44.99		31.29	3.75	2.17	30.27	6.12	0.26		3741.3		19.03	1.33	155.9	58.2	0.21
35	43.41		32.29	5.04	1.74	29.59	5.88	0.21		1699.2		18.02	2.18	125.7	49.7	0.17
36	40.14		31.90	3.97	1.91	29.01	6.02	0.21		2982.7		18.75	2.33	135.1	53.6	0.17
37	41.30		31.81	5.54	1.93	28.85	5.41	0.19		1651.7		17.82	2.65	139.1	55.7	0.18
38	43.90		32.73	4.56	2.08	31.34	6.27	0.23		3026.6		18.65	1.47	152.4	58.8	0.20
39	60.43		34.90	6.28	2.26	32.08	6.59	0.19		1630.5		18.13	0.95	190.4	71.3	0.26
40	40.71		32.24	4.97	1.43	29.79	4.99	0.23		1029.4		17.67	2.55	101.7	38.7	0.13
41	42.69		34.76	5.63	2.10	29.45	5.17	0.21		2003.1		18.10	2.15	158.7	60.3	0.21
42	61.52		32.37	5.47	2.22	30.16	6.52	0.23		2230.0		18.30	0.53	185.4	67.8	0.25
43	56.76		32.84	5.36	2.25	32.40	6.61	0.21		2405.2		18.37	0.49	182.7	69.0	0.24
44	57.81		33.58	5.24	2.25	32.20	6.58	0.23		2609.8		18.50	0.41	184.9	68.0	0.25
45	52.62		34.49	4.55	2.26	31.66	6.56	0.27		3337.9		18.85	0.69	177.9	64.2	0.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 ID	R; (%)	T; (%)	S; (%)	Pmax; (%)	Pm; (%)	M; (%)
4	7.8	1.6	4.7	1.1	0.4	0.4
6	5.7	3.8	4.3	3.1	12.2	0.3
11	12.2	1.0	6.1	2.5	14.7	0.7
18	5.6	0.8	1.0	1.9	5.3	0.1
23	9.1	0.7	4.3	1.1	1.0	0.1
33	0.9	2.9	1.9	1.3	2.9	0.9

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.

Acknowledgement

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.