- Email: [email protected]
shell method for calculating springback of anisotropic sheet metals undergoing axisymmetric loading
International Journal of Plasticity 16 (2000) 677±700 www.elsevier.com/locate/ijplas
A hybrid membrane/shell method for calculating springback of anisotropic sheet metals undergoing axisymmetric loading Farhang Pourboghrata,*, Michael E. Karabinb, Richard C. Beckerb, Kwansoo Chungc a
Michigan State University, Department of Mechanical Engineering, East Lansing, MI 48824, USA b Alcoa Technical Center, Alcoa Center, PA 15069-0001, USA c Seoul National University, Department of Fiber and Polymer Science, Seoul, South Korea Received in final revised form 10 October 1999
Abstract To reduce the computation cost of ®nite element analyses aiding die design for sheet metal stamping, a hybrid membrane/shell method was developed to determine the springback of anisotropic sheet metal undergoing axisymmetric loading. The hybrid membrane/shell method uses a membrane model to analyze the stamping operation. The bending/unbending strains and stresses varying through thickness are calculated analytically from the incremental shape determined by the membrane analysis. These bending strains and stresses and the ®nal membrane shape are used with a shell ®nite element model to unload the sheet and calculate springback. The accuracy of the springback prediction with the hybrid method was veri®ed against the springback of 2036-T4 aluminum and a DQAK steel sheet drawn into a cup. It was found that, in comparison with a full shell model, a minimum of 50% CPU time saving and a comparable accuracy was achieved when the hybrid method was used to predict springback. # 2000 Elsevier Science Ltd. All rights reserved. Keywords: Axisymmetric loading; Anisotropic sheet metals; Hybrid membrane/shell method
1. Introduction The ®nite element analysis of forming processes needs improvement in computation time for practical industrial applications. One popular scheme to reduce computation time in a sheet metal forming analysis is to assume a thin sheet to be a membrane, neglecting the variation of deformation in the thickness direction (Frey * Corresponding author. 0749-6419/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S0749-6419(99)00067-4
678
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
Nomenclature d1±3 D E h k1 ; k2 K L n~ P R t T; M u; w " l1; l2
principal strain-rate components strain rate tensor Young's modulus plastic hardening parameter principal curvatures of the shell strength coecient fourth order elastic tensor strain hardening exponent unit normal to the ¯ow potential surface normal anisotropy parameter thickness of the shell tension and bending moment displacements of the mid-surface of the shell eective strain principal stretches of the shell (updated Lagrangian) Poisson ratio
: ~ ~ ~ y
Jaumann rate of stress tensor material time derivative of stress tensor stress tensor initial yield strength of material eective stress time
r
Subscripts and superscripts total Lagrangian parameter ( )o o () reference state ( )o initial state
( )
unit vector
and Wenner, 1987; Sklad and Siekirk, 1990; Saran and Wagoner, 1991; Wenner, 1992). Numerical calculations based on membrane elements are cost-eective to analyze strain localization, especially if deformation gradient in the thickness direction is not so large. In bending dominant problems (e.g. ¯anging), however, the deformation gradient and, therefore, the stress gradient in the thickness direction become signi®cant (Stoughton, 1985; Choudhry and Lee, 1989). In such cases, membrane solutions are not valid. Furthermore, when external forces and moments are unloaded after bending (and unbending) materials undergo springback, which is a physical consequence of the through-thickness stress gradient developed during bending. The
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
679
springback is an important measure to analyze for process optimization (Kara®llis and Boyce, 1992). Membrane calculations do not properly predict springback because the stress gradient is not included in the calculation. To account for the bending in sheet forming, shell elements have been developed by adding bending capability to a membrane element (Wang and Tang, 1986; Kubli, et al., 1991; Wagoner et al., 1995; Esche et al, 1999). Another alternative to account for the bending is to superpose the bending eect onto the membrane calculations incrementally (Pourboghrat and Chandorkar, 1992; Pourboghrat and Chu, 1995; Pourboghrat et al., 1998, 1999). In this method, the bending eect is added analytically after membrane solutions are obtained separately: a sequential hybrid of numerical and analytical solutions. In the shell calculations, solutions are obtained by solving nonlinear equilibrium equations of forces and moments simultaneously; therefore, the sequential hybrid method is more cost-eective. In this work, a method based on the hybrid scheme is further developed for the axisymmetric loading condition. The hybrid method consists of a membrane analysis for loading, an intermediate analytical calculation of the through-thickness stress gradients due to bending and unbending, followed by a shell analysis for unloading. In this proposed hybrid method, the incremental bending/unbending deformation strains are calculated from the previously calculated incremental shapes from the membrane ®nite element analysis (e.g. ABAQUS-implicit in this paper). The stress is then updated considering the Jaumann rate of stress calculated from an elastoplastic constitutive model based on the superposed total strain increments (stretching + bending + unbending). As for the superposition of the bending/unbending on membrane stretches, monotonic loading is assumed in which true (or logarithmic) strains are proportional incrementally (during a small ®nite interval). However, the whole strain path is non-proportional. For such deformation, stress is obtained from the incremental theory of elasto-plasticity based on minimum plastic work paths during each increment (Hill, 1948; Nadai, 1963; Chung and Richmond, 1993, 1994). To calculate springback, an ABAQUS model is constructed using shell elements where the initial shape and plastic strains are those of the sheet at the end of the loading step and initial stresses are equal to the total through-thickness stresses calculated analytically. After performing an equilibrium check, ABAQUS proceeds to unload the sheet by removing the forming tools. In this paper, Hill's 1948 yield function, accounting for the normal anisotropy, coupled with isotropic hardening is used to calculate stresses in reverse loading. 2. Constitutive equations During loading, Hooke's law is used to calculate stress below the elastic limit, i.e. 5y , where y is the initial yield stress of the sheet obtained from a uniaxial tensile test. Beyond the elastic limit, i.e. > y , the co-rotational time derivative of stress (Jaumann stress rate) is calculated, for a given strain rate, from an elasto-plastic constitutive equation (Becker, 1992):
680
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
2 6 r 6 L ~ 4~ P h~
3 L : PPL 7 ~ ~~~ 7:D 5 ~ : ~ P :L :P ~ ~ ~
1
r Here and D De Dp are the Jaumann rate of stress and strain rate tensors, ~ tensor, is the material ¯ow strength, h @=@" is ~ respectively, ~ is the~ stress ~ the plastic hardening L is the fourth order elastic tensor and
parameter,
q
P @ =@ = @ =@ , where @ =@ @ =@ : @ =@ , is the second order ten~ ~ ~ ~ ~ ~ sor representing the unit normal to the ¯ow potential. The plastic strain rate, associated with Eq. (1), is calculated from the following expression: r
: P :L :PÿP : ~ ~ " ~ ~ ~ P:L:P P: ~ ~ ~ ~ ~
2
ÿ The fourth order elastic tensor L Lijkl used in this work is the standard tensor for ~ only two independent components. the isotropic elasticity, which has The plasticity algorithm is implemented in a semi-implicit manner for computational eciency. The plastic strain rate is adjusted such that Eq. (5) is satis®ed to within a tolerance, but P in Eqs. (1) and (2) is evaluated only at the beginning of the ~ gradient approximation has the same accuracy as the increment. This forward backward gradient (®rst order) but it is not unconditionally stable. Hence, the time step must be restricted. While avoiding the computation of P saves little for quad~ ratic yield functions, the savings could be substantial for complex yield surfaces (e.g Barlat et al., 1997). It is assumed that the material's ¯ow strength, obtained from a uniaxial tensile test, can be characterized in the plastic range by a general power-law equation: K " "o n
3
ÿ where " d" is the eective plastic strain, K is the strength coecient, n is the strain-hardening exponent and "o is a ®tting parameter. Using Eq. (3), the plastic hardening parameter h can be calculated as: h kn " "o nÿ1
4
As for the yield function, Hill's (1948) quadratic model for rigid plastic sheet metals with normal anisotropy is used: 2 ÿ 3 2 3 ÿ 1 2 R 1 ÿ 2 2 ÿ R 1 2 0
5
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
681
where R "width ="thickness is the normal anisotropy parameter obtained from a uniaxial tensile test. Also, in this paper, we have only considered isotropic hardening where it is assumed that the shape of the initial yield surface remains unchanged but the size grows radially as material work hardens. Discussions on the kinematic and combination type isotropic and kinematic hardening models are referred to the previous work (Pourboghrat et al., 1998). In the membrane analysis as well as in the shell model used for unloading, the yield function must also be available in the FEM code. In order to evaluate any nonstandard yield functions, user material (UMAT) subroutines can be written for use with ABAQUS. To represent the initial planar anisotropy of aluminum sheet alloys for future 3D sheet forming applications, non-quadratic (Barlat et al., 1997) yield functions are available. 3. Simpli®cations for axisymmetric loading For a planar isotropic material, when an initially circular blank is loaded axisymmetrically, the in-plane shear strain and stress vanish, i.e. "12 0; 12 0 (where subscripts 1 and 2 correspond to the radial and circumferential directions, respectively). Also, for a very thin sheet, it is reasonable to assume that a state of plane r stress exists in the sheet, i.e. 3 3 0 (where subscript 3 corresponds to the thickness direction). Furthermore, by adopting Kirchho's assumption that plane sections remain plane and normal to the middle surface of the sheet, shear strains in the thickness direction can be neglected; i.e. "13 "23 0. Then, by using the above assumptions, the yield function in Eq. (5) can be simpli®ed to the following expression: 22 22 R 1 ÿ 2 2 ÿ R 1 2 0
6
where 1 and 2 are principal stresses which are the radial and the circumferential inplane stresses, respectively. ÿ The strain rate tensor, D dij , simpli®es to only having three principal components, d1 d11 , d2 d22~,and d3 d33 . But, only two of these components are independent, since the strain rate in the thickness direction, d3, can be found as a r function of d1 and d2. By writing Eq. (1) in the indicial form and then setting 3 0 (plane stress) it is possible to obtain the plane stress moduli and an expression for d3 as a function of d1 and d2 (Osias, 1972). The ®nal form of Eq. (1) after the simpli®cation becomes: r
K d
7
where K is the plane stress elastic±plastic moduli and the Greek symbols vary between 1 and 2.
682
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
4. Kinematics of axisymmetric deformation The deformation of the mid-surface of an element in the plane containing the meridian curve will be considered for the axisymmetric case. It is noted that what follows is developed for the analytical portion of the hybrid membrane/shell method and that the treatment of kinematics in ABAQUS or other FEM code is independent. The deformation of the element in the circumferential direction will be obtained from the theory of shells of revolution (Flugge, 1973) and Figs. 1±3 show the midsurface of an element of the sheet at the initial time o , reference time o and at current time as it bends and stretches. Wang and Tang (1986) brie¯y discussed the kinematics of deformation for the case of axisymmetric, based on total Lagrangian formulation (Fig. 3). Choudhry and Lee (1989, 1994) discussed the kinematics of deformation, based on the updated Lagrangian formulation for the case of plane strain (see Fig. 2). In this paper, derivation of principal curvatures, stretches and strains of a shell element undergoing an axisymmetric deformation will be discussed in detail using both total and updated Lagrangian formulations. 4.1. Principal stretches (total Lagrangian formulation) Fig. 3 shows a ¯at element at the initial time o and at the current time o as it bends and stretches. By referring to Fig. 3, the following relationship between the current and the initial position vectors and the total displacement vector for the mid-surface of the element can be written:
Fig. 1. A schematic of a shell of revolution (Flugge, 1973).
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
683
Fig. 2. Shell mid-surface at reference time o and current time (updated Lagrangian).
Fig. 3. Shell element at original time o and current time (total Lagrangian).
x XU ~ ~ ~
8
where the initial position vector X and the total displacement vector U are, de®ned ~ in the global coordinate system, ~
684
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
X1 ; X 0 ~
u U w ~
Therefore, the position vector x becomes: ~ X1 u x w
9
10
Consider the unit tangent vector, a^, of the current mid-surface to calculate the ~ i.e principal stretch in the radial direction; dx dx dS dx 1 1 u;X1 1 ~ ~ ~ a^ w;X1 ds dS ds dX1 lo1 lo1 ~
11
where lo1 ds=dS. Here, s and S are travel lengths of the current and initial surfaces from center, respectively, and the comma stands for the dierentiation with respect to X1. Then, by calculating the magnitude of vector a^ from Eq. (11) and setting it ~ equal to unity, lo1 is
q
a^ a^ a^ 1:0 12
~ ~ ~ lo1
hÿ
i12 2 1 u;X1 w;2X1
13
The principal stretch in the circumferential direction, lo2 , is calculated from the change in the radial position of the element; i.e. lo2
X1 u u 1 X1 X1
14
4.2. Principal curvatures and stretches (updated Lagrangian) After bending and stretching, the principal mid-surface curvatures, k1 and k2, of a shell element at the current con®guration, o , can be calculated from the known information about the element at the reference con®guration (time o , see Fig. 2), as follows: r R uA^ wN^
15
where R is the reference con®guration at time o and u and w are incremental dis~ placements de®ned in Fig. 2. In Eq. (15), the unit tangent vector A^ and the unit ~ principal normal vector N^ to the mid-surface of the reference are de®ned as ~
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
685
@R A^ ~ R;s @S ~ ~
16
A;s N^ ~ K1 ~
17
where K1 is the curvature. Therefore, N^ :A^ 0 ~ ~
18
ÿK1 A^ N^ ;s ~ ~
19
and
To calculate the current curvature, k1, the unit tangent vector, a^ a= a and ~ ~ ~
the unit principal normal vector of the mid-surface of the shell, n^ n= n , are ~ ~ ~ calculated. Using Eq. (15), the tangent vector a is found as: ~ @r 20 a ~ r;s R;s u;s A^ uA^ ;s w;s N^ wN^ ;s ~ ~ ~ ~ ~ ~ @S ~ By substituting from Eqs. (16), (17) and (19) into (20), and re-arranging, the following expression results: @r ÿ ÿ a ~ r;s 1 u;s ÿ K1 w A^ w;s K1 u N^ cA^ dN^ ~ ~ ~ ~ ~ @S ~
21
The principal incremental stretch of the mid-surface in the radial direction is calculated from the magnitude of the base vector in Eq. (21):
q p 2 2
22 l1
a~ a~ :a~ c d The current length of the mid-surface of the shell in the radial direction, ds, is calculated from the reference length, dS, and l1 as follows: ds l1 dS The unit principal normal vector of the surface of the current shell, n^, is ~ ^ ^ a^ ;s a^ ;s ÿdA cN ~ ~ n^ ~ ~ l1 ~ a~^ ;s a~^ ;s
23
24
686
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
which, from Eqs. (21), (22) and (24), shows that n^:a^ 0. The current principal cur~ ~ vature of the shell, k1, can now be found as: 1 1 k1 ÿa^:n^ ;s ÿr;s :n^ ;s ÿ 2 r;s :n^ ;s ÿ 2 a:n^ ;s ~ ~ ~ ~ l1 ~ ~ l1 ~ ~ where a is given by Eq. (20) and n^ ;S can be derived from Eq. (24): ~ ~ dn^ l1 : ÿd;s A^ c;s N^ ÿ K1 dN^ ÿ K1 cA^ ÿ l1;S : ÿdA^ cN^ ~ ~ ~ ~ ~ n^ ;s ~ dS ~ l2
25
26
1
In Eq. (26), l1;S is assumed to vanish within an element and the above expression simpli®es as: n^ ;S ~
ÿ ÿ d;S K1 c A^ ÿc;s K1 d N^ ~ ~ ÿ l1
27
By substituting from Eqs. (27), (20) and (20) into Eq. (25), the current centerline curvature of the shell, 0, can be found: k1
cd;s ÿ dc;s K1 l21 l31
28
It should be noted that it is possible to recover Eq. (13), for total Lagrangian formulation, from Eq. (22). This is done by setting the reference curvature K1 equal to zero for a ¯at sheet, replacing S with X1 and treating the incremental displacements u and w as total displacements in Eqs. (21) and (22). Similarly, the current curvature of the shell, k1, can be calculated from the initial con®guration by setting K1 equal to zero, replacing l1 with lo1 and S with X1 and treating u and w as total displacements in Eq. (28) ÿ ÿw;x1 u;x1 x1 1 u;x1 w;x1 x1 29 k1 3 lo1 The current principal stretch increment of the mid-surface of the shell in the circumferential direction, l2 , is found from the geometry (in Fig. 2) as: l2
r1 R1
30
The principal curvature of the shell in the circumferential direction, k2, is found from Fig. 1 (Flugge, 1973), as follows:
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
k2
1 sin 1 dw w;s ÿ ÿ ÿ r2 R R ds R
687
31
where w and s were de®ned previously and R is de®ned as (see Fig. 1): R X1 u
32
Since it is known that ds lo1 dX1 , Eq. (31) becomes k2 ÿ
1 dw 1 dw dX1 1 dw w;X ÿ ÿ o ÿ o 1 R ds R dX1 ds l1 R dX1 l1 R
33
Finally, by substituting from Eq. (32) into (33), we will get the following expression for the principal curvature in the circumferential direction: k2 ÿ
w;X1 X1 u:lo1
34
4.3. Strain rate The principal strain rates, d1 and d2 for any through-the-thickness surface are calculated from the principal stretches, l1 and l1, and curvatures, k1 and k2, of the mid-surface of the shell. The principal stretches of any through-the-thickness surfaces l 1; 2, are equal to: l z l 1 zk
35
and through-the-thickness principal strains are equal to: " z ln l z ln l ln 1 zk
36
The strain rates are then calculated by taking the time derivative of Eq. (36) as: h i : : : d ln l z z:k k :z l 37 d l 1 zk dt : where in Eqs. (35)±(37) takes on the values 1 and 2. In Eq. (37), l are the rate of principal stretches l1 and l2 given by Eqs. (22) and (30). The rate of the centerline : curvatures, k , are calculated as the dierence between the current and reference values: : k ÿ ok k
38
688
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
In Eq. (37), z is the current location of a through-the-thickness surface, originally a distance away from the mid-surface. Based on the plastic incompressibility condition and ignoring elastic volume change, the following relationship between z and for each through-the-thickness layer can be obtained: z2 z3 l1 :l2 : z k1 k2 : k1 :k2 : 2 3
39
where ÿto =244to =2 and ÿt=24z4t=2 and parameters to and t are the original and the current thickness of the shell respectively. The current thickness of the shell is updated, after each deformation increment as follows: t ot ed3
40
where o t is the thickness of the shell at the reference state. 5. Stresses In a co-rotational coordinate system, the spins relative to the co-rotational frame are zero and the material derivative of stress is equal to: D r : ~ D ~ ~
41
r
where is given by Eq. (7). The updated current total stress is obtained by adding : ~ the incremental stress, : ,1 to the stress tensor at the reference state, o as ~ ~ ~ follows: o ~ ~ ~
42
where is an increment of time and the reference state is at the beginning of the time increment. 6. The analytical hybrid membrane/shell method To calculate the bending deformation, the hybrid method utilizes the intermediate sheet shape information (x and y coordinates of each node in the sheet for every : : 1 Here we assume 1 =2 1 = 2 because the principal directions are materially ®xed for axisymmetric deformations and the incremental deformations and the incremental deformation theory is supplied.
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
689
increment) generated by the membrane code in a post processing fashion (Pourboghrat and Chu, 1995; Pourboghrat et al., 1998, 1999). The intermediate membrane shapes are used to calculate the bending strain increment for each integration point by an analytical model, as explained in Section 4. The bending stress increments are calculated for each integration point using the elasto-plastic constitutive equation described in Section 3. The total (bending + stretching) deformation at the end of each increment is obtained by integrating the stresses and strains over the increments, as described in sections 4 and 5. The pair of tension components calculated by integrating through-the-thickness hybrid shell stresses are T
t=2 ÿt=2
z:dz 1; 2
43
which are the same as the in-plane membrane forces. Eq. (43) decides the location of the neutral surface. The pair of bending moments that are generated by the hybrid shell stresses are M
=2 ÿ =2
z:z:dz
44
which are assumed to be balanced by additional external contact forces between the sheet and the tooling prior to springback calculation (the contact forces in the membrane analysis balance only tension-curvature forces). Springback is calculated with a shell model (e.g. ABAQUS, 1998), using the membrane shape, the hybrid shell stresses and the forming tools in their ®nal position. Prior to the unloading, the equilibrium of the hybrid shell stresses with the external contact forces are enforced (by ABAQUS), with minimal adjustments made by ABAQUS to the membrane solution. 7. Springback results and discussions In the punch-stretching test, a ¯at sheet is initially clamped between a ¯at blank holder (binder) and the bottom die with a speci®ed blank holding force, as shown in Fig. 4. Then, the hemispherical punch bends and stretches the sheet by pushing the sheet into the die cavity. In the punch-stretching test, the blank holding force is very large. There is very little draw-in and the sheet mostly stretches. Fig. 4 shows the simulated middle surface for the axisymmetric punch-stretching test using membrane element in the ®nite element code ABAQUS, (1998). Due to the symmetry, only half of the sheet was used for the ®nite element simulation. Table 1 shows the friction coecient and the material properties of the aluminum alloy and the DQAK steel sheet used in the ®nite element simulations. Table 2 shows the
690
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
tooling speci®cations, the sheet dimensions and the blank holding force used in the punch-stretching test (Stevenson, 1993). In Fig. 4, the deformed shapes of the sheet correspond to the punch displacements of 7.62, 12.70 and 20.3 mm, respectively. Fig. 5a and b, for the aluminum alloy sheet, shows the dierence between the
Table 1 Material properties of the sheet alloys and friction coecient used in modeling Material type
K (Mpa)
n
E (Gpa)
R
Y (Mpa)
"o
friction coecient
AL 2036-T4 DQAK Steel
604.0 550.0
0.214 0.248
69.0 210.0
0.84 2.10
196.0 188.0
0.33 0.30
0.0052 0.0132
0.1 0.1
Table 2 Tooling geometry and boundary conditions Case Sheet thickness Punch radius, Die radius, Die cavity Punch displacement Blank holding (mm) Rp (mm) Rd (mm) opening (mm) (mm) force (KN) A B C
0.68 0.68 0.68
101.6 101.6 101.6
15.0 15.0 15.0
105.0 105.0 105.0
7.6 12.7 20.3
150.0 150.0 150.0
Fig. 4. The simulated middle surface for the axisymmetric punch-stretching test using membrane element (ABAQUS). The deformed shapes correspond to punch depths of 7.62, 12.7 and 20.32 mm, respectively.
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
691
Fig. 5. Simulated axisymmetric cup stretch/drawing test and deformed shapes of the sheet at 7.62 and 20.32 mm punch depths, respectively (with ABAQUS code).
simulated sheet shapes by membrane and shell elements at punch displacements of 7.62 and 20.3 mm, respectively. Fig. 5a shows a small dierence in the two shapes around the die corner, while Fig. 5b shows no dierence in the two shapes due to increased tension in the sheet. The punch-stretching experiments performed by Stevenson (1993) to measure springback were done in two steps. In this paper, similar to the experiments, springback of the sheet was also calculated in two steps using both shell (with one integration point on the mid-surface and ®ve through the thickness) and the hybrid membrane/shell methods (with 11 through-the-thickness layers). In the ®rst step of the springback simulation, only the punch was removed and the blank holder was left engaged (used as a reference surface). The vertical displacement of the sheet at the punch pole was then measured. Figs. 6a±c and 7a±c show the comparison
692
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
Fig. 6. Loaded and unloaded sheet shapes (SPB1):hybrid vs. shell vs. experiment at punch displacements of (a) 7.6 mm (b) 12.7 mm and (c) 20.3 mm for aluminum alloy sheet.
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
693
Fig. 7. Loaded and unloaded sheet shapes (SPB1):hybrid vs. shell vs. experiment at punch displacements of (a) 7.6 mm (b) 12.7 mm and (c) 20.3 mm for DQAK steel sheet.
694
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
between the measured springback by Stevenson (1993), and numerically predicted springback (SPB1) with shell and hybrid methods for aluminum and steel alloy sheets, respectively. It is interesting to note, from Fig. 7c, that the steel sheet undergoes very little springback after being deformed by the punch to a depth of 20.3 mm. The reason for this is that the sheet is highly stretched and the residual stresses in the sheet are negligible. On the other hand, for the same deformation, the aluminum sheet springback by about 1mm (Fig. 6c). Both shell and hybrid models predict similar springback results. In the second step of the springback simulation, both the blank holder and the punch were removed and by holding the last node under the blank holder ®xed (used as a reference point), the sheet was allowed to come to a complete unloading. Figs. 8a±c and 9a±c show the comparison between the shell and the hybrid method's unloaded shape predictions (SPB2). The exact comparison of predicted unloaded sheet shape with the experimental results, for the second step (SPB2), was dicult since unlike the ®rst step (SPB1) there was no known ®xed reference surface for the experimental measurements. In the experiments (Stevenson, 1993), the maximum height in the deformed sheet was measured using the coordinate measuring machine (CMM). To compare predicted and experimental values, the maximum heights in the deformed sheet for the shell and the hybrid models were calculated from Figs. 8a±c and 9a±c, as Hmax= YmaxÿYmin, and then compared against the experimental data. The last three columns in Tables 3 and 4 show the comparison between the predicted and measured maximum height data for aluminum and DQAK steel sheets. 7.1. Shell versus hybrid Figs 6a±9a consistently show that: (1) the aluminum sheet has a larger springback than the steel sheet, (2) the unloaded shapes of the sheet predicted by the hybrid membrane/shell method, for the punch travel of 7.62 mm, are not in good agreement with those predicted by the conventional shell model. The observed dierences in the shapes are caused by a larger springback predicted around the die corner (at X67.5 mm, Figs. 8a and 9a) by the hybrid method. These discrepancies in predicted unloaded shapes, between the hybrid and the shell models, are substantially improved for the punch travel of 12.7 mm (see Figs. 8b and 9b) and are completely eliminated for the punch travel of 20.3mm (see Figs. 8c and 9c). By referring to Fig. 10a, it can be seen that around the die radius area (nodes 160± 170), at the punch displacement of 7.6 mm, the thickness strains ("t 0:1%) predicted by the hybrid method are smaller than the initial yield strain of the sheet ("y 0:28%). In this situation the hybrid method, unlike the shell model, predicts a larger springback due to elastic thickness strains at the die radius area (see Figs. 8a and 9a). But, as the magnitude of the thickness strain and geometric stiening increases, with the deeper punch penetrations, the discrepancy in the unloaded sheet shape between the hybrid and shell models disappears. Fig. 10c shows that for the aluminum sheet, at the punch displacement of 20.3 mm, the thickness strains are plastic ("t ÿ0:5%) around the die radius area and therefore both hybrid and shell models predict similar springback (see Figs. 8c and 9c).
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
695
Fig. 8. Loaded and unloaded sheet shapes (SPB2):hybrid vs. shell at punch displacements of (a) 7.6 mm (b) 12.7 mm and (c) 20.3 mm for aluminum alloy sheet.
696
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
Fig. 9. Loaded and unloaded sheet shapes (SPB2) hybrid vs. shell at punch displacements of (a) 7.6 mm (b) 12.7 mm and (c) 20.3 mm for DQAK steel sheet.
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
697
Table 3 Springback of aluminum 2036-T4 alloy sheet Max. punch height (mm)
SPB1 hybrid (mm)
SPB1 shell (mm)
SPB1 experiment (mm)
SPB2 hybrid Hmax (mm)
SPB2 shell Hmax (mm)
SPB2 experiment (mm)
ÿ7.62 ÿ12.70 ÿ20.30
ÿ5.69 ÿ11.37 ÿ19.39
ÿ5.51 ÿ11.30 ÿ19.36
ÿ5.26 ÿ10.87 ÿ18.82
ÿ4.04ÿ0.59=ÿ4.63 ÿ9.50ÿ0.59=ÿ10.09 ÿ18.26ÿ0.32=ÿ18.58
ÿ4.72ÿ0.01=ÿ4.73 ÿ9.66ÿ0.43=ÿ10.09 ÿ18.26ÿ0.29=ÿ18.55
ÿ5.23 ÿ11.20 ÿ19.28
Table 4 Springback of draw quality aluminum killed (DQAK) steel sheet Max. punch height (mm)
SPB1 hybrid (mm)
SPB1 shell (mm)
SPB1 experiment (mm)
SPB2 hybrid Hmax (mm)
SPB2 shell Hmax (mm)
SPB2 experiment (mm)
ÿ7.62 ÿ12.70 ÿ20.30
ÿ6.95 ÿ12.22 ÿ19.96
ÿ6.83 ÿ12.19 ÿ19.96
ÿ6.95 ÿ12.25 ÿ20.00
ÿ6.03ÿ0.45=ÿ6.48 ÿ11.56ÿ0.27 =ÿ11.83 ÿ19.58ÿ0.14=ÿ19.72
ÿ6.27ÿ0.20=ÿ6.47 ÿ11.58ÿ0.25=ÿ11.83 ÿ19.57ÿ0.15=ÿ19.72
ÿ6.39 ÿ11.70 ÿ19.45
Fig. 10b and d show that when the hybrid method uses the shell solution, the error in predicted thickness around the die corner diminishes. Therefore, the main source of error in springback prediction by the hybrid method comes from the use of the membrane solution. However, the magnitude of these errors is small (Fig. 8a and 9a) and it diminishes as the sheet further deforms and stretches (Figs. 8b and 9b and c). 7.2. Comparison of CPU time between shell and hybrid methods Table 5 shows a comparison between the total CPU seconds taken by the hybrid membrane/shell method and a full shell model to analyze the punch stretching and a two-step springback calculation problem. The hybrid method took about 30 CPU s, on a HP-755 workstation, to compute bending stresses for 100 increments of membrane solutions. Each membrane solution (i.e. sheet shape) was comprised of 200 membrane elements. The shell models, with one integration point on the mid-surface and both 5 and 11 integration points through the thickness, required a much larger CPU time to analyze the same problem, as shown in Table 5. The amount of the CPU time required for calculating springback was almost the same for both hybrid and shell models. This is because the hybrid method also used a shell model to perform the unloading steps. It can be seen that the hybrid membrane/shell model has at least a 50% advantage in CPU time over the conventional shell analyses for this simple problem.
698
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
Fig. 10. (a) Predicted sheet thickness:hybrid using membrane shapes vs. shell at punch displacement of 7.6 mm for aluminum alloy sheet; (b) predicted sheet thickness:hybrid using shell shapes vs. shell at punch displacement of 7.6 mm for aluminum alloy sheet; (c). Predicted sheet thickness:hybrid using membrane shapes vs. shell at punch displacement of 20.3 mm for aluminum alloy sheet; (d) Predicted sheet thickness:hybrid using shell shapes vs. shell at punch displacement of 20.3 mm for aluminum alloy sheet.
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
699
Table 5 Comparison of the CPU time for dierent methods Method
CPU (s)
Number of iterations
Hybrid membrane/shell Full shell, ®ve integration points Full shell, 11 integration points
821 1706 2173
171 191 230
8. Conclusions To reduce the cost of the ®nite element analyses for sheet forming, a hybrid membrane/shell method was developed to study the springback of anisotropic sheet metals undergoing axisymmetric loading. The goal was to apply the hybrid method incrementally in order to determine the bending strain prediction accurately. In this model, the incremental bending/unbending strains and stresses were calculated only from the information obtained from ABAQUS about the shape of the membrane sheet. An ABAQUS model with shell elements was used to unload the ®nal shape of the sheet with initial stresses and strains that were calculated analytically. For veri®cation, the hybrid method was applied to predict springback in an axisymmetric punch-stretching test for 2036-T4 aluminum and DQAK steel materials. The hybrid membrane/shell springback predictions were in good agreement with the ABAQUS full shell analysis. It was observed that there are some dierences between the experimental springback data and those predicted by the numerical models. These discrepancies were attributed to neglecting to model planar anisotropy, particularly for aluminum 2036-T4 alloy, and the possible Bauschinger eect. It is intended to apply the hybrid membrane/shell method to three-dimensional sheet forming. In that case, the eect of modeling planar anisotropy and the Bauschinger eect on predicted springback results for punch stretching will be assessed once again. A minimum of 50% CPU time saving was achieved by using the hybrid method to predict the springback of an axisymmetric cup instead of a full shell model. It is expected to save even more CPU time for the three dimensional problems. Acknowledgements The authors wish to sincerely thank Dr. Owen Richmond of Alcoa for the support they have received during the course of this work. Farhang Pourboghrat would like also to thank Dr. Robin Stevenson of the Engineering Mechanics Department at GM for his helpful discussions and providing experimental samples and data for this work. References ABAQUS, 1998. User's and Theory Manuals, version 5.8. Hibbitt Karlsson Sorensen Inc, RI. Barlat, F., Maeda, Y., Chung, K., Yanagawa, M., Brem, J.C., Hayashida, Y., Lege, D.J., Matsui, K., Murtha, S.J., Hattori, S., Becker, R.C., Makosey, S., 1997. Yield function development for aluminum alloy sheets. J. Mech. Phys. Solids 45, 1727.
700
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
Becker, R.C., 1992. An ecient semi-implicit integration algorithm for anisotropic ¯ow potentials, unpublished work. Choudhry, S., Lee, J.K., 1989, Numerical simulation of thin sheet metal forming processes including bending eects, In: Durham, D., Saigal, A. (Eds.), Proc. microstructural development and control in materials processing (ASME/WAM, ASME MD-Vol. 14), San Francisco, CA, pp. 11±18. Choudhry, S., Lee, J.K., 1994. Dynamic plane-strain ®nite element simulation of industrial sheet-metal forming processes, Int. J. Mech. Sci. 36 (3), 189±207. Chung, K., Richmond, O., 1993. A deformation theory of plasticity based on minimum work path. Int. J. Plasticity. Chung, K., Richmond, O., 1994. The mechanics of ideal forming. J. of Applied Mechanics 61, 176±181. Esche, S.K., Kinzel, G.K., Lee, J.K., 1999. An axisymmetric membrane element with bending stiness for static implicit sheet metal forming simulation. ASME Journal of Applied Mechanics 66, 153±164. Flugge, W., 1973. Stresses in Shells, 2nd Edition. Springer-Verlag, New York. Frey, W.H., Wenner, M.L., 1987. Development and application of a one-dimensional ®nite element code for sheet metal forming analysis. In: Samanta, S.K. et al. (Eds.), Interdisciplinary Issues in Materials Processing and Manufacturing, Vol. 1. Hill, R., 1948. Theory of the yielding and plastic ¯ow of anisotropic metals. Royal Soc. London Proc. 193A, 281. Kara®llis, A.P., Boyce, M.C., 1992. Tooling design in sheet metal forming using springback calculations. Int. J. Mech. Sci. 34 (2), 113±131. Kubli, W., Anderheggen, E., Reissner, J., 1991. Nonlinear solver with uncoupled bending and stretching deformation for simulating thin sheet metal forming. Proc. FE-Simulation of 3-D Sheet Metal Forming Processes In Automotive Industry, 14±16 May, Zurich, Switzerland. Nadai, A., 1963. Theory of ¯ow and fracture of solids, Vol. 2. McGraw-Hill, New York, p.74. Osias, J., 1972. Finite deformation of elasto-plastic solids. Ph.D. thesis, Carnegie Mellon University. Pourboghrat, F., Chandorkar, K., 1992. Springback calculation for plane strain sheet forming using ®nite element membrane solution. Symposium on Numerical Methods for Simulation of Industrial Metal Forming Processes, ASME WAM'92, 8±13 November, Anaheim. Pourboghrat, F., Chu, E., 1995. Springback in plane strain stretch/draw sheet forming. Int. J. Mech. Sci. 36 (3), 327±341. Pourboghrat, F., Chung, K., Richmond, O., 1998. A hybrid membrane/shell method for rapid estimation of springback in anisotropic sheet metals. ASME Journal of Applied Mechanics 65 (3), 671±684. Pourboghrat, F., Yoon, J.W., Chung, K., Barlat, F., 1999. A 3D hybrid membrane/shell method to predict the springback of anisotropic sheet metals. The 7th International Symposium on Plasticity and Its Current Applications, 5±13 January 1999, Cancun, Mexico. Saran, M.J., Wagoner, R.H., 1991. A user's manual and veri®cation examples for section form (Sheet-S/ version 3), section analysis of sheet stampings. Report no. ERC/NSM-S-91-17. Sklad, M.P., Siekirk, J.F., 1990. Proc. 16th IDDRG Congress, ASM International, Bolange, Sweden, pp. 295±303. Stevenson, R., 1993. Springback in simple axisymmetric stampings. Metallurgical Transactions 24A, 925±933. Stoughton, T.B., 1985. Finite element modeling of 1008 ak sheet steel stretched over a rectangular punch with bending eects. In: Wang, N.M., Tang, S.C. (Eds.), Computer Modeling Of Sheet Metal Forming Processes. The Metallurgical Society, Warrendale, PA, pp. 143±159. Wagoner, R.H., Zhou, D., Sriram, S., 1995. The element bending group method. Proc. International Conference on Computational Engineering Science, 30 July±3 August, HI, USA. Wang, N.M., Tang, S.C., 1986. Analysis of bending eects in sheet forming operations. Int. J. Numer. Methods Eng. 25, 253±267. Wenner, M.L., 1992. Elementary solutions and process sensitivities for plane-strain sheet-metal forming. J. of Applied Mechanics 59, 23±28.
A hybrid membrane/shell method for calculating springback of anisotropic sheet metals undergoing axisymmetric loading Farhang Pourboghrata,*, Michael E. Karabinb, Richard C. Beckerb, Kwansoo Chungc a
Michigan State University, Department of Mechanical Engineering, East Lansing, MI 48824, USA b Alcoa Technical Center, Alcoa Center, PA 15069-0001, USA c Seoul National University, Department of Fiber and Polymer Science, Seoul, South Korea Received in final revised form 10 October 1999
Abstract To reduce the computation cost of ®nite element analyses aiding die design for sheet metal stamping, a hybrid membrane/shell method was developed to determine the springback of anisotropic sheet metal undergoing axisymmetric loading. The hybrid membrane/shell method uses a membrane model to analyze the stamping operation. The bending/unbending strains and stresses varying through thickness are calculated analytically from the incremental shape determined by the membrane analysis. These bending strains and stresses and the ®nal membrane shape are used with a shell ®nite element model to unload the sheet and calculate springback. The accuracy of the springback prediction with the hybrid method was veri®ed against the springback of 2036-T4 aluminum and a DQAK steel sheet drawn into a cup. It was found that, in comparison with a full shell model, a minimum of 50% CPU time saving and a comparable accuracy was achieved when the hybrid method was used to predict springback. # 2000 Elsevier Science Ltd. All rights reserved. Keywords: Axisymmetric loading; Anisotropic sheet metals; Hybrid membrane/shell method
1. Introduction The ®nite element analysis of forming processes needs improvement in computation time for practical industrial applications. One popular scheme to reduce computation time in a sheet metal forming analysis is to assume a thin sheet to be a membrane, neglecting the variation of deformation in the thickness direction (Frey * Corresponding author. 0749-6419/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S0749-6419(99)00067-4
678
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
Nomenclature d1±3 D E h k1 ; k2 K L n~ P R t T; M u; w " l1; l2
principal strain-rate components strain rate tensor Young's modulus plastic hardening parameter principal curvatures of the shell strength coecient fourth order elastic tensor strain hardening exponent unit normal to the ¯ow potential surface normal anisotropy parameter thickness of the shell tension and bending moment displacements of the mid-surface of the shell eective strain principal stretches of the shell (updated Lagrangian) Poisson ratio
: ~ ~ ~ y
Jaumann rate of stress tensor material time derivative of stress tensor stress tensor initial yield strength of material eective stress time
r
Subscripts and superscripts total Lagrangian parameter ( )o o () reference state ( )o initial state
( )
unit vector
and Wenner, 1987; Sklad and Siekirk, 1990; Saran and Wagoner, 1991; Wenner, 1992). Numerical calculations based on membrane elements are cost-eective to analyze strain localization, especially if deformation gradient in the thickness direction is not so large. In bending dominant problems (e.g. ¯anging), however, the deformation gradient and, therefore, the stress gradient in the thickness direction become signi®cant (Stoughton, 1985; Choudhry and Lee, 1989). In such cases, membrane solutions are not valid. Furthermore, when external forces and moments are unloaded after bending (and unbending) materials undergo springback, which is a physical consequence of the through-thickness stress gradient developed during bending. The
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
679
springback is an important measure to analyze for process optimization (Kara®llis and Boyce, 1992). Membrane calculations do not properly predict springback because the stress gradient is not included in the calculation. To account for the bending in sheet forming, shell elements have been developed by adding bending capability to a membrane element (Wang and Tang, 1986; Kubli, et al., 1991; Wagoner et al., 1995; Esche et al, 1999). Another alternative to account for the bending is to superpose the bending eect onto the membrane calculations incrementally (Pourboghrat and Chandorkar, 1992; Pourboghrat and Chu, 1995; Pourboghrat et al., 1998, 1999). In this method, the bending eect is added analytically after membrane solutions are obtained separately: a sequential hybrid of numerical and analytical solutions. In the shell calculations, solutions are obtained by solving nonlinear equilibrium equations of forces and moments simultaneously; therefore, the sequential hybrid method is more cost-eective. In this work, a method based on the hybrid scheme is further developed for the axisymmetric loading condition. The hybrid method consists of a membrane analysis for loading, an intermediate analytical calculation of the through-thickness stress gradients due to bending and unbending, followed by a shell analysis for unloading. In this proposed hybrid method, the incremental bending/unbending deformation strains are calculated from the previously calculated incremental shapes from the membrane ®nite element analysis (e.g. ABAQUS-implicit in this paper). The stress is then updated considering the Jaumann rate of stress calculated from an elastoplastic constitutive model based on the superposed total strain increments (stretching + bending + unbending). As for the superposition of the bending/unbending on membrane stretches, monotonic loading is assumed in which true (or logarithmic) strains are proportional incrementally (during a small ®nite interval). However, the whole strain path is non-proportional. For such deformation, stress is obtained from the incremental theory of elasto-plasticity based on minimum plastic work paths during each increment (Hill, 1948; Nadai, 1963; Chung and Richmond, 1993, 1994). To calculate springback, an ABAQUS model is constructed using shell elements where the initial shape and plastic strains are those of the sheet at the end of the loading step and initial stresses are equal to the total through-thickness stresses calculated analytically. After performing an equilibrium check, ABAQUS proceeds to unload the sheet by removing the forming tools. In this paper, Hill's 1948 yield function, accounting for the normal anisotropy, coupled with isotropic hardening is used to calculate stresses in reverse loading. 2. Constitutive equations During loading, Hooke's law is used to calculate stress below the elastic limit, i.e. 5y , where y is the initial yield stress of the sheet obtained from a uniaxial tensile test. Beyond the elastic limit, i.e. > y , the co-rotational time derivative of stress (Jaumann stress rate) is calculated, for a given strain rate, from an elasto-plastic constitutive equation (Becker, 1992):
680
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
2 6 r 6 L ~ 4~ P h~
3 L : PPL 7 ~ ~~~ 7:D 5 ~ : ~ P :L :P ~ ~ ~
1
r Here and D De Dp are the Jaumann rate of stress and strain rate tensors, ~ tensor, is the material ¯ow strength, h @=@" is ~ respectively, ~ is the~ stress ~ the plastic hardening L is the fourth order elastic tensor and
parameter,
q
P @ =@ = @ =@ , where @ =@ @ =@ : @ =@ , is the second order ten~ ~ ~ ~ ~ ~ sor representing the unit normal to the ¯ow potential. The plastic strain rate, associated with Eq. (1), is calculated from the following expression: r
: P :L :PÿP : ~ ~ " ~ ~ ~ P:L:P P: ~ ~ ~ ~ ~
2
ÿ The fourth order elastic tensor L Lijkl used in this work is the standard tensor for ~ only two independent components. the isotropic elasticity, which has The plasticity algorithm is implemented in a semi-implicit manner for computational eciency. The plastic strain rate is adjusted such that Eq. (5) is satis®ed to within a tolerance, but P in Eqs. (1) and (2) is evaluated only at the beginning of the ~ gradient approximation has the same accuracy as the increment. This forward backward gradient (®rst order) but it is not unconditionally stable. Hence, the time step must be restricted. While avoiding the computation of P saves little for quad~ ratic yield functions, the savings could be substantial for complex yield surfaces (e.g Barlat et al., 1997). It is assumed that the material's ¯ow strength, obtained from a uniaxial tensile test, can be characterized in the plastic range by a general power-law equation: K " "o n
3
ÿ where " d" is the eective plastic strain, K is the strength coecient, n is the strain-hardening exponent and "o is a ®tting parameter. Using Eq. (3), the plastic hardening parameter h can be calculated as: h kn " "o nÿ1
4
As for the yield function, Hill's (1948) quadratic model for rigid plastic sheet metals with normal anisotropy is used: 2 ÿ 3 2 3 ÿ 1 2 R 1 ÿ 2 2 ÿ R 1 2 0
5
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
681
where R "width ="thickness is the normal anisotropy parameter obtained from a uniaxial tensile test. Also, in this paper, we have only considered isotropic hardening where it is assumed that the shape of the initial yield surface remains unchanged but the size grows radially as material work hardens. Discussions on the kinematic and combination type isotropic and kinematic hardening models are referred to the previous work (Pourboghrat et al., 1998). In the membrane analysis as well as in the shell model used for unloading, the yield function must also be available in the FEM code. In order to evaluate any nonstandard yield functions, user material (UMAT) subroutines can be written for use with ABAQUS. To represent the initial planar anisotropy of aluminum sheet alloys for future 3D sheet forming applications, non-quadratic (Barlat et al., 1997) yield functions are available. 3. Simpli®cations for axisymmetric loading For a planar isotropic material, when an initially circular blank is loaded axisymmetrically, the in-plane shear strain and stress vanish, i.e. "12 0; 12 0 (where subscripts 1 and 2 correspond to the radial and circumferential directions, respectively). Also, for a very thin sheet, it is reasonable to assume that a state of plane r stress exists in the sheet, i.e. 3 3 0 (where subscript 3 corresponds to the thickness direction). Furthermore, by adopting Kirchho's assumption that plane sections remain plane and normal to the middle surface of the sheet, shear strains in the thickness direction can be neglected; i.e. "13 "23 0. Then, by using the above assumptions, the yield function in Eq. (5) can be simpli®ed to the following expression: 22 22 R 1 ÿ 2 2 ÿ R 1 2 0
6
where 1 and 2 are principal stresses which are the radial and the circumferential inplane stresses, respectively. ÿ The strain rate tensor, D dij , simpli®es to only having three principal components, d1 d11 , d2 d22~,and d3 d33 . But, only two of these components are independent, since the strain rate in the thickness direction, d3, can be found as a r function of d1 and d2. By writing Eq. (1) in the indicial form and then setting 3 0 (plane stress) it is possible to obtain the plane stress moduli and an expression for d3 as a function of d1 and d2 (Osias, 1972). The ®nal form of Eq. (1) after the simpli®cation becomes: r
K d
7
where K is the plane stress elastic±plastic moduli and the Greek symbols vary between 1 and 2.
682
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
4. Kinematics of axisymmetric deformation The deformation of the mid-surface of an element in the plane containing the meridian curve will be considered for the axisymmetric case. It is noted that what follows is developed for the analytical portion of the hybrid membrane/shell method and that the treatment of kinematics in ABAQUS or other FEM code is independent. The deformation of the element in the circumferential direction will be obtained from the theory of shells of revolution (Flugge, 1973) and Figs. 1±3 show the midsurface of an element of the sheet at the initial time o , reference time o and at current time as it bends and stretches. Wang and Tang (1986) brie¯y discussed the kinematics of deformation for the case of axisymmetric, based on total Lagrangian formulation (Fig. 3). Choudhry and Lee (1989, 1994) discussed the kinematics of deformation, based on the updated Lagrangian formulation for the case of plane strain (see Fig. 2). In this paper, derivation of principal curvatures, stretches and strains of a shell element undergoing an axisymmetric deformation will be discussed in detail using both total and updated Lagrangian formulations. 4.1. Principal stretches (total Lagrangian formulation) Fig. 3 shows a ¯at element at the initial time o and at the current time o as it bends and stretches. By referring to Fig. 3, the following relationship between the current and the initial position vectors and the total displacement vector for the mid-surface of the element can be written:
Fig. 1. A schematic of a shell of revolution (Flugge, 1973).
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
683
Fig. 2. Shell mid-surface at reference time o and current time (updated Lagrangian).
Fig. 3. Shell element at original time o and current time (total Lagrangian).
x XU ~ ~ ~
8
where the initial position vector X and the total displacement vector U are, de®ned ~ in the global coordinate system, ~
684
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
X1 ; X 0 ~
u U w ~
Therefore, the position vector x becomes: ~ X1 u x w
9
10
Consider the unit tangent vector, a^, of the current mid-surface to calculate the ~ i.e principal stretch in the radial direction; dx dx dS dx 1 1 u;X1 1 ~ ~ ~ a^ w;X1 ds dS ds dX1 lo1 lo1 ~
11
where lo1 ds=dS. Here, s and S are travel lengths of the current and initial surfaces from center, respectively, and the comma stands for the dierentiation with respect to X1. Then, by calculating the magnitude of vector a^ from Eq. (11) and setting it ~ equal to unity, lo1 is
q
a^ a^ a^ 1:0 12
~ ~ ~ lo1
hÿ
i12 2 1 u;X1 w;2X1
13
The principal stretch in the circumferential direction, lo2 , is calculated from the change in the radial position of the element; i.e. lo2
X1 u u 1 X1 X1
14
4.2. Principal curvatures and stretches (updated Lagrangian) After bending and stretching, the principal mid-surface curvatures, k1 and k2, of a shell element at the current con®guration, o , can be calculated from the known information about the element at the reference con®guration (time o , see Fig. 2), as follows: r R uA^ wN^
15
where R is the reference con®guration at time o and u and w are incremental dis~ placements de®ned in Fig. 2. In Eq. (15), the unit tangent vector A^ and the unit ~ principal normal vector N^ to the mid-surface of the reference are de®ned as ~
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
685
@R A^ ~ R;s @S ~ ~
16
A;s N^ ~ K1 ~
17
where K1 is the curvature. Therefore, N^ :A^ 0 ~ ~
18
ÿK1 A^ N^ ;s ~ ~
19
and
To calculate the current curvature, k1, the unit tangent vector, a^ a= a and ~ ~ ~
the unit principal normal vector of the mid-surface of the shell, n^ n= n , are ~ ~ ~ calculated. Using Eq. (15), the tangent vector a is found as: ~ @r 20 a ~ r;s R;s u;s A^ uA^ ;s w;s N^ wN^ ;s ~ ~ ~ ~ ~ ~ @S ~ By substituting from Eqs. (16), (17) and (19) into (20), and re-arranging, the following expression results: @r ÿ ÿ a ~ r;s 1 u;s ÿ K1 w A^ w;s K1 u N^ cA^ dN^ ~ ~ ~ ~ ~ @S ~
21
The principal incremental stretch of the mid-surface in the radial direction is calculated from the magnitude of the base vector in Eq. (21):
q p 2 2
22 l1
a~ a~ :a~ c d The current length of the mid-surface of the shell in the radial direction, ds, is calculated from the reference length, dS, and l1 as follows: ds l1 dS The unit principal normal vector of the surface of the current shell, n^, is ~ ^ ^ a^ ;s a^ ;s ÿdA cN ~ ~ n^ ~ ~ l1 ~ a~^ ;s a~^ ;s
23
24
686
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
which, from Eqs. (21), (22) and (24), shows that n^:a^ 0. The current principal cur~ ~ vature of the shell, k1, can now be found as: 1 1 k1 ÿa^:n^ ;s ÿr;s :n^ ;s ÿ 2 r;s :n^ ;s ÿ 2 a:n^ ;s ~ ~ ~ ~ l1 ~ ~ l1 ~ ~ where a is given by Eq. (20) and n^ ;S can be derived from Eq. (24): ~ ~ dn^ l1 : ÿd;s A^ c;s N^ ÿ K1 dN^ ÿ K1 cA^ ÿ l1;S : ÿdA^ cN^ ~ ~ ~ ~ ~ n^ ;s ~ dS ~ l2
25
26
1
In Eq. (26), l1;S is assumed to vanish within an element and the above expression simpli®es as: n^ ;S ~
ÿ ÿ d;S K1 c A^ ÿc;s K1 d N^ ~ ~ ÿ l1
27
By substituting from Eqs. (27), (20) and (20) into Eq. (25), the current centerline curvature of the shell, 0, can be found: k1
cd;s ÿ dc;s K1 l21 l31
28
It should be noted that it is possible to recover Eq. (13), for total Lagrangian formulation, from Eq. (22). This is done by setting the reference curvature K1 equal to zero for a ¯at sheet, replacing S with X1 and treating the incremental displacements u and w as total displacements in Eqs. (21) and (22). Similarly, the current curvature of the shell, k1, can be calculated from the initial con®guration by setting K1 equal to zero, replacing l1 with lo1 and S with X1 and treating u and w as total displacements in Eq. (28) ÿ ÿw;x1 u;x1 x1 1 u;x1 w;x1 x1 29 k1 3 lo1 The current principal stretch increment of the mid-surface of the shell in the circumferential direction, l2 , is found from the geometry (in Fig. 2) as: l2
r1 R1
30
The principal curvature of the shell in the circumferential direction, k2, is found from Fig. 1 (Flugge, 1973), as follows:
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
k2
1 sin 1 dw w;s ÿ ÿ ÿ r2 R R ds R
687
31
where w and s were de®ned previously and R is de®ned as (see Fig. 1): R X1 u
32
Since it is known that ds lo1 dX1 , Eq. (31) becomes k2 ÿ
1 dw 1 dw dX1 1 dw w;X ÿ ÿ o ÿ o 1 R ds R dX1 ds l1 R dX1 l1 R
33
Finally, by substituting from Eq. (32) into (33), we will get the following expression for the principal curvature in the circumferential direction: k2 ÿ
w;X1 X1 u:lo1
34
4.3. Strain rate The principal strain rates, d1 and d2 for any through-the-thickness surface are calculated from the principal stretches, l1 and l1, and curvatures, k1 and k2, of the mid-surface of the shell. The principal stretches of any through-the-thickness surfaces l 1; 2, are equal to: l z l 1 zk
35
and through-the-thickness principal strains are equal to: " z ln l z ln l ln 1 zk
36
The strain rates are then calculated by taking the time derivative of Eq. (36) as: h i : : : d ln l z z:k k :z l 37 d l 1 zk dt : where in Eqs. (35)±(37) takes on the values 1 and 2. In Eq. (37), l are the rate of principal stretches l1 and l2 given by Eqs. (22) and (30). The rate of the centerline : curvatures, k , are calculated as the dierence between the current and reference values: : k ÿ ok k
38
688
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
In Eq. (37), z is the current location of a through-the-thickness surface, originally a distance away from the mid-surface. Based on the plastic incompressibility condition and ignoring elastic volume change, the following relationship between z and for each through-the-thickness layer can be obtained: z2 z3 l1 :l2 : z k1 k2 : k1 :k2 : 2 3
39
where ÿto =244to =2 and ÿt=24z4t=2 and parameters to and t are the original and the current thickness of the shell respectively. The current thickness of the shell is updated, after each deformation increment as follows: t ot ed3
40
where o t is the thickness of the shell at the reference state. 5. Stresses In a co-rotational coordinate system, the spins relative to the co-rotational frame are zero and the material derivative of stress is equal to: D r : ~ D ~ ~
41
r
where is given by Eq. (7). The updated current total stress is obtained by adding : ~ the incremental stress, : ,1 to the stress tensor at the reference state, o as ~ ~ ~ follows: o ~ ~ ~
42
where is an increment of time and the reference state is at the beginning of the time increment. 6. The analytical hybrid membrane/shell method To calculate the bending deformation, the hybrid method utilizes the intermediate sheet shape information (x and y coordinates of each node in the sheet for every : : 1 Here we assume 1 =2 1 = 2 because the principal directions are materially ®xed for axisymmetric deformations and the incremental deformations and the incremental deformation theory is supplied.
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
689
increment) generated by the membrane code in a post processing fashion (Pourboghrat and Chu, 1995; Pourboghrat et al., 1998, 1999). The intermediate membrane shapes are used to calculate the bending strain increment for each integration point by an analytical model, as explained in Section 4. The bending stress increments are calculated for each integration point using the elasto-plastic constitutive equation described in Section 3. The total (bending + stretching) deformation at the end of each increment is obtained by integrating the stresses and strains over the increments, as described in sections 4 and 5. The pair of tension components calculated by integrating through-the-thickness hybrid shell stresses are T
t=2 ÿt=2
z:dz 1; 2
43
which are the same as the in-plane membrane forces. Eq. (43) decides the location of the neutral surface. The pair of bending moments that are generated by the hybrid shell stresses are M
=2 ÿ =2
z:z:dz
44
which are assumed to be balanced by additional external contact forces between the sheet and the tooling prior to springback calculation (the contact forces in the membrane analysis balance only tension-curvature forces). Springback is calculated with a shell model (e.g. ABAQUS, 1998), using the membrane shape, the hybrid shell stresses and the forming tools in their ®nal position. Prior to the unloading, the equilibrium of the hybrid shell stresses with the external contact forces are enforced (by ABAQUS), with minimal adjustments made by ABAQUS to the membrane solution. 7. Springback results and discussions In the punch-stretching test, a ¯at sheet is initially clamped between a ¯at blank holder (binder) and the bottom die with a speci®ed blank holding force, as shown in Fig. 4. Then, the hemispherical punch bends and stretches the sheet by pushing the sheet into the die cavity. In the punch-stretching test, the blank holding force is very large. There is very little draw-in and the sheet mostly stretches. Fig. 4 shows the simulated middle surface for the axisymmetric punch-stretching test using membrane element in the ®nite element code ABAQUS, (1998). Due to the symmetry, only half of the sheet was used for the ®nite element simulation. Table 1 shows the friction coecient and the material properties of the aluminum alloy and the DQAK steel sheet used in the ®nite element simulations. Table 2 shows the
690
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
tooling speci®cations, the sheet dimensions and the blank holding force used in the punch-stretching test (Stevenson, 1993). In Fig. 4, the deformed shapes of the sheet correspond to the punch displacements of 7.62, 12.70 and 20.3 mm, respectively. Fig. 5a and b, for the aluminum alloy sheet, shows the dierence between the
Table 1 Material properties of the sheet alloys and friction coecient used in modeling Material type
K (Mpa)
n
E (Gpa)
R
Y (Mpa)
"o
friction coecient
AL 2036-T4 DQAK Steel
604.0 550.0
0.214 0.248
69.0 210.0
0.84 2.10
196.0 188.0
0.33 0.30
0.0052 0.0132
0.1 0.1
Table 2 Tooling geometry and boundary conditions Case Sheet thickness Punch radius, Die radius, Die cavity Punch displacement Blank holding (mm) Rp (mm) Rd (mm) opening (mm) (mm) force (KN) A B C
0.68 0.68 0.68
101.6 101.6 101.6
15.0 15.0 15.0
105.0 105.0 105.0
7.6 12.7 20.3
150.0 150.0 150.0
Fig. 4. The simulated middle surface for the axisymmetric punch-stretching test using membrane element (ABAQUS). The deformed shapes correspond to punch depths of 7.62, 12.7 and 20.32 mm, respectively.
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
691
Fig. 5. Simulated axisymmetric cup stretch/drawing test and deformed shapes of the sheet at 7.62 and 20.32 mm punch depths, respectively (with ABAQUS code).
simulated sheet shapes by membrane and shell elements at punch displacements of 7.62 and 20.3 mm, respectively. Fig. 5a shows a small dierence in the two shapes around the die corner, while Fig. 5b shows no dierence in the two shapes due to increased tension in the sheet. The punch-stretching experiments performed by Stevenson (1993) to measure springback were done in two steps. In this paper, similar to the experiments, springback of the sheet was also calculated in two steps using both shell (with one integration point on the mid-surface and ®ve through the thickness) and the hybrid membrane/shell methods (with 11 through-the-thickness layers). In the ®rst step of the springback simulation, only the punch was removed and the blank holder was left engaged (used as a reference surface). The vertical displacement of the sheet at the punch pole was then measured. Figs. 6a±c and 7a±c show the comparison
692
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
Fig. 6. Loaded and unloaded sheet shapes (SPB1):hybrid vs. shell vs. experiment at punch displacements of (a) 7.6 mm (b) 12.7 mm and (c) 20.3 mm for aluminum alloy sheet.
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
693
Fig. 7. Loaded and unloaded sheet shapes (SPB1):hybrid vs. shell vs. experiment at punch displacements of (a) 7.6 mm (b) 12.7 mm and (c) 20.3 mm for DQAK steel sheet.
694
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
between the measured springback by Stevenson (1993), and numerically predicted springback (SPB1) with shell and hybrid methods for aluminum and steel alloy sheets, respectively. It is interesting to note, from Fig. 7c, that the steel sheet undergoes very little springback after being deformed by the punch to a depth of 20.3 mm. The reason for this is that the sheet is highly stretched and the residual stresses in the sheet are negligible. On the other hand, for the same deformation, the aluminum sheet springback by about 1mm (Fig. 6c). Both shell and hybrid models predict similar springback results. In the second step of the springback simulation, both the blank holder and the punch were removed and by holding the last node under the blank holder ®xed (used as a reference point), the sheet was allowed to come to a complete unloading. Figs. 8a±c and 9a±c show the comparison between the shell and the hybrid method's unloaded shape predictions (SPB2). The exact comparison of predicted unloaded sheet shape with the experimental results, for the second step (SPB2), was dicult since unlike the ®rst step (SPB1) there was no known ®xed reference surface for the experimental measurements. In the experiments (Stevenson, 1993), the maximum height in the deformed sheet was measured using the coordinate measuring machine (CMM). To compare predicted and experimental values, the maximum heights in the deformed sheet for the shell and the hybrid models were calculated from Figs. 8a±c and 9a±c, as Hmax= YmaxÿYmin, and then compared against the experimental data. The last three columns in Tables 3 and 4 show the comparison between the predicted and measured maximum height data for aluminum and DQAK steel sheets. 7.1. Shell versus hybrid Figs 6a±9a consistently show that: (1) the aluminum sheet has a larger springback than the steel sheet, (2) the unloaded shapes of the sheet predicted by the hybrid membrane/shell method, for the punch travel of 7.62 mm, are not in good agreement with those predicted by the conventional shell model. The observed dierences in the shapes are caused by a larger springback predicted around the die corner (at X67.5 mm, Figs. 8a and 9a) by the hybrid method. These discrepancies in predicted unloaded shapes, between the hybrid and the shell models, are substantially improved for the punch travel of 12.7 mm (see Figs. 8b and 9b) and are completely eliminated for the punch travel of 20.3mm (see Figs. 8c and 9c). By referring to Fig. 10a, it can be seen that around the die radius area (nodes 160± 170), at the punch displacement of 7.6 mm, the thickness strains ("t 0:1%) predicted by the hybrid method are smaller than the initial yield strain of the sheet ("y 0:28%). In this situation the hybrid method, unlike the shell model, predicts a larger springback due to elastic thickness strains at the die radius area (see Figs. 8a and 9a). But, as the magnitude of the thickness strain and geometric stiening increases, with the deeper punch penetrations, the discrepancy in the unloaded sheet shape between the hybrid and shell models disappears. Fig. 10c shows that for the aluminum sheet, at the punch displacement of 20.3 mm, the thickness strains are plastic ("t ÿ0:5%) around the die radius area and therefore both hybrid and shell models predict similar springback (see Figs. 8c and 9c).
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
695
Fig. 8. Loaded and unloaded sheet shapes (SPB2):hybrid vs. shell at punch displacements of (a) 7.6 mm (b) 12.7 mm and (c) 20.3 mm for aluminum alloy sheet.
696
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
Fig. 9. Loaded and unloaded sheet shapes (SPB2) hybrid vs. shell at punch displacements of (a) 7.6 mm (b) 12.7 mm and (c) 20.3 mm for DQAK steel sheet.
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
697
Table 3 Springback of aluminum 2036-T4 alloy sheet Max. punch height (mm)
SPB1 hybrid (mm)
SPB1 shell (mm)
SPB1 experiment (mm)
SPB2 hybrid Hmax (mm)
SPB2 shell Hmax (mm)
SPB2 experiment (mm)
ÿ7.62 ÿ12.70 ÿ20.30
ÿ5.69 ÿ11.37 ÿ19.39
ÿ5.51 ÿ11.30 ÿ19.36
ÿ5.26 ÿ10.87 ÿ18.82
ÿ4.04ÿ0.59=ÿ4.63 ÿ9.50ÿ0.59=ÿ10.09 ÿ18.26ÿ0.32=ÿ18.58
ÿ4.72ÿ0.01=ÿ4.73 ÿ9.66ÿ0.43=ÿ10.09 ÿ18.26ÿ0.29=ÿ18.55
ÿ5.23 ÿ11.20 ÿ19.28
Table 4 Springback of draw quality aluminum killed (DQAK) steel sheet Max. punch height (mm)
SPB1 hybrid (mm)
SPB1 shell (mm)
SPB1 experiment (mm)
SPB2 hybrid Hmax (mm)
SPB2 shell Hmax (mm)
SPB2 experiment (mm)
ÿ7.62 ÿ12.70 ÿ20.30
ÿ6.95 ÿ12.22 ÿ19.96
ÿ6.83 ÿ12.19 ÿ19.96
ÿ6.95 ÿ12.25 ÿ20.00
ÿ6.03ÿ0.45=ÿ6.48 ÿ11.56ÿ0.27 =ÿ11.83 ÿ19.58ÿ0.14=ÿ19.72
ÿ6.27ÿ0.20=ÿ6.47 ÿ11.58ÿ0.25=ÿ11.83 ÿ19.57ÿ0.15=ÿ19.72
ÿ6.39 ÿ11.70 ÿ19.45
Fig. 10b and d show that when the hybrid method uses the shell solution, the error in predicted thickness around the die corner diminishes. Therefore, the main source of error in springback prediction by the hybrid method comes from the use of the membrane solution. However, the magnitude of these errors is small (Fig. 8a and 9a) and it diminishes as the sheet further deforms and stretches (Figs. 8b and 9b and c). 7.2. Comparison of CPU time between shell and hybrid methods Table 5 shows a comparison between the total CPU seconds taken by the hybrid membrane/shell method and a full shell model to analyze the punch stretching and a two-step springback calculation problem. The hybrid method took about 30 CPU s, on a HP-755 workstation, to compute bending stresses for 100 increments of membrane solutions. Each membrane solution (i.e. sheet shape) was comprised of 200 membrane elements. The shell models, with one integration point on the mid-surface and both 5 and 11 integration points through the thickness, required a much larger CPU time to analyze the same problem, as shown in Table 5. The amount of the CPU time required for calculating springback was almost the same for both hybrid and shell models. This is because the hybrid method also used a shell model to perform the unloading steps. It can be seen that the hybrid membrane/shell model has at least a 50% advantage in CPU time over the conventional shell analyses for this simple problem.
698
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
Fig. 10. (a) Predicted sheet thickness:hybrid using membrane shapes vs. shell at punch displacement of 7.6 mm for aluminum alloy sheet; (b) predicted sheet thickness:hybrid using shell shapes vs. shell at punch displacement of 7.6 mm for aluminum alloy sheet; (c). Predicted sheet thickness:hybrid using membrane shapes vs. shell at punch displacement of 20.3 mm for aluminum alloy sheet; (d) Predicted sheet thickness:hybrid using shell shapes vs. shell at punch displacement of 20.3 mm for aluminum alloy sheet.
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
699
Table 5 Comparison of the CPU time for dierent methods Method
CPU (s)
Number of iterations
Hybrid membrane/shell Full shell, ®ve integration points Full shell, 11 integration points
821 1706 2173
171 191 230
8. Conclusions To reduce the cost of the ®nite element analyses for sheet forming, a hybrid membrane/shell method was developed to study the springback of anisotropic sheet metals undergoing axisymmetric loading. The goal was to apply the hybrid method incrementally in order to determine the bending strain prediction accurately. In this model, the incremental bending/unbending strains and stresses were calculated only from the information obtained from ABAQUS about the shape of the membrane sheet. An ABAQUS model with shell elements was used to unload the ®nal shape of the sheet with initial stresses and strains that were calculated analytically. For veri®cation, the hybrid method was applied to predict springback in an axisymmetric punch-stretching test for 2036-T4 aluminum and DQAK steel materials. The hybrid membrane/shell springback predictions were in good agreement with the ABAQUS full shell analysis. It was observed that there are some dierences between the experimental springback data and those predicted by the numerical models. These discrepancies were attributed to neglecting to model planar anisotropy, particularly for aluminum 2036-T4 alloy, and the possible Bauschinger eect. It is intended to apply the hybrid membrane/shell method to three-dimensional sheet forming. In that case, the eect of modeling planar anisotropy and the Bauschinger eect on predicted springback results for punch stretching will be assessed once again. A minimum of 50% CPU time saving was achieved by using the hybrid method to predict the springback of an axisymmetric cup instead of a full shell model. It is expected to save even more CPU time for the three dimensional problems. Acknowledgements The authors wish to sincerely thank Dr. Owen Richmond of Alcoa for the support they have received during the course of this work. Farhang Pourboghrat would like also to thank Dr. Robin Stevenson of the Engineering Mechanics Department at GM for his helpful discussions and providing experimental samples and data for this work. References ABAQUS, 1998. User's and Theory Manuals, version 5.8. Hibbitt Karlsson Sorensen Inc, RI. Barlat, F., Maeda, Y., Chung, K., Yanagawa, M., Brem, J.C., Hayashida, Y., Lege, D.J., Matsui, K., Murtha, S.J., Hattori, S., Becker, R.C., Makosey, S., 1997. Yield function development for aluminum alloy sheets. J. Mech. Phys. Solids 45, 1727.
700
F. Pourboghrat et al. / International Journal of Plasticity 16 (2000) 677±700
Becker, R.C., 1992. An ecient semi-implicit integration algorithm for anisotropic ¯ow potentials, unpublished work. Choudhry, S., Lee, J.K., 1989, Numerical simulation of thin sheet metal forming processes including bending eects, In: Durham, D., Saigal, A. (Eds.), Proc. microstructural development and control in materials processing (ASME/WAM, ASME MD-Vol. 14), San Francisco, CA, pp. 11±18. Choudhry, S., Lee, J.K., 1994. Dynamic plane-strain ®nite element simulation of industrial sheet-metal forming processes, Int. J. Mech. Sci. 36 (3), 189±207. Chung, K., Richmond, O., 1993. A deformation theory of plasticity based on minimum work path. Int. J. Plasticity. Chung, K., Richmond, O., 1994. The mechanics of ideal forming. J. of Applied Mechanics 61, 176±181. Esche, S.K., Kinzel, G.K., Lee, J.K., 1999. An axisymmetric membrane element with bending stiness for static implicit sheet metal forming simulation. ASME Journal of Applied Mechanics 66, 153±164. Flugge, W., 1973. Stresses in Shells, 2nd Edition. Springer-Verlag, New York. Frey, W.H., Wenner, M.L., 1987. Development and application of a one-dimensional ®nite element code for sheet metal forming analysis. In: Samanta, S.K. et al. (Eds.), Interdisciplinary Issues in Materials Processing and Manufacturing, Vol. 1. Hill, R., 1948. Theory of the yielding and plastic ¯ow of anisotropic metals. Royal Soc. London Proc. 193A, 281. Kara®llis, A.P., Boyce, M.C., 1992. Tooling design in sheet metal forming using springback calculations. Int. J. Mech. Sci. 34 (2), 113±131. Kubli, W., Anderheggen, E., Reissner, J., 1991. Nonlinear solver with uncoupled bending and stretching deformation for simulating thin sheet metal forming. Proc. FE-Simulation of 3-D Sheet Metal Forming Processes In Automotive Industry, 14±16 May, Zurich, Switzerland. Nadai, A., 1963. Theory of ¯ow and fracture of solids, Vol. 2. McGraw-Hill, New York, p.74. Osias, J., 1972. Finite deformation of elasto-plastic solids. Ph.D. thesis, Carnegie Mellon University. Pourboghrat, F., Chandorkar, K., 1992. Springback calculation for plane strain sheet forming using ®nite element membrane solution. Symposium on Numerical Methods for Simulation of Industrial Metal Forming Processes, ASME WAM'92, 8±13 November, Anaheim. Pourboghrat, F., Chu, E., 1995. Springback in plane strain stretch/draw sheet forming. Int. J. Mech. Sci. 36 (3), 327±341. Pourboghrat, F., Chung, K., Richmond, O., 1998. A hybrid membrane/shell method for rapid estimation of springback in anisotropic sheet metals. ASME Journal of Applied Mechanics 65 (3), 671±684. Pourboghrat, F., Yoon, J.W., Chung, K., Barlat, F., 1999. A 3D hybrid membrane/shell method to predict the springback of anisotropic sheet metals. The 7th International Symposium on Plasticity and Its Current Applications, 5±13 January 1999, Cancun, Mexico. Saran, M.J., Wagoner, R.H., 1991. A user's manual and veri®cation examples for section form (Sheet-S/ version 3), section analysis of sheet stampings. Report no. ERC/NSM-S-91-17. Sklad, M.P., Siekirk, J.F., 1990. Proc. 16th IDDRG Congress, ASM International, Bolange, Sweden, pp. 295±303. Stevenson, R., 1993. Springback in simple axisymmetric stampings. Metallurgical Transactions 24A, 925±933. Stoughton, T.B., 1985. Finite element modeling of 1008 ak sheet steel stretched over a rectangular punch with bending eects. In: Wang, N.M., Tang, S.C. (Eds.), Computer Modeling Of Sheet Metal Forming Processes. The Metallurgical Society, Warrendale, PA, pp. 143±159. Wagoner, R.H., Zhou, D., Sriram, S., 1995. The element bending group method. Proc. International Conference on Computational Engineering Science, 30 July±3 August, HI, USA. Wang, N.M., Tang, S.C., 1986. Analysis of bending eects in sheet forming operations. Int. J. Numer. Methods Eng. 25, 253±267. Wenner, M.L., 1992. Elementary solutions and process sensitivities for plane-strain sheet-metal forming. J. of Applied Mechanics 59, 23±28.