Plane strain transversely anisotropic analysis in sheet metal forming simulation using 6-component Barlat yield function

In most FEM codes, the isotropic-elastic and transversely anisotropic-elastoplastic model using Hill’s yield function has been widely adopted in 3D shell elements (modified to meet the plane stress condition) and 3D solid elements. However, when the 4-node quadrilateral plane strain or axisymmetric element is used for 2D sheet metal forming simulation, the above transversely anisotropic Hill model is not available in some FEM code like Ls-Dyna. A novel approach for explicit analysis of transversely anisotropic 2D sheet metal forming using 6-component Barlat yield function is elaborated in detail in this paper, the related formula between the material anisotropic coefficients in Barlat yield function and the Lankford parameters are derived directly. Numerical 2D results obtained from the novel approach fit well with the 3D solution.


Introduction
In order to accurately simulate sheet metal forming processes, it is essential to describe correctly the material constitutive behaviors. Since most sheet metals exhibit anisotropic material behaviors, the use of appropriate anisotropic yield criterion is important to predict material behaviors accurately. Moreover, anisotropy has an important effect on the strain distribution in sheet metal forming process, and it is closely related to thinning and formability of sheet metal, so the anisotropy of the material should be properly considered to capture the realistic material behaviors.
The influence of plastic anisotropy on sheet metal forming has been studied with the help of FEM codes combined with appropriate anisotropic yield functions. Many such functions have been proposed. The quadratic yield function by Hill (1948) has long been one of the popular choices to represent planar anisotropy and has been widely used in FEM forming simulation. Several non-quadratic criteria were developed by Hill (1979Hill ( , 1990, Hershey (1954), Hosford (1972), Bassani (1977), Gotoh (1977), Logan and Hosford (1980), Barlat and Lian (1989), Karafillis and Boyce (1993), Bron and Besson (2004), Banabic et al. (2005) and Barlat et al. (1991Barlat et al. ( , 1997Barlat et al. ( , 2003Barlat et al. ( , 2005. In general, Hill (1948) has been useful for explaining phenomena associated to anisotropic plasticity particularly for steels, the others can be used to improve the yielding description of aluminum alloys. But in many circumstances Yld89 (Barlat and Lian 1989), Yld91 (Barlat et al. 1991) can be used for steels or aluminum alloys. In particular, the yield criteria Yld89 for planar anisotropy, and Hill (1990) have three stress components and are applicable to plane stress condition. The analytical forms of these criteria are relative simple. The criteria Yld91 account for six stress components and can be applied to general 3D elasto-plastic continuum codes. Barlat et al. (2007) have made detailed discussion about the relation between these yield criteria. In fact, Yld91 is a particular case of Yld2004-18p (Barlat et al. 2005), Yld89 is a particular case of Yld2000-2d (Barlat et al. 2003), the yield function proposed by Banabic et al. (2005) is identical to Yld2000-2d. In some special conditions, Yld91 can reduce to Hill (1948), Mises or Tresca yield function.
Section analysis provides a faster and more efficient alternative procedure for analyzing complex part shapes in some cases. In this procedure, a cross section is selected from the tooling along a direction of interest. The problem is then analyzed in 2D, assuming plane strain or axisymmetric conditions. Although most sections in sheet forming do not completely satisfy such assumptions, there are still many local sections in a complicated sheet forming process which can be successfully simulated by 2D section analysis. In general, this method can provide a quicker analysis and designers can modify local geometry of tools with experience at preliminary design stages.
The present study attempts to conduct an explicit section analysis of 2D sheet metal forming using the 4-node quadrilateral plane strain element. In most FEM codes, the isotropic-elastic, transversely anisotropic-elastoplastic model proposed by Hill (1948) has been widely used in 3D shell elements (modified to meet the plane-stress condition) and 3D solid elements. However, the above transversely anisotropic model is not available for simulating 2D sheet metal forming process using 4-node quadrilateral plane strain or axisymmetric element in some FEM code like Ls-Dyna (Hallquist 2007). In this paper, a novel approach for the explicit analysis of transversely anisotropic 2D sheet metal forming using 6-component Barlat yield function (Barlat et al. 1991) is proposed, the related formula and parameters are derived directly in detail, the numerical 2D results obtained from the novel approach fit well with the 3D solution.
2 Fundamental theory

Basic finite element equation
The sheet metal forming process can be treated as a dynamic contacting problem. In the forming process, the whole system should satisfy the following finite element equation: where F ext is the vector of external force, F c the vector of contact force, F int the vector of internal force, and M the mass matrix. In the explicit algorithm, Eq. (1) can be solved as follows: If a concentrated mass matrix is assumed, it is not necessary to solve the assembled equations since the mass matrix becomes diagonal. The acceleration vector € U n at t n can be solved by Eq.

The general Hill orthotropic anisotropic yield criteria
In 1948, Hill proposed a orthotropic anisotropic yield function (Hill 1948) according to Mises criteria: where F, Q, H, L, M and N are anisotropic constants relating with the material yield behaviors. When F = Q = H = L/3 = M/3 = N/3, (5) reduces to Mises criteria. For the general anisotropic material behaviors, they meets the follow condition: where r s 1 ; r s 2 ; r s 3 ; r s 23 ; r s 31 and r s 12 are the tension or shear yield stresses in corresponding directions.

Transversely anisotropic criteria
For the plane strain problems, Hill's orthotropic yield function can be simplified: For transversely anisotropic sheet material, The following relation can be derived from (6).
The Lankford parameter is determined as follows:
After some manipulations, Eq. (20) can be described explicitly in details, given in Appendix.

The relation between material anisotropic coefficients
In general, the Lankford parameters R in three directions (i.e., 0°, 45°and 90°) should be available when sheet metals are provided. On the other hand, the six material anisotropic coefficients a, b, c, f, g and h are generally not available. In order to use the above derived model, it is necessary to obtain their quantitative values from the Lankford parameters. The arbitrary Lankford parameters R / can be described as: where _ e p t is the plastic strain rate which is perpendicular to tensile axis of the uniaxial tensile sheet metal specimen, _ e p z the plastic strain rate which is along the thickness direction of sheet metal, and / the angle between the tensile axis and the rolling direction.
The stress status of a body unit in the uniaxial tensile sheet metal specimen is shown in Fig. 1. From Fig. 1 it can be obtained that: r xy ¼ r 11 sin / cos /; r xx ¼ r 11 cos 2 /; r yy ¼ r 11 sin 2 / _ e p 11 ¼ cos 2 /_ e p xx þ sin 2 /_ e p yy þ 2 sin / cos /_ e p xy ; _ e p 22 ¼ sin 2 /_ e p xx þ cos 2 /_ e p yy À 2 sin / cos /_ e p xy Assuming that the sheet metal material is incompressible, the Lankford parameter can be derived: Substituting formula (22) into Eq. (23), it can be obtained: According to the related flow criteria, the relation between R / and the yield function can be derived: Applying formula (20) In order to obtain the relation explicitly between R / and the six material anisotropic coefficients a, b, c, f, g and h, for simplicity and not losing generality, it can be obtained by assuming m = 2: W ¼ 4 cos 2 2h þ p 6 þ 4 cos 2 2h À 3p 6 þ 4 cos 2 2h þ 5p 6 ¼ 6 or xy sin 2/ À oI 2 or xx sin 2 / À oI 2 or yy cos 2 / oI 2 Applying formula (12) and (22), it can be obtained: Substituting three angles (/ = 0°, 45°, 90°) into Eq. (25), it can be obtained: The three Eqs. (29)-(31), contain four unknown coefficients a, b, c and h. They can't be solved deterministically. However, it is noted that the coefficients f, g and h represent the influences of material anisotropy applying to the stress component (r yz , r zx , r xy ). For isotropy materials, the six coefficients a, b, c, f, g and h are all equal to 1. When the anisotropy of the sheet metal is described, it can approximately be made that f = g = h = 1. It is now possible to solve a, b and c in the assembled Eqs. (29)-(31). For transversely anisotropic sheet metal, it's easy to see that: The three coefficients a, b and c can be solved in the Eq. (32).

Numerical example
The tool geometry data of sheet metal drawing and bending is defined in Fig. 2, the plate is L 9 W 9 H = 210 9 40 9 0.65 mm, and the hoderforce is 32 kN. The material property of the sheet is the isotropic-elastic transversely anisotropic-elastoplastic, whose parameters are listed below: where l is Coulomb friction coefficient. For transversely anisotropic sheet metal, we assume f = g = h = 1 as discussed earlier. For BCC materials, when Yld91 is used, let m = 6, but the three coefficients a, b and c can still be solved from (32) for simplicity and almost not losing accuracy. For this example, they are a = b = 0.866, c = 1. 688. A 320 9 4, mesh of 320 elements along the length direction and 4 elements across the thickness are used. The plane strain element in Ls-Dyna code is used. The deformed shapes and distributions of equivalent stress at punch depth 60 mm are shown in Fig. 3. It can be seen from Fig. 3 that the peak equivalent stress appears near the corner of the die. In Fig. 4, the 3D solution using the same Barlat model and the 3D BT shell element model is provided by Ls-Dyna code. Figure 5 shows the punch force as a function of the time for the 2D and 3D results. From Figs. 3 and 4, it can be concluded that the peak equivalent stress and equivalent stress distribution of the 2D results are almost same as those of the 3D results. From Fig. 5, it shows that there is a good agreement of punch force between 2D and 3D results. Fig. 1 The stress status of a body unit in the uniaxial tensile sheet metal specimen Fig. 2 Initial sheet shape and tool data Plane strain transversely anisotropic analysis 331 Therefore, the novel approach for transversely anisotropic section analysis is successful and valid.

Conclusions
The explicit analysis of transversely anisotropic plane strain sheet metal forming is implemented in Ls-Dyna code using 6-component Barlat yield function. The related formula and parameters has been derived explicitly. The transversely anisotropic property of sheet metal can be described by the three parameters a, b, and c. Numerical results demonstrated the validity of this approach.