Hole-Expansion: Sensitivity of Failure Prediction on Plastic Anisotropy Modeling

Abstract

The influence of yield function parameters on hole-expansion (HE) predictions are investigated for an anisotropic AA6022-T4 aluminum sheet. The HE experiment is performed in a fully-instrumented double-action hydraulic press with a flat-headed punch. Full strain fields are measured by a stereo-type digital image correlation (DIC) system. The stress state gradually changes from uniaxial to plane-strain tension to biaxial tension in the radial direction. Besides HE, the plastic anisotropy of AA6022-T4 is characterized by uniaxial tension and plane-strain tension experiments. Uniaxial tension is considered as the most important, since it is the stress state along the hoop direction in the hole. For the finite element (FE) simulation, the Yld2000-2d non-quadratic anisotropic yield function is used with two different parameter sets, calibrated by: (1) uniaxial tension only (termed Calib1) and, (2) both uniaxial and plane-strain tension (Calib2). The strain field predictions show a good agreement with the experiments only for Calib2, which takes into account plane-strain as well uniaxial tension. This indicates the importance of biaxial modes, and in particular plane-strain tension, for the adopted yield function to produce accurate HE simulations.

1. Introduction

The hole-expansion (HE) of a thin sheet using a flat-headed punch is akin to the expansion of a circular hole in an equibiaxial tensile field. The stress state at the hole periphery is uniaxial tension in a local coordinate system tangent to the hole, while farther away it is equibiaxial tension, with all the intermediate states occurring in-between. Furthermore, if the sheet is anisotropic, the strain fields around the expanding hole clearly reflect the corresponding orientational dependence of the material. By monitoring the evolution of these fields, e.g., using stereo-type digital image correlation (DIC) [1,2], the HE experiment can be used to probe the anisotropy of the sheet. Inversely, it can be used to validate the structure and calibration of an anisotropic yield function that has been selected to model the sheet behavior [2,3].

The isotropic HE case was treated analytically by Taylor [4] and Budiansky and Mangasarian [5]. Perhaps the first investigation of the anisotropic case was by [6], who used the yield function proposed by Parmar and Mellor [7] and discovered that the strain field predictions were very sensitive to the anisotropic parameters. Further solutions were provided by Durban [8,9,10] for the case of normal anisotropy (i.e., isotropy in the plane of the sheet). The effect of two non-quadratic yield criteria [11,12] was discussed by Cohen et al. [13].

Beyond these works, a variety of anisotropic yield criteria, e.g., Hill ’48 [14], Barlat ’89 [15], and Yld2000-2D [16], or crystal plasticity, have been used in numerical simulations of HE [3,17,18,19,20,21,22,23,24]. In some of these works, anisotropy was calibrated based on uniaxial tension tests at angles to the rolling direction (RD). When rotated to the material axes, these tests included a shear stress component, but they were unsuitable in probing the specifics of the yield locus in the plane-strain and equibiaxial tension regions, which, as discussed earlier, are dominant in the majority of the HE specimen.

In this paper, the effect of the calibration of the anisotropic yield function on the HE strain field predictions is examined. The material selected is AA6022-T4, as in our earlier works [2]. Two alternative data-sets for calibration are adopted, leading to two distinct representations of the yield locus. These are then introduced in finite element (FE) simulations of HE. The predictions of global (i.e., punch force-displacement) and local (i.e., strain fields and their evolution with punch displacement) responses are compared to the experiments.

The results of this work can aid industry in selecting and calibrating an anisotropic yield function appropriate to each material at hand. In particular, comparing the predictions to the HE experiment reveals the ability of the yield function to capture the resulting strain fields, including their local features. This is in turn necessary for accurate predictions of failure in the form of necking and the resulting tearing, which is a prerequisite for reliable virtual process design. Furthermore, this work can be used as a 1st step in assessing the forming limits that can be expected in hole stretching and flanging.

2. Plastic Anisotropy

2.1. Review of Anisotropic Yield Function

The non-quadratic anisotropic yield function Yld2000-2d [16] is constructed based on the non-quadratic isotropic one [11] formulated as:

$${\left|s_1-s_2\right|}^m + {\left|2s_1+s_2\right|}^m + {\left|s_1+2s_2\right|}^m=2{\overline{\sigma }}^m$$

(1)

where $s_{\mathrm{i}=1~3}$ are the principal components of deviatoric stress $\mathit{s}$ of Cauchy stress $\mathit{\sigma }$, $\overline{\sigma }$ is equivalent stress, and $m$ is an exponent based on the crystal structure, i.e., $m$ = 6 for body-centered cubic (BCC) and $m$ = 8 for face-centered cubic (FCC) materials [25]. To introduce plastic anisotropy in the isotropic function, the Cauchy stress $\mathit{\sigma }$ is linearly transformed into ${\mathit{S}}^{\mathbf{\prime }}$ and ${\mathit{S}}^{\mathbf{″}}$ as follows:

$${\mathit{S}}^{\mathbf{\prime }}={\mathit{C}}^{\mathbf{\prime }}\mathit{s}={\mathit{C}}^{\mathbf{\prime }}\mathit{T}\mathit{\sigma }={\mathit{L}}^{\mathbf{\prime }}\mathit{\sigma } \mathrm{and} {\mathit{S}}^{\mathbf{″}}={\mathit{C}}^{\mathbf{″}}\mathit{s}={\mathit{C}}^{\mathbf{″}}\mathit{T}\mathit{\sigma }={\mathit{L}}^{\mathbf{″}}\mathit{\sigma }$$

(2)

by two transformation operators ${\mathit{L}}^{\mathbf{\prime }}$ and ${\mathit{L}}^{\mathbf{″}}$. In plane-stress components form, Equation (2) can be written as:

$$\left\{\begin{array}{c}{S}_x^{\prime }\\ S_y^{\prime }\\ S_{xy}^{\prime }\end{array}\right\}= \left[\begin{array}{ccc}{L}_{11}^{\prime }& L_{12}^{\prime }& 0\\ L_{21}^{\prime }& L_{22}^{\prime }& 0\\ 0& 0& L_{66}^{\prime }\end{array}\right] \left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\sigma }_{xy}\end{array}\right\} \mathrm{and} \left\{\begin{array}{c}{S}_x^{″}\\ S_y^{″}\\ S_{xy}^{″}\end{array}\right\}= \left[\begin{array}{ccc}{L}_{11}^{″}& L_{12}^{″}& 0\\ L_{21}^{″}& L_{22}^{″}& 0\\ 0& 0& L_{66}^{″}\end{array}\right] \left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\sigma }_{xy}\end{array}\right\}$$

(3)

Each operator ${\mathit{L}}^{\mathbf{\prime }}$ and ${\mathit{L}}^{\mathbf{″}}$ is determined by anisotropic material parameters, i.e., ${\alpha }_{i=1~8}$, as:

$$\left\{\begin{array}{c}{L}_{11}^{\prime }\\ L_{12}^{\prime }\\ \begin{array}{c}{L}_{21}^{\prime }\\ L_{22}^{\prime }\\ L_{66}^{\prime }\end{array}\end{array}\right\}= \left[\begin{array}{ccc} 2/3& 0& 0\\ -1/3& 0& 0\\ 0& -1/3& 0\\ 0& 2/3& 0\\ 0& 0& 1\end{array}\right] \left\{\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ {\alpha }_7\end{array}\right\} \mathrm{and}$$

(4a)

$$\left\{\begin{array}{c}{L}_{11}^{″}\\ L_{12}^{″}\\ \begin{array}{c}{L}_{21}^{″}\\ L_{22}^{″}\\ L_{66}^{″}\end{array}\end{array}\right\}=\frac{1}{9} \left[\begin{array}{ccccc}-2& 2& 8& -2& 0\\ 1& -4& -4& 4& 0\\ 4& -4& -4& 1& 0\\ -2& 3& 2& -2& 0\\ 0& 0& 1& 0& 9\end{array}\right] \left\{\begin{array}{c}{\alpha }_3\\ {\alpha }_4\\ {\alpha }_5\\ {\alpha }_6\\ {\alpha }_8\end{array}\right\}$$

(4b)

As a result, by the transformation of Equation (1), Yld2000-2d anisotropic yield function is expressed as follows:

$${\left|S_1^{\prime }-S_2^{\prime }\right|}^m + {\left|2S_1^{″}+S_2^{″}\right|}^m + {\left|S_1^{″}+2S_2^{″}\right|}^m=2{\overline{\sigma }}^m$$

(5)

where $S_1^{\prime }, S_2^{\prime }, S_1^{″},S_2^{″}$ are the principal stress components of ${\mathit{S}}^{\mathbf{\prime }}$ and ${\mathit{S}}^{\mathbf{″}}$.

2.2. Calibration of Yld2000-2d Parameters

As in our earlier work [2], the material of this study is the automotive aluminum alloy AA6022-T4, received as a sheet of 1 mm thickness. Parameter calibration of Yld2000-2d can involve various experimental results, such as uniaxial tension, plane-strain tension, simple shear, in-plane biaxial tension, disk compression, etc. Depending on the applications of interest, different experiments can be used for the calibration, which results in different anisotropic parameters and yield loci. For example, in the original paper [16], the parameters were calibrated based on uniaxial tension in every 45° from the RD and equibiaxial tension. Four flow stresses and four strain ratios, in total eight experimental data, were used to calibrate the eight material parameters ${\alpha }_{i=1~8}$ using Newton–Raphson method. Ha et al. [26,27] examined the ductile fracture of automotive aluminum sheets, for which plane-strain tension was taken into account in the parameter calibration to capture the deformation in a localized zone before fracture. In HE [2], although the stress state varies from uniaxial to equibiaxial tension along the radial direction from the edge, anisotropy in uniaxial tension is often considered the most critical for parameter calibration due to the high strain concentration near the hole edge. As a result, the effect of various biaxial stress states in the surrounding material, which potentially lead to fracture, is often ignored in the material modeling and analysis of the numerical simulation.

In this work, the role of plastic anisotropy in biaxial tension, especially plane-strain tension, on the HE simulation is investigated by using two different sets of Yld2000-2d parameters. Each parameter set is calibrated with different experiments: Calib1 is based on uniaxial tension in every 15° from the RD, and Calib2 includes plane-strain tension in every 45° from the RD in addition to the same experiments for Calib1. Experimental data for elasto-plastic behavior is taken from the authors’ previous work [2]. The values are summarized in

Table 1. Summary of mechanical properties [2].

Elastic Properties
Young’s Modulus	E = 70 GPa		Poisson’s Ratio			ν = 0.3
Anisotropic Plastic Properties (Wp = 20 MJ/m3)
Uniaxial tension		RD		30°	60°		90°
	$\sigma /\overline{\sigma }$	1.000		1.016	0.976		0.954
	r-value	0.793		0.465	0.352		0.510
Plane-strain tension		RD		45°	90°
	$\sigma /\overline{\sigma }$	1.096		1.019	0.978

. Based on them, anisotropic parameters of Yld2000-2d are calibrated independently. The eight anisotropic parameters ${\alpha }_{i=1~8}$ are determined using the Levenberg–Marquardt nonlinear least-square algorithm. For the plane-strain tension, only the major stress component along the loading direction can be identified in the experiment (using the procedure detailed in [28,29,30,31]) and the minor one, transverse to the loading direction, is determined by the optimization with ${\alpha }_{i=1~8}$, with the plane-strain assumption (i.e., zero minor strain, transverse to the loading) [28,29,30]. The calibrated parameters for Calib1 and Calib2 are listed in

Table 2. Yld2000-2D parameters for Calib1 and Calib2.

m = 8	${\mathsf{\alpha }}_{\mathbf{1}}$	${\mathsf{\alpha }}_{\mathbf{2}}$	${\mathsf{\alpha }}_{\mathbf{3}}$	${\mathsf{\alpha }}_{\mathbf{4}}$	${\mathsf{\alpha }}_{\mathbf{5}}$	${\mathsf{\alpha }}_{\mathbf{6}}$	${\mathsf{\alpha }}_{\mathbf{7}}$	${\mathsf{\alpha }}_{\mathbf{8}}$
Calib1	0.970	1.054	1.253	1.128	1.065	1.253	0.940	0.909
Calib2	0.968	1.022	1.067	1.091	1.014	0.977	0.914	1.064

.

For the strain hardening, the same parameters as the Voce isotropic model are used with [2], i.e., $\overline{\mathsf{\sigma }}=430.52-213.71·\exp (-8.88·\overline{\mathsf{\epsilon }})$. As detailed, that hardening curve was identified by a special procedure, i.e., by matching the internal and external work inside the diffuse neck of a tensile specimen [32,33].

Figure 1 shows a comparison of the experiment and the prediction of Yld2000-2d for the normalized flow stresses and r-values in the uniaxial tension in every 15°. Both parameter sets, i.e., Calib1 and Calib2, capture the plastic anisotropy in the uniaxial tension very well and their predictions are almost identical. This is a natural output of the parameter calibration since both parameter sets are calibrated using the uniaxial tension.

Nevertheless, the corresponding yield loci of Calb1 and Calb2 show a significant difference in Figure 2, especially near plane-strain and biaxial stress states. This is because the experiment of plane-strain tension is not considered in the parameter calibration of Calib1 and, as a result, the plane-strain condition is predicted by Calib1 to be far away from the experiment. It should be noted that the plane-strain condition of Calib1 (blue empty square symbol) is located much closer to the uniaxial tension, while Calib2 is located near to the equibiaxial tension (red filled square symbol).

A quantitative assessment of the performance of Calib1 and Calib2 in capturing all the experiments is shown in Figure 3, using the representation proposed by [34]. See also [27,29,30,35]. For each experiment, the orientation of the stress state with respect to a reference direction is computed as $\cos \mathsf{\omega }$. Then, the normalized stress with respect to the equivalent stress (Figure 3a) and the angle between the normal of the predicted yield locus and the experiment (Figure 3b) are plotted against $\cos \mathsf{\omega }$. The deviations of these quantities from 1 are a quantitative measure of the ability by the given yield function to capture the experiments. As a comparison, von Mises is also evaluated with Yld2000-2d (Calib1 and Calib2); the latter is consistently much better than the former for all experiments considered in this work, see Figure 3. Both Calib1 and Calib2 of Yld2000-2d show a good agreement in the uniaxial tension condition, but Calib1 greatly misses the plane-strain tension condition.

3. Hole-Expansion

3.1. Experiment

The HE experiment was performed using a specimen with a circular hole of 35 mm diameter. The hole was prepared by end-milling to a good surface finish. The experiment was performed in a fully-instrumented hydraulic press of 260 kN capacity [2]. A flat-headed punch of 100 mm diameter with a 12 mm punch nose-radius expanded the hole as it traveled. To prevent the drawing of the flange and instead promote the expansion of the hole, a lock-ring was used. The tool design is shown in Figure 4.

A stereo-type DIC, i.e., 3D-DIC, system was used to measure the surface strain field throughout the process. Using incompressibility, the in-plane principal strains were then converted to thickness strains throughout the HE experiment. In particular, the thickness strain evolution was reported at different levels of punch displacement from a hoop with a 40 mm initial diameter, i.e., 2.5 mm inland from the hole edge, to avoid problematic correlation near the hole edge.

Detailed descriptions of the experimental set-up, procedure, and findings have been reported elsewhere [2]. Here, for completeness it suffices to summarize the main findings. The punch force-displacement was almost linear, with a sudden drop indicating the occurrence of rupture. The continuous monitoring of the specimen by the stereo-type DIC throughout HE revealed that the initially circular hole progressively became out-of-round, and that the strain fields varied along the circumference, reflecting the anisotropy of the sheet. To that effect, a two-fold symmetric pattern appeared early on and progressively intensified throughout the HE experiment. In this AA6022-T4 specimen, the greatest strain development occurred at around ±45° and ±135° from the RD. The maximum circumferential strain at the onset of rupture was 0.37, which is significantly higher that the uniform elongation in uniaxial tension (~0.18). The corresponding true thickness strain was −0.28 at that location. The reason that these strains were significantly larger than those experienced in uniaxial tension, even though uniaxial tension conditions prevail at the hole edge, is that in HE each concentric “ring” of the sheet is supported by another, less strained one, further away from the hole. The strain fields also varied in the radial direction, changing from uniaxial tension (at the hole circumference) to plane-strain to equibiaxial tension (at the punch nose-radius).

3.2. Finite Element Analysis

To investigate the effect of plastic anisotropy on the HE behavior, a FE simulation was performed using an Abaqus/Standard 2019 (implicit solver). The same tooling with the experiment (Figure 4), but in a quarter size, was modeled using non-deformable, analytical rigid body elements. The blank was constructed using deformable, four-node shell elements with reduced integration (S4R), for which one element every 1.5° was assigned along the hole circumference. Nine integration points were applied through the thickness direction based on a sensitivity evaluation [18]. The tools and the meshed blank used for the FE simulation are shown in Figure 5. Surface-to-surface contact was assumed for all contacts with the blank, and the Coulomb friction law was applied with a constant coefficient $\mu$ = 0.22 [2]. Constitutive models for elasto-plastic material behavior were implemented into a user-defined material subroutine (UMAT) [36], i.e., Hooke’s law, Yld2000-2d anisotropic yield function, and the Voce isotropic strain hardening law.

The simulation results with Calib1 and Calib2 were compared with the experiment regarding the punch force-displacement curve and the thickness strain variation around the hole. The former, force-displacement, is shown in Figure 6 with two experiments in black and red dots, with Calib1 in a blue solid line and Calib2 in a red solid line. The force level in Calib2 was slightly higher than in Calib1, but both predictions agreed reasonably well with the experimental force-displacement curve. This indicated that the plastic anisotropy did not have a significant effect on the structural behavior, such as force-displacement, in the HE simulation, and is in concert with our earlier findings in a variety of problems [2,28,37,38,39].

In contrast, the strain variation near the hole periphery showed a much higher sensitivity to the specifics of the description of plastic anisotropy. Figure 7a,b presents the thickness strain variation at different levels of punch displacement $\delta$, i.e., $\delta /{\delta }_{max}$ = 0.5, 0.65, 0.8, 0.9, and 1, where ${\delta }_{max}$ is the punch displacement at the onset of fracture, respectively, for Calib1 and Calib2. Included are two experiments, to establish the repeatability of these results (for clarity, only the final strains are shown for the 2nd experiment). For both parameters, the thickness strain variation had the greatest thinning at 45°, which indeed lead to the first rupture in the experiment; furthermore, the average thickness strain at five punch displacement levels was reasonably predicted. Nevertheless, the local accuracy of the thickness strain prediction was significantly improved in Calib2 compared to Calib1. This evidently indicates that the plastic anisotropy near the plane-strain condition plays an important role in the prediction of the thickness strain variation in the HE simulation. This will be further discussed in the following section.

4. Discussion

The influence of the accurate description of plastic anisotropy, especially for the plane-strain condition, in the HE simulation is assessed based on the radial (Figure 8a) and thickness (Figure 8b) strain distributions at the onset of fracture. In the radial strain distribution (Figure 8a), the contours for the plane-strain condition, identified in the figure at the boundary of the purple and blue colors, are observed to exist at different locations for Calib1 and Calib2: in Calib1, this location is much closer to the hole where the material deforms under uniaxial tension, while in Calib2 it is closer to the punch radius.

Considering the poor prediction of Calib1 (see Figure 3 and Figure 7), the proximity of the plane-strain condition to the hole indicates that it can interfere with the uniaxial tension deformation at the edge, much more than in the prediction of Calib2. In other words, an accurate description of the material behavior in biaxial stress states, especially near plane-strain, is important for the accurate prediction of the strain variation around the hole which can lead to the fracture. The thickness strain distribution (Figure 8b) can be understood in the same vein, with the thickness strain variation shown in Figure 4 at ${\delta }_{max}$: Calib1 predicts a weak variation of the thickness strain along the circumference and, as a result, it doesn’t capture the strain concentration near 45°, while Calib2 predicts the experiment well.

Beyond the strain state, the stress states at the onset of fracture are analyzed in the majority of the flat region of the deforming specimen. Data is extracted from 1586 nodes between the hole edge (r = 17.5 mm) and near the punch contact (r = 33.2 mm) as seen in Figure 9, in which the same color indicates the same radius from the center in the undeformed coordinate. Both Calib1 and Calib2 are plotted in the plane-stress space (${\sigma }_{xx},{\sigma }_{yy},{\sigma }_{xy}$) and its projection on the plane where ${\sigma }_{xy}$ = 0 in Figure 10 and Figure 11, respectively. In both cases, the stress states gradually change from uniaxial tension in every orientation between the RD and TD to equibiaxial tension as the radial distance increases from the hole edge to near the punch radius. It should be noted that the stress states beyond uniaxial tension, e.g., simple shear, do not exist in HE, which indicates that the anisotropy in simple shear does not play a critical role in the HE simulation. In this regard, the stress states between uniaxial tension and plane-strain tension in different orientations play the most important role in the parameter calibration in HE analysis. Similar to the observation in Figure 8a, the poorly described plane-strain conditions in Calib1 result in the plane-strain condition being closely located to uniaxial tension (near the yellow dots in Figure 9), in comparison to Calib2 (near the orange dots).

In closing, it should be noted that the HE experiment includes many stress states in the 1st quadrant of the plane stress yield locus. As such, it can be used as a way to inversely identify the constitutive behavior of a given sheet. Exploring this idea is left as future work by the authors.

5. Conclusions

The role of plane-strain condition on thickness strain variation in HE is investigated through numerical simulation using two different material parameter sets of the anisotropic yield function Yld2000-2d. Both parameter sets can capture the anisotropy in uniaxial tension, which is the major stress state along the hole edge. The structural response (i.e., punch force-displacement) is also well-predicted by both sets. However, only Yld2k-Calib2, which is calibrated including the plane-strain tension condition, shows a good agreement with the experiment in the thinning predictions. The comparison emphasizes the importance of a more holistic representation of plastic anisotropy, e.g., by incorporating the plane-strain tension condition, for accurate predictions of hole-expansion simulations.