Prediction of springback in sheet forming by a new finite strain model with nonlinear kinematic and isotropic hardening

Journal of Materials Processing Technology 209 (2009) 4062–4075

Contents lists available at ScienceDirect

Journal of Materials Processing Technology journal homepage: www.elsevier.com/locate/jmatprotec

Prediction of springback in sheet forming by a new finite strain model with nonlinear kinematic and isotropic hardening I.N. Vladimirov ∗ , M.P. Pietryga, S. Reese Institute of Solid Mechanics, Braunschweig University of Technology, D-38106 Braunschweig, Germany

a r t i c l e

i n f o

Article history: Received 9 June 2008 Received in revised form 11 September 2008 Accepted 18 September 2008 Keywords: Springback Finite strains Kinematic hardening Isotropic hardening Exponential map Experiments

a b s t r a c t The paper discusses the application of a finite strain model to predict springback in sheet forming. The concept combines both isotropic and kinematic hardening and utilizes a new algorithm based on the exponential map. A special focus is put on the simulation of the draw-bend test of DP600 dual phase steel sheets and the comparison of the numerical results with experimental data. Further, we investigate the sensitivity of the results with respect to geometry parameters as the tool radius and the sheet thickness as well as the ratio between isotropic and kinematic hardening. In addition, the model is applied to the S-rail benchmark problem and the thermoforming of thermoplastic blends, where large elastic strains occur. The simulation results match the experiments very well. © 2008 Elsevier B.V. All rights reserved.

1. Introduction Springback denotes the property of sheet metal parts to change shape upon unloading after forming. This shape change is elastically driven since upon removal of the deformation load a stress reduction takes place and the total strain is decreased by the amount of elastic strain which results in springback. The phenomenon of springback is a growing concern in the sheet metal forming industry because it leads to major assembly problems and complicates the design of the die. Thus, prediction and compensation of springback are essential in the design stage of forming processes to achieve precise final part shapes to avoid problems in the assembly. Due to the fact that traditional trial-and-error methods for springback compensation are time-consuming, while empirical adjustments are not applicable to complex geometries, the simulation of springback has become essential for the tool and process design. There are three major approaches in analyzing springback: the analytical method, the semi-analytical method, and the finite element analysis. The classical analytical approach is based upon simplifications of tool descriptions and material behaviour. Early analytical solutions for springback were established for plane stress pure bending with perfect elastoplasticity (Gardiner, 1957),

∗ Corresponding author. E-mail address: [email protected] (I.N. Vladimirov). 0924-0136/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jmatprotec.2008.09.027

and were extended to plane strain pure bending (Queener and De Angelis, 1968). Further developments include solutions for plane strain bending with superimposed tensile force (Wang et al., 1993) and for plane stress bending of integrated circuit leadframes (Chan and Wang, 1999). More recently, Jeunechamps et al. (2006) presented a closed form technique to predict springback in creep age-forming and investigated the influence of geometric parameters on springback of aluminium plates. Moon et al. (2008) developed an analytical model for the prediction of sidewall curl during the draw-bend test, based on the moment-curvature relationship. In general, the analytical approach is not very effective in the springback analysis of real forming processes since it assumes simplified process conditions and material properties. On the other hand, semi-analytical methods based upon the hybrid membrane/shell method are being used for the springback analysis of draw bending. Pourboghrat and Chu (1995) described a method to predict springback for plane strain stretch/draw operation. Further, a semi-analytical approach for springback prediction of an anisotropic aluminium alloy was developed by Pourboghrat et al. (1998). Recently, Lee et al. (2007) proposed a numerical procedure to evaluate springback in a 2D draw-bend test based on the hybrid method which superposes bending effects onto membrane solutions. With the fast increase in computational power and solution techniques, the finite element method (FEM) is increasingly being used for the prediction and simulation of springback. The major

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

advantage of finite element methods over analytical and semianalytical approaches is the possibility of modelling complicated tool descriptions and realistic material behaviour. Springback simulations are more sensitive to numerical tolerances than forming simulations, including effects of element type, integration scheme, and unloading scheme. Li et al. (2002) investigated the springback of sheets formed by the draw-bend test including parametric studies on process and numerical parameters such as back force, tool radius-to-sheet thickness ratio, friction, mesh size, number of integration points. They also observed that the presence of the Bauschinger effect significantly influences springback and thus must be taken into account. The Bauschinger effect is characterized by a reduced yield stress upon reverse loading after plastic deformation has occurred during the initial loading. The presence of the Bauschinger effect has also been addressed by Zhao and Lee (2001) and the use of isotropic hardening alone has been identified as a source of inaccurate springback prediction. In the work of Papeleux and Ponthot (2002) a maximum difference of about 7◦ in the springback angle predicted by isotropic and kinematic hardening models is reported, exclusive use of isotropic hardening leading to an overestimation of springback. Chung et al. (2005) evaluated the springback of automotive sheets based on a modification of the Chaboche combined isotropic/kinematic hardening model (Chaboche, 1989) and the nonquadratic anisotropic yield function of (Barlat et al., 2003). The formulation of the numerical model is described in (Chung et al., 2005), while the identification of the material parameters is discussed in Lee et al. (2005a), and the verification with experiments is presented in Lee et al. (2005b). Further, Ragai et al. (2005) assessed the influence of sheet anisotropy on the springback of draw-bend specimens both by means of experiments and finite element simulations. Accurate prediction of springback by either analytical, semianalytical or finite element approaches depends on the internal stress distribution within the sheet metal. The stress response calculated by the material model being used significantly affects the computed springback. When sheet metal workpieces are unloaded and removed from the tooling after forming, they experience elastic unloading and springback. During the reverse loading the material usually exhibits the Bauschinger effect which is represented by a translation of the yield surface in the stress space. It becomes clear that use of isotropic hardening alone, represented by an expansion of the yield surface, cannot properly simulate the Bauschinger effect and the springback. Thus, accurate simulation of springback requires the use of an appropriate constitutive model capable of describing the cyclic hardening behaviour and the Bauschinger effect. Large strain kinematic/isotropic hardening models applied to describe springback are discussed in a number of recent works. Yoshida and Uemori (2003) developed a two-surface combined hardening model based on the assumption of small elastic and large plastic strains. A hypoelastic–plastic model taking into account isotropic, kinematic and rotational hardening, to describe the rotation of the anisotropy axes, is introduced in Choi et al. (2006a) and applied to forming simulations in Choi et al. (2006b). Recently, Oliveira et al. (2007) studied the influence of the work-hardening modelling on the springback prediction by using the Swift power law, the Voce type saturation law, the Chaboche combined hardening model and Teodosiu’s microstructural model (Haddadi et al., 2006). In the work of Svendsen et al. (2006) an anisotropic large deformation model with combined hardening is applied to predict the springback in the draw-bend test. In a recent work the authors (Vladimirov et al., 2008) develop a finite strain model with combined nonlinear kinematic/isotopic hardening. For the integration of the evolution equations a new algorithm based on the exponential map has been introduced. This time integration scheme has the

4063

advantage that both the symmetry of the internal variables as well as plastic incompressibility are preserved. The comparison with two modified backward Euler schemes demonstrates higher accuracy and better performance in favour of the exponential algorithm, in particular for large time steps. The goal of the present paper is to investigate the applicability of the combined hardening model of Vladimirov et al. (2008) for the accurate prediction of springback. The paper is organized as follows. In Section 2 the derivation of the constitutive equations of the model is briefly summarized. Section 3 includes a short summary of the important features of the exponential algorithm used for the integration of the evolution equations for the internal variables. In Section 4 we investigate the springback prediction capabilities of the model for sheet metal specimens after forming. First, the draw-bend test of DP600 dual phase steel sheets is simulated and the results of various simulations are compared to experimental data. Herein, different tool radius-to-sheet thickness ratios are considered and the influence of isotropic, kinematic and combined hardening is examined. A further validation of the model is carried out by means of the S-rail benchmark problem of NUMISHEET96 (Lee et al., 1996). Finally, the model that incorporates both large elastic and large plastic strains is applied to the simulation of the thermoforming process of thermoplastic blends which exhibit significant elastic strains before plastification. 2. Material model A thorough derivation of the constitutive equations of the model from a thermodynamic framework has been presented in a recent publication by the authors (see Vladimirov et al., 2008). To enable a better understanding of the model, we repeat the most important steps of the derivation in the following. The starting point is the classical multiplicative decomposition of the deformation gradient F = Fe Fp into elastic (Fe ) and plastic (Fp ) parts. By means of the introduction of the additional multiplicative split of the plastic part of the deformation gradient into elastic and inelastic parts, i.e. Fp = Fpe Fpi , a continuum mechanical extension of the rheological model of Armstrong–Frederick kinematic hardening can be achieved. The multiplicative split of the plastic deformation gradient is physically motivated. According to Lion (2000), the elastic part of the plastic deformation gradient Fpe results from dislocation-induced lattice rotations and stretches on the microscale, whereas Fpi relates to local plastic deformations stemming from inelastic slip on crystallographic slip systems. The Helmholtz free energy per unit undeformed volume is assumed to have the form =

e (Ce ) +

kin (Cpe ) +

iso ()

(1)

−1 is the elastic right Cauchy–Green where Ce = FTe Fe = F−T p CFp

−1 deformation tensor and Cpe = FTpe Fpe = F−T pi Cp Fpi the elastic part of the plastic right Cauchy–Green tensor. The first part e represents the macroscopic elastic material properties. The second term kin corresponds to the elastic energy stored in dislocation fields due to kinematic hardening. The third term iso describes the additional amount of stored energy due to isotropic hardening where  is the isotropic hardening variable (accumulated plastic strain). Utilizing the Clausius–Duhem form of the second law of thermodynamics − ˙ + S · (1/2) C˙ ≥ 0 and assuming that the free energy e is an isotropic function of the elastic right Cauchy–Green tensor Ce , and also that kin is an isotropic function of Cpe , one obtains a relation for the second Piola–Kirchhoff stress tensor

S = 2F−1 p

∂ e −T F ∂Ce p

(2)

4064

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075 Table 1 Material parameters

and evolution equations for the plastic deformation (flow rule) dp = ˙

∂ MD − D = ˙ ∂M ||MD − D ||

(3)

as well as for kinematic and isotropic hardening, respectively: dpi



b = ˙ MD , c kin

˙ =

2 ˙  3

(4)

Here, dp is the symmetric part of the plastic velocity gradient lp = F˙ p F−1 p . Analogously, dpi denotes the symmetric part of lpi := F˙ p F−1 . The yield function  is of von Mises type and reads p i

i



D

D

 = ||M −  || −

2 (y − R) 3

(5)

where R = −∂ iso /∂ is the “drag” stress and the superscript D denotes the deviator, i.e. AD = A − (1/3)(tr A)1. The quantities M = 2Ce (∂ e /∂Ce ) and Mkin = 2Cpe (∂ kin /∂Cpe ) represent symmetric Mandel-type stress tensors. The tensor  := 2Fpe (∂ kin /∂Cpe )FTpe is the back stress in the intermediate configuration and ˙ is the plastic multiplier. A dot between two second-order tensors describes their scalar product. Finally, the Kuhn–Tucker conditions ˙ ≥ 0,  ≤ ˙ 0,  = 0 complete the set of constitutive equations. The derivation of the material model is suitably carried out in the intermediate configuration. For the numerical implementation of the constitutive equations it is, however, more appropriate to work in the undeformed or reference configuration. Carrying out the corresponding pull-back of the tensors M, Mkin and  and utilizing the relation C˙ p = 2FTp dp Fp enables the alternative representation of the constitutive model in the reference configuration: • Stress tensors S = 2F−1 p

∂ e −T F , ∂Ce p

X = 2 F−1 pi

∂ kin −T F , ∂Cpe pi

Y = CS − Cp X,

• Evolution equations

D

D T

C˙ pi

,

Y · (Y )

b = 2˙ YD Cp , c kin i

Shear modulus [MPa] Lamé constant [MPa] Initial yield stress [MPa] Rate of change of kinematic hardening modulus with plastic deformation [–] Initial kinematic hardening modulus [MPa] Rate of change of elastic range size with plastic deformation [–] Maximum change in size of elastic range [MPa]

Vladimirov et al. (2008) proposed a new type of the wellknown exponential map algorithm which automatically preserves the plastic volume and the symmetry of the internal variables. This integration scheme relies on the use of the spectral decomposition to evaluate the exponential tensor functions in closed form. In the latter study we compared the exponential map with spectral decomposition to two modified backward Euler schemes that preserve plastic incompressibility. The results showed higher accuracy and improved performance when using large time steps in favour of the exponential algorithm. Consider first the plastic flow rule, written here in the more general format ˙ ˙ C˙ p = f(C, Cp , Cpi ) = g(C, Cp , Cpi )Cp

(6)

where the short hand notation f is defined by YD Cp

f = 2

D T

= fC−1 p Cp = g Cp

(7)

YD · (Y )

Since C˙ p is a symmetric tensor, it follows that also f is symmetric. Due to the fact that the tensor g is not symmetric, direct use of the exponential map for the design of an implicit integration formula ¯ :=  t) ˙ ( ¯ Cp = exp(g)C pn

(8)

would require the use of the series representation exp A =  ∞ 1 k

Ykin = Cp X

YD Cp C˙ p = 2˙ 

 y b c ˇ Q

 ˙ =

2 ˙  3

A . However, it can be shown that the series representation k=0 k! ¯ can be replaced by the expression Up exp(U ¯ −1 fU−1 )U−1 of exp(g) p p p √ where Up = Cp denotes the plastic right stretch tensor. The final form of the integrated plastic flow rule then reads −1 −1 ¯ −1 −1 ¯ −1 −1 −1 Cp C−1 pn Cp = Up exp(Up fUp )Up ⇒ Cpn = Up exp(Up fUp )Up

(9)

• Yield function



 =

D

D T

Y · (Y ) −

2 (y − R), 3

−ˇ

R = −Q (1 − e

)

• Kuhn–Tucker conditions ˙ ≥ 0,

 ≤ 0,

˙  =0

Here, X is the back stress tensor in the reference configuration, Y and Ykin are non-symmetric stress-like quantities. The internal variables of the model are Cp = FTp Fp , Cpi = FTpi Fpi and  which describe the evolution of plastic deformation, the evolution of kinematic hardening, and the evolution of isotropic hardening, respectively. There is a total number of seven material parameters (see Table 1). 3. Exponential map algorithm In this section, a short summary of the algorithm used for the numerical integration of the evolution equations is presented.

An important advantage of this format lies in the fact that using −1 the spectral decomposition for U−1 p fUp , Eq. (9) can be computed in closed form. The integration of the evolution equation for kinematic hardening is performed analogously. Without derivation the integration formula for the kinematic hardening rule is given by −1 −1 −1 ¯ −1 C−1 pi n = Upi exp(Upi fkin Upi )Upi

(10)

The final form of the two discretized evolution equations, together with the yield criterion, is summarized below in a residuum format: r1 r2 r3

−1 ¯ −1 −1 −1 = −C−1 pn + Up exp(Up fUp )Up = 0 −1 ¯ −1 fkin U−1 )U−1 = 0 = −C−1 + U exp( U pi n pi pi pi pi ==0

(11)

Due to the symmetry of the internal variables, the symmetric tensors r1 and r2 can be represented in Voigt form as sixdimensional arrays. Thus, the resulting system of equations consists of 13 nonlinear scalar equations which have to be solved iteratively by means of e.g. the Newton method. This local iteration is performed in every time step until the root-mean-square value of the

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

4065

total residuum becomes smaller than a prescribed tolerance. The −1 variables that have to be updated during the iteration are U−1 p , Upi ¯ and . 4. Results and discussion In this section the ability of the model to predict the springback of sheet metal parts is investigated. First, the draw-bend test of DP600 dual phase steel sheets is simulated and the results of various simulations are compared to data from experiments. In addition, the model is applied to the simulation of the S-rail benchmark problem and the thermoforming of thermoplastic blends. 4.1. Draw-bending of a sheet metal strip The first example deals with the process of draw-bending (also termed “bending under tension”) of a narrow sheet metal strip. The draw-bend test was developed by Matlock and coworkers for the investigation of springback of sheet metal parts (Valance and Matlock, 1992). 4.1.1. Process description The draw-bend test is useful for springback investigation, since it represents a typical example of a sheet metal forming operation which produces springback similarly to industrial forming operations, but has the advantage of simplicity. During this process, the sheet metal strip is subjected to tension, bending and unbending. This represents a kind of cyclic loading and, as already stated in the literature (e.g. Carden et al., 2002), springback is sensitive to the Bauschinger effect, particularly for bend/unbend operations (i.e. cyclic loading). The latter statement clarifies the importance of kinematic hardening in the context of springback modelling. As pointed out by e.g. Zhao and Lee (2001), isotropic hardening alone is insufficient, and accurate prediction of springback requires the use of a material model which combines both isotropic and kinematic hardening. The experimental setup for the draw-bend test is shown in Fig. 1. The draw-bend test consists of two hydraulically actuated clamp jaws oriented perpendicularly to each other, with a fixed or rotating cylinder at the intersection of their action lines to simulate a tooling radius over which the sheet metal is drawn. The left actuator provides a constant back force while the right actuator moves at a constant speed. The test is conducted in the following manner: the strip of sheet metal is held fixed between the clamp jaws and is being drawn and bent over the cylindrical tool (roller). In this work, the case of a lubricated roller (minimum friction) is considered. After having been drawn up to a prescribed drawing distance x, the strip is unloaded and removed from the machinery to observe its spring-

Fig. 1. Draw-bend test machine and experimental setup, from Kleiner et al. (2005).

Fig. 2. Springback of DP600 sheet metal specimens, from Kleiner et al. (2005).

Fig. 3. Draw-bend simulation: simple bending to obtain the initial L-shape.

back (Fig. 2). Fig. 2 shows two DP600 steel sheet specimens after springback. The change of shape of the specimen is dominated by the bend-unbend area of the strip (so-called ‘sidewall curl’). Sidewall curl is a springback phenomenon resulting in the curvature of the sheet. The shape change associated with simple bending (near the roller) is much smaller. Figs. 3–5 show how the simulation of the draw-bend test is carried out in this study. Fig. 3 depicts the position of the sheet strip, after it has been bent over the cylindrical tool. By applying a tensile force at the right end and prescribing a displacement at the left one, the sheet is being drawn over the tool to the end position, illustrated

Fig. 4. Draw-bend simulation: end position of the sheet drawn over the tool.

4066

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

Fig. 7. Fit of material model for DP600 steel sheets with thickness 1 mm, only isotropic hardening.

4.1.2. Parameter identification Next, we validate the material model discussed in Section 2 by first identifying the material parameters based on experimental stress–strain curves and then performing finite element (FE) simulations of the draw-bend process and comparing the predicted springback with experiments. For this purpose, we use test data for a DP600 dual-phase steel. In Figs. 6–8 a flow curve, obtained by averaging a number of uniaxial tensile tests on DP600 sheet specimens with a thickness of 1 mm (in the rolling direction), is represented by means of dotted symbols. Due to the lack of compression data (no information about the Bauschinger effect) it is difficult to uniquely determine the hardening parameters by means

of only this experiment. Note that the magnitude of the Bauschinger effect is seldom reported in the literature because of the instability of metal sheets subjected to in-plane compression. It should be noted that the yield stress y , the Young’s modulus E and the Poisson’s ratio  (hence also the Lamé constants  and ) are obtained directly from measurements and do not have to be identified. We employ the following parameter identification strategy: first, the model is fitted based on the experimental flow curve by using only kinematic hardening (Fig. 6). This means that the isotropic hardening parameters ˇ and Q have zero values and only the kinematic hardening parameters b and c are varied. Next, the fitting of the model is performed by using only isotropic hardening (Fig. 7) whereas only the isotropic hardening parameters ˇ and Q are varied this time. In this way, four ‘start’ values for the hardening parameters are obtained (see Table 2). Note the differences in the model curves in Figs. 6 and 7. They both fit very well to the experimental data in the first quadrant (tensile stresses, tensile strains) but differ considerably in the compressive range. One clearly notices the Bauschinger effect in Fig. 6(straining in one direction reduces the yield stress in the opposite direction), and the stress saturation controlled by the isotropic hardening parameters in Fig. 7. Finally, both kinematic and isotropic hardening (i.e. combined hardening) are used in the fitting procedure (Fig. 8). Beginning with the start values of the four hardening parameters b, c, ˇ and Q it turns out that approximately equal amounts of kinematic and isotropic hardening are needed to obtain a good cor-

Fig. 6. Fit of material model for DP600 steel sheets with thickness 1 mm, only kinematic hardening.

Fig. 8. Fit of material model for DP600 steel sheets with thickness 1 mm, combined hardening.

Fig. 5. Draw-bend simulation: removal from tooling and springback.

in Fig. 4. This is the best way to obtain an optimal match between the modelled process and the real situation. Then, loads and constraints (except for all three constraints of three nodes on the left straight part) are removed and the strip is allowed to “spring back” (Fig. 5). Figs. 3–5 represent the distribution of the von Mises stress over the volume of the sheet metal strip in the range 0.08 ≤ V ≤ 789.6 MPa. The decrease of the stresses during the springback phase is due to the elastic unloading.

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

4067

Table 2 Material parameters for thickness 1 mm. Fit

 [MPa]

[MPa]

y [MPa]

b [–]

c [MPa]

ˇ [–]

Q [MPa]

Kinematic Isotropic Combined

72,509.4 72,509.4 72,509.4

140,753.5 140,753.5 140,753.5

325 325 325

7.75 0 12.5

2750 0 2600

0 20 8

0 420 215

Table 3 Material parameters for thickness 2 mm. Fit

Combined

 [MPa] [MPa] y [MPa] b [–] c [MPa] ˇ [–] Q [MPa]

74269.2 144169.6 400 13.5 2600 7.5 215

Table 4 Material parameters for thickness 0.6 mm.

Fig. 9. Fit of material model for DP600 steel sheets with thickness 2 mm, combined hardening.

relation between model and experiment. The finally identified set of material parameters is shown in Table 2. The same strategy is followed to identify the material parameters of DP600 sheets with a thickness of 2 mm and 0.6 mm. It should be noted that due to the difference in the sheet thicknesses (i.e. different level of rolling of the sheet specimen), different material properties are observed, mirrored by different flow curves. The two experimental uniaxial tensile stress–strain curves are presented in Figs. 9 and 10 where in both cases the combined hardening model has been used for the fitting. In both, Figs. 9 and 10, the experimental data are depicted by means of dots and the model curves are represented by continuous lines. The obtained parameters are listed in Tables 3 and 4.

Fit

Combined

 [MPa] [MPa] y [MPa] b [–] c [MPa] ˇ [–] Q [MPa]

71353.4 138509.5 350 6.5 1765 2.25 205

4.1.3. Springback prediction In the previous subsection three sets of material parameters have been identified based on experimental stress–strain curves for DP600 steel sheets with thicknesses of 0.6, 1 and 2 mm, respec-

Fig. 11. FE simulation vs. experiment: thickness 1 mm, radius 10 mm.

tively. Here, FE simulations of the draw-bend test are conducted. We consider three different sheet thicknesses and three different roller radii, namely 5, 10 and 15 mm. Thus, radius-thickness ratios R/t (R-roller radius, t-sheet thickness) of 2.5, 5, 7.5, 8.33, 10, 15, 16.67 and 25 are investigated. In other words, we perform altoTable 5 Geometry parameters.

Fig. 10. Fit of material model for DP600 steel sheets with thickness 0.6 mm, combined hardening.

Radius R [mm]

Length L [mm]

Draw distance x [mm]

5 10 15

356 432 509

95 157 195

4068

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

Fig. 12. FE simulation vs. experiment: thickness 1 mm, radius 5 mm.

gether nine draw-bend simulations. For each FE simulation the predicted springback is compared to the measured one. Table 5 summarizes the geometry parameters (sheet length L, drawing distance x) for each of the three radii of the cylindrical tool. The tensile stress applied at the right end of the sheet is equal to half the material yield stress. The width of all sheets is 50 mm. The implicit FE simulations are performed by means of the finite element software package ABAQUS/Standard (ABAQUS, 2003), where the combined hardening material model from Section 2 has been implemented as a user subroutine UMAT. After investigating the mesh sensitivity of the FE result, a converged mesh of 9600 (length × thickness × width = 300 × 8 × 4) C3D8R brick elements with reduced integration and hourglass stabilization has finally been used for all simulations. It should be noted that the influence of the hourglass stiffness parameter in C3D8R on the springback results in this study has been found to be negligible. Consider first the case t = 1 mm, R = 10 mm. Utilizing the material parameters for ‘only kinematic’ hardening, ‘only isotropic’

Fig. 13. FE simulation vs. experiment: thickness 1 mm, radius 15 mm.

Fig. 14. FE simulation vs. experiment: thickness 2 mm, radius 10 mm.

hardening, and combined hardening (Table 2), we conduct three FE draw-bend simulations. The predicted sheet strip profiles after springback are plotted in Fig. 11 together with the experimentally measured strip profile, represented by a sequence of dots. Note that due to the high number of measurement points along the sheet profile the experimental result looks as if represented by a thick line which is not the case. It is obvious that using kinematic hardening alone underpredicts the strip springback whereas by taking into account only isotropic hardening, springback is overpredicted. Note that this result is in agreement with the observation of Li et al. (2002) of pervasive overprediction of springback on the basis of an isotropic hardening assumption. On the other hand, the simulation performed with both kinematic and isotropic hardening predicts the springback very well; the agreement with the experimental data is excellent. It should be emphasized that the FE simulation of the case t = 1 mm, R = 10 mm has been used to find the best combination of isotropic and kine-

Fig. 15. FE simulation vs. experiment: thickness 2 mm, radius 5 mm.

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

4069

Fig. 16. FE simulation vs. experiment: thickness 2 mm, radius 15 mm.

Fig. 18. FE simulation vs. experiment: thickness 0.6 mm, radius 5 mm.

matic hardening. Again, the flow curves (Section 4.1.2) are not appropriate for this purpose because they do not contain data in the compressive stress range. The results of the simulations with the same sheet thickness (t = 1 mm) but with the tool radii R = 5 mm and R = 15 mm are plotted in Figs. 12 and 13, respectively. In both cases, and also from now on, the complete combined hardening model has been applied in the draw-bend simulations. In both, Figs. 12 and 13, the predicted springback (continuous line) correlates very closely with the experimentally measured one (dots). It should be emphasized that no additional fitting has been done here. These are simply the results of another finite element simulation including the previously identified material parameters. Next, the results for the sheets with the thickness t = 2 mm (material parameters from Table 3) are discussed. The case t = 2 mm, R = 10 mm has been used to uniquely identify the hardening parameters, whereas the other two cases t = 2 mm, R = 5 mm and t = 2 mm, R = 15 mm have been simulated without any additional

fitting. In Figs. 14–16 the simulated sheet profiles after springback are presented together with the experimentally measured ones for tool radii R = 10 mm, R = 5 mm and R = 15 mm, respectively. In the first two cases (Figs. 14 and 15) the match between the FE result and the experiment is excellent, only in Fig. 16 there is a slight deviation from the experimental data. The final three cases consider sheets with thicknesses of 0.6 mm and roller radii of 10, 5 and 15 mm, respectively (Figs. 17–19 ). The agreement between simulation and experiment is highly satisfactory for the case t = 0.6 mm, R = 10 mm, which is the case upon which the parameter identification is based. There is a deviation of the FE result from the experimental data for the case t = 0.6 mm, R = 5 mm (Fig. 18), although the prediction of the springback angle (slope of the straight part of the strip) agrees with the experiment. On the other hand, springback is slightly underpredicted by the FE simulation for the case t = 0.6 mm, R = 15 mm (Fig. 19).

Fig. 17. FE simulation vs. experiment: thickness 0.6 mm, radius 10 mm.

Fig. 19. FE simulation vs. experiment: thickness 0.6 mm, radius 15 mm.

4070

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

4.2. S-rail forming The NUMISHEET96 S-rail benchmark problem (Lee et al., 1996) was designed to compare numerical simulations with each other and with experiments for forming, shape distortion, wrinkling, and springback. It represents a kind of a deep drawing operation in which the tools have a curved, so-called ‘S-rail’, shape. The blank is mounted between a die and a blankholder (binder) and is formed into the curved die by means of an upwards moving punch. Fig. 20 shows the geometry and dimensions of the blank and the initial setup (note that the Y -axis is directed along the rolling direction). For a complete description of the process parameters and the geometry of the tools refer to Lee et al. (1996). Figs. 21–23 demonstrate the simulation of the S-rail forming process as performed in this study. In these three figures the distribution of the von Mises equivalent stress over the volume of the blank in the range 0.85 ≤ V ≤ 376.0 MPa is shown. We have simulated the process by means of ABAQUS/Standard and the user

material subroutine UMAT. A fine mesh of 21,648 C3D8R solid elements with reduced integration has been utilized in the simulation. During the S-rail forming process, the punch moves upwards at a constant speed and presses the sheet metal blank into the die opening (Fig. 21). After the punch has reached the prescribed travel distance of 37 mm (Fig. 22), the tools are removed and the blank is allowed to spring back (Fig. 23). This is done in the simulation by removing the contact constraints between the tools and the blank and by fixing all degrees-of-freedom of two nodes on the top surface of the blank. Due to the elastic unloading during the springback phase, a considerable decrease of the stress level is observed (Fig. 23). The identification of the material parameters of the model has been done based on a stress–strain curve for a draw quality mild interstitial free (IF) steel obtained from uniaxial tensile tests (from Lee et al., 1996). We apply again the combined hardening model with the same identification strategy as in Section 4.1(Fig. 24). The identified material parameters read:  =

Fig. 20. S-rail: blank shape and initial setup, from Lee et al. (1996); Y-axis along the rolling direction.

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

4071

Fig. 21. S-rail simulation: forming phase.

79, 231 MPa, = 118, 846 MPa, y = 152 MPa, b = 6, c = 660 MPa, ˇ = 4.5, Q = 145 MPa. It should be noted that a Coulomb’s coefficient of friction fric = 0.11, as prescribed in Lee et al. (1996), has been used in the FE simulations of the S-rail benchmark problem. Following the instructions in Lee et al. (1996), two cases have been simulated: first, with a blankholder force BHF = 200 kN, and then, with a blankholder force set equal to 10 kN. The comparison between the experimental punch force–punch displacement diagram and the one obtained from the FE simulation of the S-rail forming process for the case BHF = 200 kN is shown in Fig. 25. The

simulation result is represented by means of a thick continuous line and the experimental data are depicted by means of line-connected dots. The computed results are in good agreement with the experiment, although the numerical model underestimates the punch force slightly. Fig. 26 shows the punch force–punch displacement diagram for BHF = 10 kN. In that case, the correspondence between simulation and experiment is very good up to a punch displacement of approximately 15 mm. Afterwards, the simulation result appears to underpredict the punch force.

Fig. 22. S-rail simulation: end of forming.

4072

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

Fig. 23. S-rail simulation: springback.

Fig. 24. S-rail: fit of material model for an IF steel.

Fig. 26. S-rail: punch force–punch displacement, blankholder force BHF = 10 kN.

Now, the finite element prediction of springback is investigated. Fig. 27 shows the top surface profile of a cross-section of the S-rail blank, passing through the points I and E (see Fig. 28), before and after springback. The deviations are to be seen mainly in the areas near the flanges, the left flange experiencing the largest springback. Note that the blankholder force used in this simulation is BHF = 10 kN.

Fig. 25. S-rail: punch force–punch displacement, blankholder force BHF = 200 kN.

Fig. 27. S-rail: springback of IE top surface profile, before and after springback.

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

Fig. 28. S-rail: cross-sections through lines IE and JD.

4073

Fig. 30. S-rail: springback of IE top surface profile, blankholder force =200 kN.

4.3. Thermoforming

Fig. 29. S-rail: springback of IE top surface profile, blankholder force =10 kN.

Next, the simulated springback is compared to the springback from experiments. Based upon the simulations with a blankholder force of 10 kN and 200 kN, the springback of the top surface profile of the cross-section through the line IE is plotted versus the experimentally measured springback in Figs. 29 and 30, respectively. The computed result matches the experiment very well for the case BHF = 10 kN (Fig. 29), even though slight deviations can be seen in the areas near the radius of the blankholder. The springback prediction for the simulation with a blankholder force of 200 kN (Fig. 30) is good for the left flange, however the springback of the right flange seems to be underpredicted.

The last example deals with the simulation of the process of thermoforming of TPO (Thermo-Plastic Olefin) sheets. TPO polymer/filler blends are used extensively in the automotive industry. Thermoforming is one of the oldest and most common methods of processing plastic materials. It represents a reshaping of a thermoplastic forming material onto moulded tools at its own forming temperature. Basically, such a process consists of inserting a thermoplastic sheet in a cold state into the forming clamp area, pre-heating it to the desired temperature, and then moving the mould upwards to form the sheet. The trapped air is sucked out with the help of a vacuum system and once cooled a reverse air supply is activated to remove the plastic part from the mould. Figs. 31–33 demonstrate the simulation of the thermoforming process, as performed in this study. A circular thermoplastic sheet (only a quarter discretized in the FE model) is clamped at the nodes along its perimeter and is put into contact with an upwards moving rigid mould (Fig. 31). After the mould has reached the prescribed forming distance (Fig. 32), vacuum is simulated by applying a negative pressure on the top surface of the sheet. This results in the sheet sticking closely to the mould and obtaining its final formed shape (Fig. 33). The FE simulation is again performed by means of ABAQUS/Standard. A fine mesh of approximately 12,000

Fig. 31. Thermoforming simulation: forming phase.

4074

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

Fig. 32. Thermoforming simulation: end of forming.

(radius × circumference × thickness = 80 × 126 × 2) C3D8R reduced integration continuum elements has been utilized. It should be noted that due to the purely mechanical material model formulation, no temperature influences have been taken into account, which is, admittedly, a strong restriction. The motivation for incorporating this application into the present paper comes from the fact that thermoforming involves large elastic strains (up to 15%). The assumption of small elastic strains usually made in metal plasticity is therefore no longer justified. Thus we see here a very interesting example for a situation where large elastic strains

have to be incorporated into the modelling approach. Unlike many other concepts in the literature (e.g. Chaboche, 1989; Yoshida and Uemori, 2003; Svendsen et al., 2006) the present model fulfils this requirement and seems to be well appropriate for the thermoforming simulation of TPO sheets. Indeed, Fig. 34 shows a typical experimental flow curve (dots) of a TPO sheet, where a significant elastic deformation occurs. The model curve (continuous line), used for the identification of the material parameters, fits very well to the experimental data. In addition, based on the FE thermoforming simulation, a tendency of an increase of the sheet

Fig. 33. Thermoforming simulation: after vacuum.

I.N. Vladimirov et al. / Journal of Materials Processing Technology 209 (2009) 4062–4075

Fig. 34. Thermoforming: fit of material model for a TPO sheet.

thickness from the top to the bottom can be found, which is in qualitative agreement with experimental observations. 5. Conclusions In the first part of the paper we have briefly introduced the constitutive equations of the material model used in this study. The constitutive model incorporates arbitrarily large elastic and plastic deformations combining nonlinear Armstrong–Frederick kinematic hardening and nonlinear Voce type isotropic hardening. The main focus of the paper is the study of the applicability of the material model for springback prediction. Due to the known dependence of springback on the Bauschinger effect, we have investigated the sensitivity of springback to using only isotropic hardening, only kinematic hardening, or using both isotropic and kinematic hardening (combined hardening). The results have shown that only by means of the combined hardening model it is possible to accurately predict the sheet springback. In this work we have applied the combined hardening model to the FE simulation of the draw-bend process and of the NUMISHEET ‘96 S-Rail benchmark problem. The numerical examples show that the material model is very suitable for simulating large deformation problems. In particular, the performed forming simulations show the potential of the model to simulate springback. Comparisons with experiments display an excellent agreement between the finite element results and the experiment data. By all means, the results demonstrate that the discussed finite strain combined hardening model is very well suited for springback prediction. Acknowledgements The authors would like to thank Alexander Brosius, Dr.-Ing. and Rajarajan Govindarajan, M.Sc. (IUL, TU Dortmund, Germany) for the development of the experiments of the draw-bend test accompanying this work. We are also grateful to M.Sc. Pengjie Sun for his assistance during the preparation of this study. References Barlat, F., Brem, J.C., Yoon, J.W., Chung, K., Dick, R.E., Choi, S.H., Pourboghrat, F., Chu, E., Lege, D.J., 2003. Plane stress yield function for aluminium metal sheets. International Journal of Plasticity 19, 1297–1319. Carden, W.D., Geng, L.M., Matlock, D.K., Wagoner, R.H., 2002. Measurement of springback. International Journal of Mechanical Sciences 44, 79–101.

4075

Chaboche, J.L., 1989. Constitutive equations for cyclic plasticity and cyclic viscoplasticity. International Journal of Plasticity 5, 247–302. Chan, K.C., Wang, C., 1999. Theoretical analysis of springback in bending of integrated circuit leadframes. Journal of Materials Processing Technology 91, 111– 115. Choi, Y., Han, C.-S., Lee, J.K., Wagoner, R.H., 2006a. Modeling multi-axial deformation of planar anisotropic elasto-plastic materials. Part I. Theory. International Journal of Plasticity 22, 1745–1764. Choi, Y., Han, C.-S., Lee, J.K., Wagoner, R.H., 2006b. Modeling multi-axial deformation of planar anisotropic elasto-plastic materials. Part II. Applications. International Journal of Plasticity 22, 1765–1783. Chung, K., Lee, M.G., Kim, D., Kim, C., Wenner, M.L., Barlat, F., 2005. Spring-back evaluation of automotive sheets based on isotropic-kinematic hardening laws and non-quadratic anisotropic yield functions. Part I. theory and formulation. International Journal of Plasticity 21, 861–882. Gardiner, F.J., 1957. The springback of metals. Transactions of ASME 79, 1–9. Haddadi, H., Bouvier, S., Banu, M., Maier, C., Teodosiu, C., 2006. Towards an accurate description of the anisotropic behaviour of sheet metals under large plastic deformations: modelling, numerical analysis and identification. International Journal of Plasticity 22, 2226–2271. Hibbit, Karlsson & Sorensen, ABAQUS/Standard User’s Manual, Version 6.3 (2003). Jeunechamps, P.-P., Ho, K.C., Lin, J., Ponthot, J.-P., Dean, T.A., 2006. A closed form technique to predict springback in creep age-forming. International Journal of Mechanical Sciences 48, 621–629. Kleiner, M., Schikorra, M., Govindarajan, R., Brosius, A., 2005. Springback analysis of sheet metals regarding material hardening. In: Geiger, M., Duflou, J., Kals, H.J.J., Shirvani, B., Singh, U.P. (Eds.), Proceedings of the 11th International Conference on Sheet Metal, pp. 712–728. Lee, J.K., Kitzel, G.L., Wagoner, R.H. (Eds.), 1996. NUMISHEET. The Ohio State University. Lee, M.G., Kim, D., Kim, C., Wagoner, R.H., Wenner, M.L., Chung, K., 2005a. Springback evaluation of automotive sheets based on isotropic-kinematic hardening laws and non-quadratic anisotropic yield functions. Part II. Characterization of material properties. International Journal of Plasticity 21, 883–914. Lee, M.G., Kim, D., Kim, C., Wenner, M.L., Chung, K., 2005b. Spring-back evaluation of automotive sheets based on isotropic-kinematic hardening laws and nonquadratic anisotropic yield functions. Part III. Applications. International Journal of Plasticity 21, 915–953. Lee, M.G., Kim, D., Wagoner, R.H., Chung, K., 2007. Semi-analytic hybrid method to predict springback in the 2D draw bend test. Journal of Applied Mechanics 74, 1264–1275. Li, K.P., Carden, W.D., Wagoner, R.H., 2002. Simulation of springback. International Journal of Mechanical Sciences 44, 103–122. Lion, A., 2000. Constitutive modelling in finite thermoviscoplasticity: a physical approach based on nonlinear rheological models. International Journal of Plasticity 16, 469–494. Moon, Y.H., Kim, D.W., Van Tyne, C.J., 2008. Analytical model for prediction of sidewall curl during stretch-bend sheet metal forming. International Journal of Mechanical Sciences 50, 666–675. Oliveira, M.C., Alves, J.L., Chaparro, B.M., Menezes, L.F., 2007. Study on the influence of work-hardening modelling in springback prediction. International Journal of Plasticity 23, 516–543. Papeleux, L., Ponthot, J.-P., 2002. Finite element simulation of springback in sheet metal forming. Journal of Materials Processing Technology 125/126, 785–791. Pourboghrat, F., Chu, E., 1995. Springback in plane strain stretch/draw sheet forming. International Journal of Mechanical Sciences 36, 327–341. Pourboghrat, F., Chung, K., Richmond, O., 1998. A hybrid membrane/shell method for rapid estimation of springback in anisotropic sheet metals. Journal of Applied Mechanics 65, 671–684. Queener, C.A., De Angelis, R.J., 1968. Elastic springback and residual stresses in sheet metal formed by bending. Transactions of ASME 61, 757–768. Ragai, I., Lazim, D., Nemes, J.A., 2005. Anisotropy and springback in draw-bending of stainless steel 410: experimental and numerical study. Journal of Materials Processing Technology 166, 116–127. Svendsen, B., Levkovitch, V., Wang, J., Reusch, F., Reese, S., 2006. Application of the concept of evolving structure tensors to the modeling of initial and induced anisotropy at large deformation. Computers and Structures 84, 1077– 1085. Valance, D.W., Matlock, D.K., 1992. Application of the bending-under-tension friction test to coated sheet steels. Journal of Materials Engineering and Performance 1 (5), 689–693. Vladimirov, I.N., Pietryga, M.P., Reese, S., 2008. On the modelling of non-linear kinematic hardening at finite strains with application to springback—comparison of time integration algorithms. International Journal for Numerical Methods in Engineering 75, 1–28. Wang, C., Kinzel, G., Altan, T., 1993. Mathematical modelling of plane strain bending of sheet and plate. Journal of Materials Processing Technology 39, 279–304. Yoshida, F., Uemori, T., 2003. A model of large-strain cyclic plasticity and its application to springback simulation. International Journal of Mechanical Sciences 45, 1687–1702. Zhao, K.M., Lee, J.K., 2001. Material properties of aluminium alloy for accurate drawbend simulation. Journal of Engineering Materials and Technology 123, 287–292.