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

Ali Najafi and Masoud Rais-Rohani

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/47852

## **1. Introduction**

354 Finite Element Analysis – Applications in Mechanical Engineering

machining. Int J Mach Tool Manu. 46: 782-800.

Milling. Int J Mach Tool Manu. 35: 751-760.

modelling. J Manuf Sci E-T Asme. 120: 623-631.

Engineering. 119: 502-508.

1 and 2. 3243: 558-564.

Vib Acoust. 130.

and Engineering. 125: 645-655.

56: 581-604.

[38] M. R. Miller, G. Mulholland and C. Anderson (2003) Experimental cutting tool

[39] Y. K. Potdar and A. T. Zehnder (2003) Measurements and Simulations of Temperature and Deformation Fields in Transient Metal Cutting. Journal of Manufacturing Science

[40] N. A. Abukhshim, P. T. Mativenga and M. A. Sheikh (2006) Heat generation and temperature prediction in metal cutting: A review and implications for high speed

[41] M. A. Davies, T. Ueda, R. M'Saoubi, B. Mullany and A. L. Cooke (2007) On the measurement of temperature in material removal processes. Cirp Ann-Manuf Techn.

[42] J. Lin (1995) Inverse estimation of the tool-work interface temperature in end milling.

[43] P. Kwon, T. Schiemann and R. Kountanya (2001) An inverse estimation scheme to measure steady-state tool–chip interface temperatures using an infrared camera.

[44] G. Fang and P. Zeng (2005) Three-dimensional thermo-elastic-plastic coupled FEM simulations for metal oblique cutting processes. J Mater Process Tech. 168: 42-48. [45] J. M. Lin (1995) Inverse Estimation of the Tool-Work Interface Temperature in End

[46] M. Chen, F. H. Sun, H. L. Wang, R. W. Yuan, Z. H. Qu and S. Q. Zhang (2003) Experimental research on the dynamic characteristics of the cutting temperature in the

[47] C. E. Leshock and Y. C. Shin (1997) Investigation on Cutting Temperature in Turning by a Tool-Work Thermocouple Technique. Journal of Manufacturing Science and

[48] S. Fraser, M. H. Attia and M. O. M. Osman (1998) Modelling, identification and control of thermal deformation of machine tool structures, part 1: Concept of generalized

[49] X. Zhang, L. Zhu, D. Zhang, H. Ding, Y. Xiong, Numerical robust optimization of spindle speed for milling process with uncertainties, International Journal of Machine Tools and Manufacture (2012), http://dx.doi.org/10.1016/j.ijmachtools.2012.05.002 [50] T. K. Hasselman and J. D. Chrostowski (1997) Effects of product and experimental variability on model verification of automobile structures. P Soc Photo-Opt Ins. 3089: 612-620. [51] E. Balmes (1998) Predicted variability and differences between tests of a single structure. Imac - Proceedings of the 16th International Modal Analysis Conference, Vols

[52] M. H. Kurdi, R. T. Haftka, T. L. Schmitz and B. P. Mann (2008) A Robust Semi-Analytical Method for Calculating the Response Sensitivity of a Time Delay System. J

[54] M. H. Richardson and D. L. Formenti, Parameter estimation from frequency response

[55] Y. Altintas, E. Shamoto, P. Lee and E. Budak (1999) Analytical prediction of stability

[53] J. Nocedal and S. J. Wright (1999) Numerical optimization: Springer verlag.

measurements using rational fraction polynomials, 1982.

lobes in ball end milling. J Manuf Sci E-T Asme. 121: 586-592.

International Journal of Machine Tools and Manufacture. 35: 751-760.

International Journal of Machine Tools and Manufacture. 41: 1015-1030.

process of high-speed milling. J Mater Process Tech. 138: 468-471.

temperature distributions. J Manuf Sci E-T Asme. 125: 667-673.

Finite-element (FE) simulations are often used to predict the response characteristics of a structural component under different boundary conditions and to help explore the design space for the optimum design while minimizing the need for physical testing. It has also been used to model various manufacturing processes, especially those involving the forming process (Cheng and Kikuchi 1985; Chung et al. 1998; Li et al., 2002; van den Boogaard et al., 2003).

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

© 2012 Najafi and Rais-Rohani, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Najafi and Rais-Rohani, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

integrated material-process-product design (Olson, 1997; McDowell et al., 2007; Acar et al., 2009).

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

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

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.

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.

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

The constitutive relation of the rate dependent classical plasticity model with isotropic

�� = �� (�� � ��

��(�� �)

��) (1)

�� � �� = �� (3)

�(�� �) = |�| � ��� � ��� � � (2)

investigate the capabilities of the model in coupled process-performance simulations.

**2. Modeling of manufacturing effects** 

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 hydroformingcrush 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**

356 Finite Element Analysis – Applications in Mechanical Engineering

the manufacturing process on product design variables.

the difference in the overall response was relatively modest.

2009).

integrated material-process-product design (Olson, 1997; McDowell et al., 2007; Acar et al.,

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

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 hydroformingcrush 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

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 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

$$
\dot{\sigma} = \mathbb{C} \colon \left( \dot{\varepsilon} - \dot{\varepsilon}^{\nu p} \right) \tag{1}
$$

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

$$f(\sigma, \kappa) = |\sigma| - \left(\sigma\_{\mathcal{Y}} + H\kappa\right) \le 0 \tag{2}$$

Flow rule and hardening law

$$
\varepsilon^{\nu p} = \dot{\lambda} \, \frac{\partial f(\sigma, \kappa)}{\partial \sigma}, \dot{\kappa} = \dot{\lambda} \tag{3}
$$

$$\dot{\lambda}\{\sigma,\kappa\} = \begin{cases} \frac{1}{c} \left[ \left( \frac{|\sigma|}{\sigma\_{\mathcal{V}}} \right)^{1/p} - 1 \right] \text{ if } f(\sigma,\kappa) \ge 0\\ 0 \text{ if } f(\sigma,\kappa) < 0 \end{cases} \tag{4}$$

$$\mathbf{e}\_{n+1} = \mathbf{e}\_{n+1} - \frac{1}{3} \operatorname{tr}(\mathbf{e}\_{n+1}) I \tag{5}$$

$$\mathbf{s}\_{n+1}^t = 2\mu \{ \mathbf{e}\_{n+1} - \mathbf{e}\_n^p \} \tag{6}$$

$$(\blacksquare)\_{n+1} = (\blacksquare)\_{n+1}^{\mathbb{C}} \tag{7}$$

$$f\_{n+1}^{\mathbf{f}} \approx f\_n^{\mathbf{f}} + \frac{\partial f}{\partial \Delta \lambda} \Delta \lambda\_{n+1} \tag{8}$$

$$
\kappa\_{n+1} = \kappa\_n + \sqrt{\frac{2}{3}} \,\Delta\lambda\_{n+1} \tag{9}
$$

$$
\varepsilon\_{n+1}^{\upsilon p} = \varepsilon\_n^{\upsilon p} + \Delta \lambda\_{n+1} \mathfrak{n}\_{n+1} \tag{10}
$$

$$
\Delta \sigma\_{n+1} = \overline{\lambda} \operatorname{tr} \{ \Delta \varepsilon\_{n+1} \} I + \Delta \mathbf{s}\_{n+1}^t - 2\mu \, \Delta \lambda\_{n+1} \, \mathfrak{n}\_{n+1} \tag{11}
$$

$$
\sigma\_{n+1} = \sigma\_n + \Delta \sigma\_{n+1} \tag{12}
$$

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

�� � � �(�� � )

� � �(�� � )

> � and ��

(13)

(14)

� are the principal major

Rupture and thinning are two responses that can be calculated from the results of the deepdrawing 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

> � � �(�� � ))�

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

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 �� and �� 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

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

� �(�� � �� �� )�

� ���

to behave linearly in both compressive and tensile plastic strains as shown in Fig. 3.

extracted from the forming limit diagram (FLD) using the following equation

� �(��

���

�

�

���� ��

�

���

where �(��) is the equation representing FLD curve and ��

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

**3.2. Springback simulation** 

� ����

�

���

constant for all the elements, �� is the shell thickness assigned to the elements.

�

**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 deepdrawing 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\{ \sum\_{l=1}^{I} R\_l^{\; 2} = \sum\_{l=1}^{I} (\varepsilon\_1^l - \phi(\varepsilon\_2^l))^2 \qquad \varepsilon\_1^l > \phi(\varepsilon\_2^l) \right. \tag{13}$$

$$\varepsilon\_1^l \le \phi(\varepsilon\_2^l)$$

where �(��) is the equation representing FLD curve and �� � and �� � 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 = \sum\_{l=1}^{N} T\_l^2 = \sum\_{l=1}^{N} (\frac{t\_l - t\_o}{t\_o})^2 \tag{14}$$

where �� and �� 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, �� is the shell thickness assigned to the elements.

#### **3.2. Springback simulation**

360 Finite Element Analysis – Applications in Mechanical Engineering

**Figure 1.** Sheet-forming FE model of a single hat section

*uy* = 0

Die

simulation

*ux* = 0

*uy* = 0

**Figure 2.** Formed hat section with the boundary conditions used for the subsequent springback

*uz* = 0 *ux* = 0

Blank

Holder

applied to preserve the constant gripping force throughout the drawing process.

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

*uz* = 0

*uy* = 0

Punch

Holder

Die

*uy* = 0

*x* 

*z* 

*y* 

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 = Max\left(\sqrt{(X\_l - X\_o)^2 + (Y\_l - Y\_o)^2 + (Z\_l - Z\_o)^2}\right) \quad \forall l = 1, NN \tag{15}$$

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

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

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

Tube Rigid Wall

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

Applied Displacement

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

**3.4. Crush simulation** 

stages discussed previously.

**Figure 5.** General setup for crush simulation

fine mesh density that results in relatively thick shell elements.

*δeff*

where �� �� � 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**

362 Finite Element Analysis – Applications in Mechanical Engineering

the mesh.

springback simulations.

**3.3. Joining and trimming** 

welding process is not modeled here.

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

where �� �� � 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

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

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

� � �����(�����)� + (�����)� + (�����)��������� � �� �� (15)

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

$$P\_m = \frac{1}{\delta\_{eff}} \int\_0^t F(t)D(t)dt\tag{16}$$

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

�� = 4300 s��

�� =1s��

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

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

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-tosurface contact between rigid wall and tube are specified. A classical multi-linear kinematic

**3.6. Effect of manufacturing process on product performance** 

hardening material model is used for this comparison.

captured by finding the plastic strain.

where �(�) is the instantaneous contact force normal to the rigid wall surface and �(�) is the instantaneous cross-head axial displacement of the rigid wall. The final cross-head displacement of the rigid wall represents ���� and it is assumed to be 125 mm, or 50% of the tube length. The maximum crush force is the largest value of �(�) 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 inhouse 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.

364 Finite Element Analysis – Applications in Mechanical Engineering

equation for finding the mean crush force is expressed as

house FORTRAN code to calculate the mean crush force values.

(Abaqus/Explicit)

1. Deep-Drawing Simulation

**Figure 6.** Sequence of coupled nonlinear FE simulations

induced in the tube model in step 3.

**3.5. Material model** 

�� <sup>=</sup> <sup>1</sup> ����

under the crush force-displacement curve by the effective crush distance. The mathematical

�

�

(Abaqus/Explicit) 2. Springback Simulation

(Abaqus/Standard)

3. Joining &

where �(�) is the instantaneous contact force normal to the rigid wall surface and �(�) is the instantaneous cross-head axial displacement of the rigid wall. The final cross-head displacement of the rigid wall represents ���� and it is assumed to be 125 mm, or 50% of the tube length. The maximum crush force is the largest value of �(�) 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-

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

Trimming 4. Crush Simulation

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

� �(�)�(�)��

(16)

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-tosurface contact between rigid wall and tube are specified. A classical multi-linear kinematic hardening material model is used for this comparison.

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

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.

Corner Radius (mm)

Thicknes s (mm)

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

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,

Holding Force (kN)

Punch Velocity (m/s)

Friction Coefficien t

Height (mm)

**Table 1.** The values assigned to design variables for sensitivity analysis

Design Variable Width (mm)

**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.

366 Finite Element Analysis – Applications in Mechanical Engineering

kN to 18 kN when the history effects are considered.

**Figure 9.** Description of geometric design variables for a 250-mm long tube

Height

Corner Radius

**4. Sensitivity analysis** 

/ die set friction coefficients.

**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

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

Width

Thickness

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


**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.

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

Decrease (-15%) Increase (+15%)

Thickness Corner Radius

Height Width

**Figure 10.** General sensitivity of rupture (a), thinning (b), springback (c), maximum crush force (d), and

(e)

Holding Force

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 processperformance optimization problem as considered here, both manufacturing- and performancelevel attributes are optimized. When faced with competing objectives, the resulting multiobjective optimization problem becomes one of finding not just one but a collection of nondominated design points that form the Pareto frontier. By specifying a particular target value

� (�(�) � ��)�, ( �(�) � ��)�, (�(�) � ��)�,

20 35 mm 2.5 7.5 mm 1.0 2.5 mm 10 50 kN 2 10 m / s 0.1 0.35

s.t. 40 70 mm

*x x x x x x x*

where design variables *x*1 to *x*<sup>7</sup> represent width, height, corner radius, thickness, holding force, punch velocity, and friction coefficient, respectively with each having both lowerand 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).

 

 (��(�) � ���)�, (����(�) � ���)�, (�(��,�,��) � ��)��

(17)

for each objective, the multi-objective optimization problem is expressed as

mean crush force (e) to the process and product design variables

Friction Coefficient

Punch Velocity

**5. Multi-objective design optimization** 


**Mean Crush Force**

min ����,�,��

**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**

368 Finite Element Analysis – Applications in Mechanical Engineering




> -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2

**Maximum Crush Force**

Friction Coefficient

Punch Velocity Holding Force

**Springback**

Friction Coefficient

**Thinning**

**Rupture**

Friction Coefficient

Friction Coefficient

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

Decrease (-15%) Increase (+15%)

seem to be greater than those of the manufacturing process parameters.

Punch Velocity

Punch Velocity

Punch Velocity Holding Force

Holding Force

Holding Force

Thickness Corner Radius

Decrease (-15%) Increase (+15%)

(a)

Thickness Corner Radius

Decrease (-15%) Increase (+15%)

(b)

Thickness Corner Radius

Decrease (-15%) Increase (+15%)

(c)

Thickness Corner Radius

(d)

Height Width

Height Width

Height Width

Height Width

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 processperformance optimization problem as considered here, both manufacturing- and performancelevel attributes are optimized. When faced with competing objectives, the resulting multiobjective optimization problem becomes one of finding not just one but a collection of nondominated 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\_{\mathbf{x}}\min\_{\mathbf{x},\mathbf{x}\_{1},\ldots,\mathbf{x}\_{7}}\left\{ (P\_{\mathbf{m}}(\mathbf{x}) - TR)^{2}, (T(\mathbf{x}) - TT)^{2}, (S(\mathbf{x}) - TS)^{2}, \ldots, \mathbf{x}\_{7} \right\}$$

$$\text{s.t.} \quad 40 \le \mathbf{x}\_{1} \le 70 \text{ mm}$$

$$\begin{aligned} \text{s.t.} \quad 40 &\le \mathbf{x}\_{1} \le 70 \text{ mm} \\ 20 &\le \mathbf{x}\_{2} \le 35 \text{ mm} \\ 2.5 &\le \mathbf{x}\_{3} \le 7.5 \text{ mm} \\ 1.0 &\le \mathbf{x}\_{4} \le 2.5 \text{ mm} \\ 10 &\le \mathbf{x}\_{5} \le 50 \text{ kN} \\ 2 &\le \mathbf{x}\_{6} \le 10 \text{ m}/\text{s} \\ 0.1 &\le \mathbf{x}\_{7} \le 0.35 \end{aligned} \tag{17}$$

where design variables *x*1 to *x*<sup>7</sup> represent width, height, corner radius, thickness, holding force, punch velocity, and friction coefficient, respectively with each having both lowerand 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 ��(�) at a design point defined by the normalized design variable vector *y* is found as

$$\hat{f}(\mathbf{y}) = \sum\_{l=1}^{M} \lambda\_l \phi(\|\mathbf{y} - \mathbf{y}\_l\|) \tag{18}$$

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

(mm)

*x*2 (mm)

*x*3 (mm)

*x*4 (mm)

�(���� <sup>−</sup> ����) (19)

� is the response predicted by the model that excludes

*x*5 (kN)

*x*6 (m/s) *x*<sup>7</sup>

Point *x*<sup>1</sup> (mm)

*x*2 (mm)

*x*3 (mm)

*x*4 (mm)

*x*5 (kN) *x*6

(m/s) *<sup>x</sup>*7 Point *<sup>x</sup>*<sup>1</sup>

**Table 2.** Selected training points and design variable values in the process-product design space

����� <sup>−</sup> ��

�

���

�� �

where *K* is the number of training points, �� is the actual response obtained from the coupled

the contribution of the ith point with ���� and ���� as the maximum and minimum values of

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

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

����� � �<sup>1</sup>

FE simulations (expect for mass), and ��

as the final indicator of metamodel accuracy.

the response, respectively.

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

where *yi* is the vector of normalized design variables at the *ith* training point, ����� ����� � = � �� is the Euclidian norm or distance from design point � to the *i*th 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 �(�) = √�� � �� and Gaussian �(�) = ���(����) 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 deepdrawing, 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 designvariable (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

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

normalized design variable vector *y* is found as

of 0 < *c* ≤ 1 depending on the selected response.

response is shown in bold.

metamodel. NRMSE is calculated using

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 ��(�) at a design point defined by the

> ��(�) =����(‖����‖ �

where *yi* is the vector of normalized design variables at the *ith* training point, ����� ����� � = � �� is the Euclidian norm or distance from design point � to the *i*th 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 �(�) = √�� � �� and Gaussian �(�) = ���(����) basis functions, where *c* is a tuning parameter that can vary in the range

In the coupled process-performance FE simulations, it is necessary to perform the deepdrawing, 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 designvariable (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

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

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

generated as test points to measure the accuracy of the surrogate models.

) (18)

���


**Table 2.** Selected training points and design variable values in the process-product design space

$$NRMSE = \sqrt{\frac{1}{K} \sum\_{l=1}^{K} \left(f\_l - \widehat{f\_l}\right)^2} \left/ \left(f\_{max} - f\_{mln}\right)\right. \tag{19}$$

where *K* is the number of training points, �� is the actual response obtained from the coupled FE simulations (expect for mass), and �� � is the response predicted by the model that excludes the contribution of the ith point with ���� and ���� 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.

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

372 Finite Element Analysis – Applications in Mechanical Engineering


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

The surrogate-based design optimization problem in Eq. (17) is solved using the multiobjective 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

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

10. 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

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

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

 normalizing the fitness values by using share fitness values. 4. Using a stochastic method to select parents for the next generation.

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

 assigning a fitness value based on Pareto ranking. calculating the niche count of each solution.

calculating the shared fitness value of each solution.

**6. Optimization results and discussion** 

application of MOGA to this problem are as follows:

5. Performing crossover and mutation operations.

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

6. Establishing a new population. 7. Evaluating the population attributes.

2. A random initial population is generated.

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.


aG = Gaussian, M = Multiquadric

**Table 4.** Metamodel type, parameter, and approximation error for each response

**Table 3.** Results from sequential coupled process-performance simulations at the training points

## **6. Optimization results and discussion**

372 Finite Element Analysis – Applications in Mechanical Engineering

(kN)

*Pm* (kN) *M*

**Table 3.** Results from sequential coupled process-performance simulations at the training points

using a logarithmic function as presented in Table 4.

Typea

**Table 4.** Metamodel type, parameter, and approximation error for each response

Response RBF

aG = Gaussian, M = Multiquadric

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

Response

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(1x107 *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

Transformation *<sup>c</sup>*

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

(kg) Point *<sup>R</sup> <sup>T</sup> <sup>S</sup> Pmax*

(kN)

NRMSE (%)

Average Error (%)

*Pm* (kN)

*M* (kg)

Point *R T S Pmax*

The surrogate-based design optimization problem in Eq. (17) is solved using the multiobjective 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:

	- assigning a rank to each solution based on non-dominated front (Sun et al., 2010).
	- assigning a fitness value based on Pareto ranking.
	- calculating the niche count of each solution.
	- calculating the shared fitness value of each solution.
	- normalizing the fitness values by using share fitness values.

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


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

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

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

**Figure 11.** Distribution of Pareto optimal set in performance (a) and process (b) criteria subspaces

approximation techniques used in solving the design optimization problem.

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

(a) (b)

process and product design considerations.

average values.


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.

374 Finite Element Analysis – Applications in Mechanical Engineering

**Table 5.** Design points on the Pareto optimum frontier

*x*4 (mm)

*x*3 (mm)

Pareto ID

*x*1 (mm)

*x*2 (mm)

Design Variables Responses

*x*6

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

(m/s) *<sup>x</sup>*<sup>7</sup> *<sup>R</sup> <sup>T</sup> <sup>S</sup> Pmax*

(kN)

*Pm*  (kN)

*M* (kg)

*x*5 (kN)

> 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.

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

Both the manufacturing and geometric design variables can have significant influence

All process responses (i.e., rupture, thinning, and springback) were greatly influenced

 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

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

Acar, E. & Rais-Rohani, M. (2008). Ensemble of Metamodels with Optimized Weight Factors.

Acar, E.; Rais-Rohani, M., Najafi, A., Marin, E. B. & Bammann, D. (2009). Concurrent Design of Product-Material Systems Using Multilevel Optimization. *Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference*,

Bayram, Y.B. & Nied, H.F. (2000). Enriched Finite Element-Penalty Function Method for Modeling Interface Cracks With Contact. *Engineering Fracture Mechanics*, Vol.65, pp.541-

Bottcher, C.S. & Frik, S. (2003). Consideration of Manufacturing Effects to Improve Crash

Cheng, J.H. & Kikuchi, N. (1985). An Analysis of Metal Forming Processes Using Large Deformation Elastic-Plastic Formulations. *Computer Methods in Applied Mechanics and* 

Chung, W.J.; Cho, J.W. & Belytschko, T. (1998). On the Dynamic Effects of Explicit FEM in

Sheet Metal Forming Analysis. *Engineering Computations*, Vol.15, pp. 750–776 Dagson, N. (2001). Inuence of the Forming Process on the Crash Response of a Roof Rail Component, Master Thesis, Department of Solid Mechanics, Linkping University,

*Structural and Multidisciplinary Optimization*, Vol.37, No.3, pp.279-294

Simulation Accuracy, 4th European LS-Dyna Users Conference, 2003

changes in the design properties with springback being the most noteworthy.

on the energy absorption behavior of the formed tube considered.

by the initial blank thickness value and corner radius.

**Author details** 

**Acknowledgement** 

were provided.

**8. References** 

557

March 2001.

Ali Najafi and Masoud Rais-Rohani *Mississippi State University, USA* 

Palm Springs, CA, USA, May 4-7, 2009

*Engineering*, Vol.49, pp.71-108

**Figure 12.** Selected design points from the Pareto set


**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.


## **Author details**

376 Finite Element Analysis – Applications in Mechanical Engineering

**Figure 12.** Selected design points from the Pareto set

*T* (%)

**Table 6.** Relative Error in Responses at Selected Pareto Points

The results of this investigation lead to the following conclusions:

information from one simulation stage to the next.

*S* (%)

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

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

 Material and component geometry variations can be modeled using a sequence of coupled nonlinear FE simulations with careful transfer of state variables and other

*Pmax* (%)

*Pm* (%)

*M* (%)

*R* (%)

Pareto ID

**7. Conclusion** 

design optimization process.

Ali Najafi and Masoud Rais-Rohani *Mississippi State University, USA* 

## **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.

## **8. References**


Fang, H. & Horstemeyer, M. F. (2006). Global Response Approximation with Radial Basis Functions. *Engineering Optimization*, Vol.38, No.4, pp.407–424

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

Najafi, A.; Marin, E.B. & Rais-Rohani, M. (2012). Concurrent Multi-Scale Crush Simulations

Noels, L.; Stainier, L. & Ponthot, J.P. (2004). Combined Implicit/Explicit Time Integration Algorithms for the Numerical Simulation of Sheet Metal Forming. *Journal of* 

Oden, J.T.; Kikuchi, N. & Song, Y.J. (1982). Penalty-Finite Element Methods for the Analysis of Stokesian Flows. *Computer Methods in Applied Mechanics and Engineering*, Vol.31,

Oden, J.T. & Pires, E.B. (1983). Numerical Analysis of Certain Contact Problems in Elasticity With Non-Classical Friction Laws. *Computers & Structures*, Vol.16, pp.481-485 Oliveira, D.; Worswick, M., Grantab R., Williams B. & Mayer R. (2006). Effect of Forming Process Variables on the Crashworthiness of Aluminium Alloy Tubes. *International* 

Olson, G. B. (1997). Computational Design of Hierarchically Structured Materials. *Science*,

Parrish, A.; Rais-Rohani, M. & Najafi, A. (2012). Crashworthiness Optimization of Vehicle Structures with Magnesium Alloy Parts, *International Journal of Crashworthiness*, Vol.17,

Rebel, G.; Park, K. C. & Felippa, C.A. (2002). A Contact Formulation Based on Localized Lagrange Multipliers: Formulation and Application to Two Dimensional Problems.

Ryou, H.; Chung K., Yoon J. & Han C. (2005). Incorporation of Sheet-Forming Effects in Crash Simulations Using Ideal Forming Theory and Hybrid Membrane and Shell Method. *Journal of Manufacturing Science and Engineering*, Vol.127, No.1, pp. 182–192

Simunovic, S. & Aramayo G. (2002). Steel Processing Properties and Their Effect on Impact Deformation of Lightweight Structures. Computer Science and Mathematics Division,

Souza Neto, E. A.; Peric, D. & Owen, D.R.J. (2008). *Computational Methods for Plasticity:* 

Sun, G.; Li, G. Gong, Z., Cui, X., Yang, X. & Li, Q. (2010). Multiobjective Robust Optimization Method for Drawbead Design in Sheet Metal Forming. *Materials & Design*,

van den Boogaard, A.H.; Meinders, T. & Huétink, J. (2003). Efficient Implicit Finite Element Analysis of Sheet Forming Processes, *International Journal for Numerical Methods in* 

Wang, N-M. & Budiansky, B. (1978). Analysis of Sheet Metal Stamping by Finite Element

Wang, G. & Shan, S. (2007). Review of Metamodeling Techniques in Support of Engineering

Wei, L. & Yuying, Y. (2008). Multi-Objective Optimization of Sheet Metal Forming Process Using Pareto-Based Genetic Algorithm. *Journal of Materials Processing Technology*,

Method. *Journal of Applied Mechanics, Transaction ASME*, Vol.45, pp. 73–82

Design Optimization, *Journal of Mechanical Design*, Vol.129, No.4, pp.370–380

*International Journal for Numerical Methods in Engineering*, Vol.54, pp. 263-297

Simo, J. C. & Hughes, T. J. R. (2000). *Computational Inelasticity*, Springer, New York

Computational Materials Science Group, American Iron and Steel Institute

with a Crystal Plasticity Model. *Thin-Walled Structures*, Vol.53, pp. 176-187

*Computational and Applied Mathematics*, Vol.168, pp. 331-339

*Journal of Impact Engineering*, Vol.32, No. 5, pp. 826–846

pp.297-326

Vol.277, No.5330, pp. 1237-42

*Theory and Applications*, John Wiley & Sons

*Engineering*, Vol.56, No.8, pp. 1083–1107

Vol.31, Issue 4, pp.1917-1929

Vol.208, No.1-3, 21, pp. 499-506

Issue 3, pp. 259-281


Najafi, A.; Marin, E.B. & Rais-Rohani, M. (2012). Concurrent Multi-Scale Crush Simulations with a Crystal Plasticity Model. *Thin-Walled Structures*, Vol.53, pp. 176-187

378 Finite Element Analysis – Applications in Mechanical Engineering

*and Structures*, Vol.83, pp.2121-2136

*Processing Technology*, Vol.197, No.1-3, pp.77-88

*Conference*, Besancon, France, pp. 509–514, 1999

University of Technology, May 2003

*Technology*, Vol.201, No.1-3, pp.431-435

Lubliner J. (2006). *Plasticity Theory*, Revised Edition

*Engineering Materials*, Vol. 340-341, pp. 21-30

*Journal of Mechanical Sciences*, Vol.44, pp.103-122

Sections. *ASME Applied Mechanics Division*, Vol.255, 591-603

Crashworthiness. *SAE Proceedings*, paper no. 980551, 1998

Issue - Genetic Algorithms and Reliability, pp.992-1007

Inc., pp.416-423, 1993

Fang, H. & Horstemeyer, M. F. (2006). Global Response Approximation with Radial Basis

Fang, H.; Rais-Rohani, M., Liu, Z. & Horstmeyer, M.F. (2005). A Comparative Study of Metamodeling Methods for Multi-objective Crashworthiness Optimization. *Computers* 

Fonseca, M. & Fleming, P. J. (1993). Genetic Algorithms for Multiobjective Optimization: Formulation, Discussion And Generalization, *Proceedings of the 5th International Conference on Genetic Algorithms*, San Francisco, CA, USA: Morgan Kaufmann Publishers

Grantab, R. (2006). Interaction between Forming and Crashworthiness of Advanced High Strength Steel S-Rails, PhD dissertation, University of Waterloo, Canada, 2006 Hu, W.; Yao, L.G. & Hua, Z. Z. (2008). Optimization of Sheet Metal Forming Processes by Adaptive Response Surface Based on Intelligent Sampling Method. *Journal of Materials* 

Kaufman, M.; Gaines, D., Kundrick, K. & Liu, S. D. (1998). Integration of Chassis Frame Forming Analysis Into Performance Models to More Accurately Evaluate

Kellicut, A.; Cowell, B., Kavikondala, K., Dutton, T., Iregbu, S. & Sturt, R. (1999). Application of the Results of Forming Simulation in Crash Models. *Proceedings of Numisheet 99* 

Konak, A.; Coit, D.W. & Smith, A. E. (2006). Multi-Objective Optimization Using Genetic Algorithms: A Tutorial. *Reliability Engineering & System Safety*, Vol.91, No.9, Special

Krusper, A. (2003). Inuences of the Forming Process on the Crash Performance - Finite Element Analysis, Master Thesis, Department of Structural Mechanics, Chalmers

Lee, S.; Kwon, Y.N., Kang, S.H., Kim, S.W. & Lee, J.H. (2008). Forming Limit of AZ31 Alloy Sheet and Strain Rate on Warm Sheet Metal Forming, *Journal of Materials Processing* 

Li, K. P.; Carden, W. P. & Wagoner, R.H. (2002). Simulation of Springback, *International* 

Lin, Y.; Krishnapur, K., Allen, J.K. & F. Mistree. (1999). Robust Design: Goal Formulations and a Comparison of Metamodeling Method, *Proceedings of DETC-99, 1999 ASME* 

Mayer, R. (2004). Theoretical Effects of Hydroforming on Crashworthiness of Straight

McDowell, D.L.; Choi, H. J., Panchal, J.H., Austin, R., Allen, J.K. & Mistree, F. (2007). Plasticity-Related Microstructure-Property Relations for Materials Design. *Key* 

Najafi, A. (2011). Coupled Sequential Process-Performance Simulation and Multi-Attribute Optimization of Structural Components Considering Manufacturing Effects, PhD

Dissertation, Computational Engineering, Mississippi State University, 2011.

*Design Engineering Technical Conferences*, Las Vegas, ND, September 12-15, 1999.

Functions. *Engineering Optimization*, Vol.38, No.4, pp.407–424

	- Williams, B.W.; Simha, C.H.M., Abedrabbo, N. Mayer, R. & Worswick, M.J. (2010). Effect of Anisotropy, Kinematic Hardening, and Strain-Rate Sensitivity on the Predicted Axial Crush Response of Hydroformed Aluminum Alloy Tubes. *International Journal of Impact Engineering*, Vol.37, No.6, pp. 652-661
	- Williams, B.W.; Oliveira, D.A., Worswick, M.J. & Mayer, R. (2005) Crashworthiness of High and Low Pressure Hydroformed Straight Section Aluminum Tubes. *Proceedings of SAE World Congress*, paper No. 2005-01-0095, 2005
	- Wriggers, P. (2006). *Computational Contact Mechanics*, Springer; 2nd edition

*Engineering*, Vol.37, No.6, pp. 652-661

*World Congress*, paper No. 2005-01-0095, 2005

Wriggers, P. (2006). *Computational Contact Mechanics*, Springer; 2nd edition

Williams, B.W.; Simha, C.H.M., Abedrabbo, N. Mayer, R. & Worswick, M.J. (2010). Effect of Anisotropy, Kinematic Hardening, and Strain-Rate Sensitivity on the Predicted Axial Crush Response of Hydroformed Aluminum Alloy Tubes. *International Journal of Impact* 

Williams, B.W.; Oliveira, D.A., Worswick, M.J. & Mayer, R. (2005) Crashworthiness of High and Low Pressure Hydroformed Straight Section Aluminum Tubes. *Proceedings of SAE* 

## *Edited by Farzad Ebrahimi*

In the past few decades, the Finite Element Analysis (FEA) has been developed into a key indispensable technology in the modeling and simulation of various engineering systems. The present book is a result of contributions of experts from international scientific community and collects original and innovative research studies on recent applications of FEA in five major topics of mechanical engineering namely, fluid mechanics and heat transfer, machine elements analysis and design, machining and product design, wave propagation and failure-analysis and structural mechanics and composite materials. It is meant to provide a small but valuable sample of contemporary research activities around the world in this field and it is expected to be useful to a large number of researchers. The introductions, data, and references in this book will help the readers know more about this topic and help them explore this exciting and fast-evolving field.

Finite Element Analysis - Applications in Mechanical Engineering

Finite Element Analysis

Applications in Mechanical Engineering

*Edited by Farzad Ebrahimi*

Photo by OlgaSalt / iStock