Calibration of Distortional Plasticity Framework and Application to U-draw Bending Simulations

Abstract

A new version of a distortional plasticity framework, the so-called homogeneous anisotropic hardening (HAH), was investigated regarding its calibration and its application to the numerical analysis of the U-draw bend test. First, the mechanical properties of a dual-phase steel sheet sample, DP780, were characterized using uniaxial tension, bulge and elastic loading-unloading tests to provide the data for the calibration of a conventional constitutive description assuming isotropic hardening. Then, tension-compression with different load reversal numbers and two-step tension tests were performed to produce strain path changes that lead to anisotropic hardening effects. The coefficients of this new model, called HAH₂₀, were determined with an optimization procedure. Several sets of coefficients were identified depending on the number of reversals considered in tension-compression. U-draw bend test simulations were also carried out to validate the HAH₂₀ model and its implementation in a finite element code, to compare the results with an older version of the HAH model family, and to assess the influence of the input data used in the coefficient calibration. In addition, the influence of the cross-loading effect on the U-draw bending predictions was examined. The main conclusion of this work is that the HAH₂₀ coefficients calibrated for the description of reverse loading depend on the number of reversals during the tension-compression tests. The results of the U-draw bending simulations indicate that this dependence may lead to significant differences in the predicted amounts of springback.

1. Introduction

Sheet metal forming has been usually regarded as an important research field in the material and manufacturing industries, in particular, automobile. This process is widely used to manufacture automotive parts because it allows the production of complex shapes with good reproducibility in a time-efficient way. During the design stage of an actual sheet metal forming process, finite element (FE) simulations have been a crucial factor reducing a considerable amount of the optimization cost and time.¹⁾ In the finite element approach, material deformation models, among other key technical and scientific innovations, have been essential factors to improve the accuracy of numerical simulations.²⁾ However, there is still a margin in which further advances in material models may significantly boost the accuracy and the speed of the simulations.

Motivated by environmental concerns, there has been an increased demand for advanced high strength steels (AHSS) or lightweight metals from the automobile manufacturers. AHSS exhibit much higher strengths and sufficient ductility compared with more conventional materials due to their complex microstructures.^3,4) However, their uncommon elasto-plastic behavior requires advanced material models for application to numerical forming process simulations, in particular, when the material is subjected to non-proportional loading.

The influence of non-linear strain paths on the plasticity of metals has been extensively investigated^5,6,7,8) because a significantly different hardening behavior occurs in comparison with the commonly used assumption in forming simulations, i.e., isotropic hardening. One of the most severe non-proportional loading, a single load reversal, leads to the Bauchinger effect, which corresponds to a much lower reloading yield stress followed by a transient stage of hardening with high rate. In addition, after a reversal, the flow stress remains usually lower compared with the monotonic loading flow stress, an effect called permanent softening. These phenomena are referred to as reverse loading effects.⁹⁾ Besides, a single strain path change involving two orthogonal stress deviators, which is called cross-loading, produces hardening fluctuations as well. Cross-loading occurs, for instance, in a sequence of two uniaxial tension tests in directions that are at about 60° from each other. This so-called two-step tension test may result in a lower reloading yield stress, called cross-loading contraction, or/and an overshooting of the monotonic flow curve, called latent hardening. These effects are referred to as cross-loading effects.¹⁰⁾ In fact, any single strain path change leads to effects that are intermediate between those of cross-loading and either monotonic or reverse loading. Since a loading history can be considered to be a succession of incremental single strain paths, the concepts of cross-loading and reverse loading effects are essential to investigate.

The different hardening behaviors observed in non-proportional loading have usually been modeled through kinematic hardening approaches, in which the yield surface translates in stress space.^11,12,13,14) These constitutive models have improved the accuracy of numerical simulation results. For instance, the finite element (FE) simulation of springback in sheet metal forming have produced better results compared with the conventional isotropic hardening assumption.¹⁵⁾ These constitutive models have demonstrated that considering non-proportional loading effect is important in sheet metal forming simulations, especially for springback predictions. In addition, another concept of hardening based on the distortion of the yield surface has been used mostly in combination with kinematic hardening.^{16,17,18,19,20)} These added features allow the prediction of complex hardening behavior during non-proportional loading. Meanwhile, the so-called homogeneous anisotropic hardening (HAH) family of models was proposed as an alternative description of reverse loading and cross-loading effects based on a pure distortional plasticity concept, that is, without a back-stress.^21,22,23,24) Several studies have indicated that the HAH model was effective in describing the complex hardening fluctuations of AHSS during non-proportional loading.^10,25,26) However, because of a number of issues in the description of strain path changes near cross-loading,^26,27) a new version of the HAH model family, called HAH₂₀, was proposed recently to eliminate these concerns.²⁴⁾ Moreover, the calibration of such rather complex constitutive model is usually an obstacle to a more wide-spread application, in particular in industry.

In this study, the HAH₂₀ model was investigated regarding the calibration method and its effectiveness, as well as its numerical performance in a well-known forming benchmark, namely, the U-draw bend test. In particular, the sensitivity of the coefficients to the experimental input data and the suitability of this model to predict springback in U-draw bending were of interest. First, the elastic chord modulus and its degradation with respect to the accumulated plastic strain were assessed from uniaxial tension loading-unloading cycles using the approach proposed by Yoshida et al.⁹⁾ Then, the conventional isotropic hardening and anisotropic yield function, which are required in the HAH framework, were calibrated using standard monotonic uniaxial tension and bulge test results. Finally, the coefficients of a previous (HAH₁₁)²¹⁾ and new (HAH₂₀)²⁴⁾ versions of the model were identified. Three calibrations were carried out for HAH₂₀ using various number of reversals in the tension-compression tests, leading to different loading histories and final effective strains. The different calibrated coefficient sets were then used for the simulations of the U-draw bend test, and predicted springback profiles were compared with each other. However, no effort was spent to provide the best possible agreement between experiments and predictions because springback depends on many parameters, in particular, the unloading elastic modulus, which was not investigated in detail in this work. The present work is a valuable step toward the establishment of a robust calibration method for the new HAH₂₀ model and to assess its performance in forming simulations for AHSS sheets.

2. Experiments

Mechanical tests were conducted on a 1.4 mm thick DP780 dual-phase steel sheet sample supplied by POSCO. All the tests conducted in this work were at least duplicated and often triplicated. In all the cases, an excellent reproducibility was observed. Most of the experiments were carried out to generate data for the calibration of the constitutive model. They are categorized below depending whether they were conducted under proportional or non-proportional loading. An additional experiment, a U-draw bend test, was conducted to assess the springback after removal of the load.

2.1. Proportional Loading Tests

2.1.1. Uniaxial Tension and Bulge Tests

Uniaxial tension and bulge tests were conducted to characterize the mechanical properties in proportional loading conditions. A 500 kN MTS tensile testing machine with two mechanical extensometers was employed for the uniaxial tension tests. The specimens were machined following the ASTM E8 standard in three directions with respect to the rolling direction, namely, 0°or RD, 45° or DD and 90° or TD to characterize material anisotropy. The tests were conducted at a strain rate of 0.002/s and the true stress-strain curves and r-values (width-to-thickness strain ratio) were measured within the uniform deformation range. An Erichsen bulge tester, Model 161, was employed for the characterization of balanced biaxial tension, namely, the true membrane stress-true thickness strain curve and the biaxial r-value, i.e. the TD-to-RD strain ratio in balanced biaxial tension. The test was analyzed based on digital image correlations (DIC) data to determine the strains and the curvature at the pole of the bulge specimen. An advantage of the bulge test is that it allows the determination of the hardening curve in a much wider strain range compared with the uniaxial tension test. Figure 1 shows the stress-strain curves measured in RD uniaxial tension and bulge tests.

Fig. 1.

Effective stress-strain curves in RD uniaxial tension and bulge tests for DP780 steel sheet sample.

Normalized flow stresses are needed to estimate the material anisotropy and calculate the yield function coefficients. A 0.2% yield point method is usually not recommended for this purpose because material anisotropy fluctuates in the early stage of plastic deformation. The anisotropy of the yield stress is only a snapshot of the material behavior at zero plastic strain but it is not representative of the behavior at larger strains. Therefore, the flow stresses at a plastic work per unit volume of 30 MPa, which corresponds to a strain of 0.04 in the uniaxial tension test, were selected as representative for this DP780 because it was observed that plastic anisotropy tends to remain unchanged after this amount. The normalized flow stresses were calculated based on the same plastic work per unit volume for all the hardening curves including the balanced biaxial flow stress. The r-values in uniaxial tension were calculated by linear approximation of the width strain as a function of the longitudinal strain in the entire uniform deformation range, assuming plastic incompressibility for the determination of the thickness strain. The biaxial r-value was calculated as the TD-to-RD strain ratio using digital image correlation data. Table 1 lists the measured mechanical properties including the normalized flow stresses and r-values of the DP780 steel sheet sample determined from the analyses of the uniaxial tension and bulge tests.

Table 1. Mechanical properties measured from uniaxial tension and bulge tests for DP780 steel sheet.

DP780	Elastic modulus (GPa)	0.2% yield stress (MPa)	Uniform elongation (true strain)	UTS (MPa)	Total elongation (true strain)	Normalized stress	r-value
RD uniaxial tension	198.0	526.13	0.115	998.4	0.184	1.000 (σ1†)	0.776 (r1†)
DD uniaxial tension	–	519.09	0.122	997.0	0.177	0.991 (σ2†)	1.103 (r2†)
TD uniaxial tension	–	531.67	0.106	1012.1	0.159	1.022 (σ3†)	0.913 (r3†)
Bulge test	–	472.32	–	–	–	1.010 (σ4†)	0.887 (r4†)

–: not measured, ^†: notations for Eq. (17).

2.1.2. Elastic Unloading-reloading Test

A uniaxial tension test was conducted up to about 6% plastic strain using the same conditions described in Section 2.1.1 but with repeated elastic unloading-reloading cycles. 15 cycles were carried out at low plastic strains with 0.5% engineering strain intervals and 11 cycles at larger strain with 1% intervals. Figure 2(a) shows the true stress-strain curve with all the unloading-reloading cycles. The elastic chord modulus, as schematically defined in Fig. 2(b), was measured based on the extensometer and load cell information. A drop of the elastic chord modulus of about 20% was observed at a strain close to 6%.

Fig. 2.

(a) True stress-strain curve in RD uniaxial tension with elastic unloading and reloading cycles. (b) Chord modulus definition.

2.2. Non-proportional Loading Tests

2.2.1. Tension-compression Test

Tension-compression tests were conducted on a custom-built tension-compression machine to characterize the reverse loading behavior.²⁸⁾ A 7 kN normal force was applied on the specimen faces to prevent buckling. The grip displacement rate was set to 0.03 mm/s. A laser extensometer was used to measure the longitudinal strain. The other experimental details for this test can be found elsewhere.²⁶⁾ The tests were conducted with three different histories leading to various amounts of accumulated strains. The first consisted in tension-compression-tension with two reversals at about 3 and –3% true strain. The second was similar but with reversals at 6 and –6% strain. Figure 3(a) represents the stress-strain response of DP780 for these two loading histories. Finally, the third case consisted in a series of tension and compression segments with the first reversal at about 2% strain and each subsequent reversal after a strain interval of 2% higher than the previous segment, leading to the stress-strain curve shown in Fig. 3(b). This last history was considered because it leads to a large accumulated strain in comparison with a more conventional test such as uniaxial tension or tension-compression with a single or two reversals.

Fig. 3.

True stress-strain curves in RD tension-compression tests for DP780 steel sheet. (a) Tension-compression-tension with reversals at ±3% and ±6% strains; (b) Tension-compression cycles with six reversals and increasing strain intervals.

2.2.2. Two-step Tension Test

A two-step uniaxial tension test was conducted with a pre-strain of 6% in the TD followed by a subsequent deformation in the DD (Fig. 4(a)), which is an optimum strain path change for the calibration of the cross-loading coefficients.²⁶⁾ First, the pre-strain TD uniaxial tension was performed on a large-scale specimen. Then, a standard tension specimen, ASTM E8, was machined in the diagonal direction (DD) from the uniformly deformed region of the pre-strained specimen as measured by the DIC technique. The DD uniaxial tension test was conducted according to the experimental procedure explained in Section 2.2.1. The engineering stress-strain curve was converted into true stress-effective strain using the plastic work equivalence. Figure 4(b) shows the stress-strain curves measured in the DD direction for the monotonic and two-step tension tests.

Fig. 4.

(a) Schematic representation of a two-step tension test consisting of TD uniaxial tension on large scale specimen and DD uniaxial tension on standard specimen. (b) True stress-effective strain curves in DD for DP780 steel sheet in monotonic loading for as-received material, and after 6% pre-strain.

2.3. U-draw Bend Test

A U-draw bend test was carried out because, due to a large springback, this test is particularly interesting for the validation of sheet metal forming simulations. The U-draw bend tester is composed of a U-shaped die with a 5 mm radius, a punch with the complementary shape, and a blank-holder, the dimensions of which are shown in Fig. 5(a). A 350 by 45 mm rectangular specimen with the length parallel to the rolling direction was maintained on the die by a holder force of 40 kN. A 2.5 mm clearance between the walls of the die and the sides of the punch was established for the 1.4 mm thick DP780 steel sheet sample. The punch moved at a velocity of 1.2 mm/s for a total displacement of 70 mm. After the punch was removed, an elastic distortion occurred due to the redistribution of the residual stress field in the specimen. The resulting springback profile was measured using a laser-measuring device, Keyence LK-H150, and the springback parameters, namely, punch and die corner angles θ₁, and θ₂, and the wall radius of curvature ρ defined in Fig. 5(b), were calculated based on this profile.

Fig. 5.

(a) Schematic and dimensions of U-draw bend test; (b) Springback parameters θ₁, θ₂ and ρ in U-draw bending.

3. Modeling

3.1. Chord Modulus Degradation

An analysis of Fig. 2 indicated that the unloading-reloading chord modulus decreases as the plastic deformation accumulates. The decrease is initially very steep but tends to saturate as the longitudinal plastic strain ε reaches a certain value. Therefore, the empirical chord modulus degradation model proposed by Yoshida et al.⁹⁾ was employed in this work to describe the elastic loss of stiffness as the plastic strain increases, i.e.   

$$E=E_0-\left(E_0-E_s\right)\left[1-\text{exp}\left(-\xi \epsilon \right)\right]$$(1)

In the above equation, E is the current chord modulus with initial and saturated values of E₀ and E_s, respectively, and ξ a coefficient that characterizes the rate of decrease as a function of the uniaxial plastic strain ε.

3.2. Isotropic Hardening

Isotropic hardening is not only one of the assumptions commonly made in numerical simulations of forming processes but it is a basic building block of the HAH₂₀ distortional plasticity model. With the isotropic hardening assumption, the yield condition of the material can be expressed as   

$$\overline{\varphi }\left(\mathbf{s}\right)={\sigma }_r(\overline{\epsilon })$$(2)

where the anisotropic effective stress φ is a function of the stress deviator s and the reference flow curve σ_r, a function of the plastic work-based effective strain ε. Among the different approximations that were assessed for the isotropic hardening representation, a combination of the conventional Swift and Voce equations with an additional linear term was selected in this work because it was the most performant, i.e.   

$${\sigma }_r\left(\overline{\epsilon }\right)=K{\left({\epsilon }_0+\overline{\epsilon }\right)}^n+\beta \overline{\epsilon }+{\sigma }_y+{\sigma }_b\left(1-exp\left(-\eta \overline{\epsilon }\right)\right).$$(3)

where K, ε₀, n, β, σ_y, σ_b, and η are the hardening coefficients.

Although the effective stress φ can be associated with any anisotropic yield function, it was defined in this study using the plane stress Yld2000-2d anisotropic yield function²⁹⁾ expressed as   

$$\overline{\varphi }\left(\mathbf{s}\right)={\left\{\frac{{\left|X_1^{\prime }-X_2^{\prime }\right|}^a+{\left|2X_2^{″}+X_1^{″}\right|}^a+{\left|2X_1^{″}+X_2^{″}\right|}^a}{2}\right\}}^{1/a}$$(4)

The Yld2000-2d is based on the isotropic yield function proposed by Hershey.³⁰⁾ The principal values $X_p^{\prime }$ and $X_q^{″}$ of two linear transformations of the stress deviator, i.e., X′=C′:s and X″=C″:s with 4th order tensors C′ and C″, were substituted for those of the stresses in the original Hershey’s isotropic formulation. For plane stress, the representations of the double dot products of s by the two tensors C′ and C″, that contain the anisotropy coefficients, reduce to   

$$\left(\begin{array}{c}{X}_{11}^{\prime }\\ X_{22}^{\prime }\\ X_{12}^{\prime }\end{array}\right)=\left[\begin{array}{ccc}{C}_{11}^{\prime }& C_{12}^{\prime }& 0\\ C_{21}^{\prime }& C_{22}^{\prime }& 0\\ 0& 0& C_{66}^{\prime }\end{array}\right]\left(\begin{array}{c}{s}_{11}\\ s_{22}\\ s_{12}\end{array}\right)$$(5)

  

$$\left(\begin{array}{c}{X}_{11}^{″}\\ X_{22}^{″}\\ X_{12}^{″}\end{array}\right)=\left[\begin{array}{ccc}{C}_{11}^{″}& C_{12}^{″}& 0\\ C_{21}^{″}& C_{22}^{″}& 0\\ 0& 0& C_{66}^{″}\end{array}\right]\left(\begin{array}{c}{s}_{11}\\ s_{22}\\ s_{12}\end{array}\right).$$(6)

Since it was shown that only three coefficients $C_{pq}^{\prime }$ are independent in the yield function, Eq. (5), $C_{12}^{\prime }$ and $C_{21}^{\prime }$ were both set to zero. Therefore, the model contains only eight coefficients, which can be expressed with another set α_k with k = 1···8. The advantage of the coefficients α_k is that they all reduce to 1 for an isotropic material. The relationships between these two sets can be expressed as follows   

$$\left(\begin{array}{c}{C}_{11}^{\prime }\\ C_{22}^{\prime }\\ C_{66}^{\prime }\end{array}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left(\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ {\alpha }_7\end{array}\right)$$(7)

  

$$\left(\begin{array}{c}{C}_{11}^{″}\\ C_{12}^{″}\\ C_{21}^{″}\\ C_{22}^{″}\\ C_{66}^{″}\end{array}\right)=\frac{1}{3}\left[\begin{array}{rrrrr}\hfill -1& \hfill 0& \hfill 4& \hfill 0& \hfill 0\\ \hfill 0& \hfill -2& \hfill 0& \hfill 2& \hfill 0\\ \hfill 2& \hfill 0& \hfill -2& \hfill 0& \hfill 0\\ \hfill 0& \hfill 4& \hfill 0& \hfill -1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 3\end{array}\right]\left(\begin{array}{c}{\alpha }_3\\ {\alpha }_4\\ {\alpha }_5\\ {\alpha }_6\\ {\alpha }_8\end{array}\right).$$(8)

For more details about this model, the reader is referred to the original work.²⁹⁾

3.3. Anisotropic Hardening

Barlat et al.²⁴⁾ proposed a new version of the homogeneous anisotropic hardening model called HAH₂₀ to describe anisotropic hardening induced by non-proportional loading. In the present article, a short explanation is provided to introduce the yield condition of the HAH₂₀ model. The general yield condition is pressure-dependent because AHSS usually exhibit a large strength-differential (S-D) effect, i.e. a difference between uniaxial tension and compression flow curves. However, in the present work, this effect was neglected because the strength of DP780 is not too high. With this assumption, the yield condition reduces to   

$$\overline{\sigma }\left(\mathbf{s},v_k\right)={\left\{\overline{\xi }{\left(\mathbf{s}\right)}^q+{\overline{\varphi }}_v\left(\mathbf{s},v_k\right)\right\}}^{1/q}={\sigma }_R\left(\overline{\epsilon }\right)$$(9)

where σ is the effective stress of the HAH₂₀ model, which depends on the deviatoric stress tensor s and a number of state variables v_k. ξ is a modified version of the yield function φ in Eq. (4) that account for cross-loading phenomena. φ_v is a function allowing the Bauschinger effect, q a constant exponent and σ_R = g_Pσ_r the monotonic reference flow stress modified by permanent softening through one of the state variable g_P. There are two possible choices for φ_v but in this study, the following was selected   

$${\overline{\varphi }}_v={\overline{\varphi }}_v\left(\mathbf{s},f_{-},f_{+},\hat{\mathbf{h}}\right)=f_{-}^q{\left|\hat{\mathbf{h}}:\mathbf{s}-\left|\hat{\mathbf{h}}:\mathbf{s}\right|\right|}^q+f_{+}^q{\left|\hat{\mathbf{h}}:\mathbf{s}+\left|\hat{\mathbf{h}}:\mathbf{s}\right|\right|}^q$$(10)

where f₋, and f₊ are two state variables that control the reverse loading distortion. f₋ and f₊ can be expressed as functions of two other state variables g₋ and g₊, namely   

$$f_{\omega }=\sqrt{\frac{3}{8}}{\left\{\frac{1}{g_{\omega }^q}-1\right\}}^{\frac{1}{q}}$$(11)

with index ω (≡ + or −). The variables g₋ and g₊ have clearer physical meanings associated with the reloading yield stress after reversal as explained for g₋ in the Appendix.

$\hat{\mathbf{h}}$ is a tensorial state variable, the so-called microstructure deviator, which is normalized, i.e.   

$$\hat{\mathbf{h}}:\hat{\mathbf{h}}=1$$(12)

$\hat{\mathbf{h}}$ serves as a memory for the loading history and controls the location of the distortion. The normalized stress tensor $\hat{\mathbf{s}}$ leading to the first plastic strain increment is taken as the initial deviator $\hat{\mathbf{h}}$. When the load changes abruptly, the double dot product between the microstructure deviator, which does not change instantly, and the new normalized stress deviator, $cos\chi =\hat{\mathbf{h}}:\hat{\mathbf{s}}$, characterizes the severity of the load change. cosχ remains equal to 1 during proportional loading, and changes to 0 for pure cross-loading and −1 for reverse loading. When cosχ>0, the material deforms in a combined mode, i.e., monotonic and pure cross-loading while, when cosχ<0, the combined mode is reversal and pure cross-loading. Therefore, in this study, the former mixed mode is called cross-loading and the latter reverse loading. After an abrupt load change, $\hat{\mathbf{h}}$ rotates towards the new stress state at a rate given by an evolution equation. Since, this evolution sometimes leads to unwanted fluctuations of the flow stress when the strain path changes are close to cross-loading³¹⁾ another deviator, ${\hat{\mathbf{h}}}^{\prime }$, was introduced in the formulation to delay this rotation in the HAH₂₀ model.

For a given microstructure deviator $\hat{\mathbf{h}}$, any stress state represented by s can be projected as the sum of two tensors, ${\mathbf{s}}_C=(\mathbf{s}:\hat{\mathbf{h}})\hat{\mathbf{h}}$ and ${\mathbf{s}}_O=\mathbf{s}-(\mathbf{s}:\hat{\mathbf{h}})\hat{\mathbf{h}}$. The first one, collinear to $\hat{\mathbf{h}}$, corresponds to the amount of monotonic or reverse loading, depending on the sign cosχ. The second, orthogonal to $\hat{\mathbf{h}}$, corresponds to the amount of cross-loading. The quantity ξ(s) in Eq. (9) is defined as follows   

$$\overline{\xi }\left(\mathbf{s}\right)={\left[\overline{\varphi }{\left({\mathbf{s}}_L\right)}^p+\overline{\varphi }{\left({\mathbf{s}}_X\right)}^p\right]}^{1/p}$$(13)

with   

$${\mathbf{s}}_L=\frac{1}{{\xi }_L\left(g_L-1\right)+1}{\mathbf{s}}_C+\frac{1}{g_L}{\mathbf{s}}_O$$(14)

and   

$${\mathbf{s}}_X={\xi }_C\frac{1-g_C}{g_L}{\mathbf{s}}_O$$(15)

Equation (13) is a modified version of the anisotropic yield function φ(s) that includes cross-loading effects. g_L and g_C are the state variables that control latent hardening and cross-loading contraction, respectively, and ξ_L is a constant coefficient. If proportional loading is used in this formulation, the material response is the same as that of the original anisotropic yield function φ(s). Note that the yield condition of the HAH₁₁ model is similar to Eqs. (9), (10), (11) except that ξ(s) is only φ(s), which indicates that cross-loading effects are not accounted for in HAH₁₁. For illustration purpose, an example of yield surface distortion, during and just after strain path change, is given in the Appendix.

The full description of the HAH₂₀ model is detailed elsewhere.²⁴⁾ However, Table A1 containing all the state variables, as well as a summary of all the corresponding evolution equations, are provided in the Appendix of this article. Moreover, all the material coefficients are also listed in Table A2 in the Appendix. Some coefficients can be used with the recommended values that are provided in this table while other must be determined using an optimization procedure, which is the topic of Section 4.1. All the coefficients that are not mentioned in the following discussion are set to the recommended values in Table A2.

3.4. Finite Element Modeling

The commercial finite element software ABAQUS was used for the simulation of the U-draw bend test in this work. A sketch of the different components used in the simulations is presented in Fig. 6. The mesh, selected through a parametric study, consisted of 0.5 mm square shell elements (S4R) with 9 integration points through the thickness. The friction coefficient was set to 0.1, which is a well-accepted value for this type of simulations. The numerical forming process was analyzed using the explicit finite element method to avoid convergence issues due to the contact between the tools and the blank. The springback after punch removal was conducted using the implicit finite element method. User-defined material subroutines (VUMAT for ABAQUS explicit and UMAT for ABAQUS implicit) implemented as in-house codes with the constitutive models described in this article were employed for the simulation of the process.³²⁾

Fig. 6.

Illustration of assembly in finite element modeling of a U-draw bend test. (Online version in color.)

4. Results and Discussion

4.1. Constitutive Model Calibration

Four material models are combined to describe the constitutive behavior, namely, elastic chord modulus degradation, anisotropic yield function, reference stress-strain curve, and anisotropic hardening. In this study, all the coefficients were determined based on the minimization of an objective function Θ, the least-square root error, characterizing the difference between experiments and predictions   

$$\mathrm{\Theta }={\displaystyle \sum }_{i=1}^n{\left[y_i{\left(x_i\right)}_{experiment}-y_i{\left(x_i\right)}_{simulation}\right]}^2$$(16)

where y_i is the property of interest as a function of the variable x_i, and n the number of data involved in the calculations. Note that this approach is used for both isotropic and anisotropic hardening models. The predicted values are all calculated using a stand-alone code based on the distortional plasticity framework that reproduces the response of the constitutive model only, given a specific loading history. The Nelder-Mead simplex method is selected for the optimization scheme because it requires only functional values, which is a significant advantage in calibrating complex constitutive models. The calibration is made in a sequential manner with several stages, that is, 1) elastic chord modulus, 2) isotropic hardening flow curve and 3) anisotropic yield function. These three stages lead to the respective determination of 3, 7 and 8 coefficients with the constitutive description selected in this study but do not present any specific difficulty. In particular, even though the coefficients are determined assuming isotropic hardening, they are not modified thereafter.³³⁾ Only then, the anisotropic hardening coefficients are determined but, again, in two sequential stages because the reverse loading coefficients are independent from those for cross-loading.

4.1.1. Elasticity

For the elastic degradation, the three unknown coefficients (E₀, E_s, ξ) are determined using the chord modulus as the variable, i.e., y_i = E_i, as a function of the uniaxial plastic strain x_i = ε_i in Eq. (16). Table 2 lists the calibrated coefficients and Fig. 7 shows the elastic chord modulus as a function of the plastic strain. As mentioned earlier, the chord modulus decreases sharply at low plastic strain and tends to saturate at strains larger than 0.06. Figure 7 indicates that the chord modulus model proposed by Yoshida and Uemori⁹⁾ captures this degradation very well. Note that, as a generalization, it is assumed that the chord modulus is a function of the effective strain in all the subsequent calculations. Although some recent researches^34,35,36,37) reported that elasticity depends on the stress state, the above assumption is selected in this work because many articles^{31,33,38,39,40,41,42,43,44,45)} indicate that it leads to accurate prediction of springback.

Table 2. Yoshida-Uemori⁹⁾ chord modulus coefficients for DP780 steel sheet.

E0 [GPa]	Es [GPa]	ξ
198.00	155.54	107.37

Fig. 7.

Experimental and approximated elastic chord modulus vs. effective strain based on unloading-reloading cycles in the RD uniaxial tension for DP780 steel sheet.

4.1.2. Isotropic Hardening

Usually the stress-strain curve measured with the bulge test is obtained in a strain range, which is about twice as large as that of uniaxial tension. However, as a drawback, the flow stress determined within the first few percent strain is not very accurate because its calculation depends on the bulge specimen curvature, which is infinite initially. In order to overcome this difficulty, the reference stress-strain curve is calibrated using the data of the RD uniaxial tension and bulge tests simultaneously. This requires the conversion of the bulge stress-strain curve to RD uniaxial tension data, based on plastic anisotropy and isotropic hardening to combine the two tests. The membrane stress of the bulge test is normalized by the uniaxial flow stress and the equivalent strain is multiplied by this normalized value in order to preserve the plastic work equivalence. After this procedure, the RD uniaxial tension data defining the hardening behavior in a strain range [0–0.1] and the processed bulge test data for strains larger than 0.1 are combined in a single data set. The combined hardening curve is approximated by the Swift-modified Voce equation using the identification procedure described above with Eq. (16). The resulting coefficients are listed in Table 3. The experimental and approximated flow curves are shown in Fig. 8, which indicates that this method allows the determination of an accurate stress-strain curve from the strain at yield to about 0.5.

Table 3. Detailed coefficients for the Swift-modified Voce equation for DP780 steel sheet.

K [MPa]	ε0	n	β [MPa]	σy [MPa]	σb [MPa]	η
738.46	0.0008	0.25	15.99	349.39	215.83	59.73

Fig. 8.

Experimental stress-strain curve with Swift-modified Voce hardening approximation for DP780 steel sheet. Combined uniaxial tension test data for strain smaller than 0.1 and bulge test data corrected for plastic anisotropy for strain larger than 0.1.

Since the yield function is defined in terms of normalized flow stresses and r-values as described in the experimental section (2.1.1), a slightly different form of the objective function Θ is chosen in this case, namely,   

$$\mathrm{\Theta }=w_1{\displaystyle \sum }_{i=1}^{n_1}{\left(\frac{{\sigma }_{i-simulation}}{{\sigma }_{i-experiment}}-1\right)}^2+w_2{\displaystyle \sum }_{i=1}^{n_2}{\left(\frac{r_{i-simulation}}{r_{i-experiment}}-1\right)}^2$$(17)

where w₁ and w₂ are weight factors, and n₁ and n₂ the number of terms considered in each summation. The set w₁ = 1 and w₂ = 0.1 would be a reasonable choice because an error of 0.05 is significant for a normalized flow stress but negligible for a r-value. In fact, since the numbers of input data and yield function coefficients are the same in this work, n₁ + n₂ = 8, these weight factors w₁ and w₂ have no influence on the results.

In Eq. (17), the four normalized flow stresses σ_i and the four r-values r_i listed in Table 1 are the properties of interest needed to calibrate the eight coefficients α_k described in Section 3.2. These properties are functions of the specific stress states considered, uniaxial and biaxial tension. The Yld2000-2d yield function exponent a is assumed to be 6 as recommended for BCC type materials,⁴⁶⁾ although other values could possibly be used if additional test results, such as in Kuwabara et al.,⁴⁷⁾ were available. Table 4 summarizes the results of this calibration showing that, since the DP780 is not very anisotropic, the coefficients of the yield function are all close to one.

Table 4. Yld2000-2d yield function coefficients for DP780 steel sheet.

α1	α2	α3	α4	α5	α6	α7	α8	a
0.97241	0.96969	0.98421	0.98559	1.00994	0.96496	1.01517	1.01556	6

4.1.3. Anisotropic Hardening

The three state variables g₋, g₊ and g_p play a role in the description of reverse loading, i.e., Bauchinger-related effect and permanent softening. Moreover, the two-state variables g_C and g_L control cross-loading contraction and latent hardening effects, respectively. The evolution equations of these state variables are described in detail in the original article²⁴⁾ but, for convenience, they are also summarized in the Appendix with all the coefficients listed in Table A.2. This table indicates that 22 coefficients are available but about half of them can be set to recommended values. Therefore, in the HAH₂₀ model, the calibration of only 10 coefficients is necessary.

The three variables g₋, g₊ and g_p require the calibration of five coefficients, namely, k₁, k₂, k₃, k₄ and k₅ using the reverse loading test results. The two variables g_C and g_L also require five coefficients C, k_C, $k_C^{\prime }$, L, and k_L using the cross-loading data. An advantage of the HAH approach is that all the coefficients do not have to be determined simultaneously.³³⁾ The five coefficients k₁, k₂, k₃, k₄ and k₅ are calibrated with the tension-compression data in this study and the other set, C, k_C, $k_C^{\prime }$, L, and k_L, using the two-step TD-DD tension data. As mentioned above, the influence of pressure is not considered in this study because the S-D effect is not very strong for DP780 but this aspect will be discussed later.

First, all the coefficients determined for elasticity and isotropic hardening, namely, elastic modulus degradation, reference stress-strain curve and anisotropic yield function are listed in Tables 2, 3 and 4, and is referred to coefficient Set 1. This set of coefficients remains the same even when incorporated in the anisotropic hardening models. Then, the coefficients k₁, k₂, k₃, k₄ and k₅ of the HAH₁₁ model are calibrated from the tension-compression cycles with six reversals as shown by Fig. 3(b), leading to coefficient Set 2. The definitions of these coefficients are slightly different from those of the HAH₂₀ model but they also characterize the Bauschinger effect. The coefficients of Set 2 are listed in Table 5.

Table 5. Identified HAH coefficients and initial values & ranges for optimization.

Set 1	Isotropic hardening assumption
Set 2 (HAH11)	k1	k2	k3	k4	k5	k
	102.15	89.61	0.54	0.78	13.05	30
Set 3-1 (HAH20)	k1	k2	k3	k4	k5	C	kC	$k_C^{\prime }$	${\xi }_C^{\prime }$	L
	84.34	97.38	0.23	0.81	5.01	0.63	164.15	29.23	6	1
Set 3-2 (HAH20)	k1	k2	k3	k4	k5	C	kC	$k_C^{\prime }$	${\xi }_C^{\prime }$	L
	82.99	84.63	0.44	0.84	13.07	0.63	166.12	29.23	6	1
Set 3-3 (HAH20)	k1	k2	k3	k4	k5	C	kC	$k_C^{\prime }$	${\xi }_C^{\prime }$	L
	78.84	85.51	0.49	0.79	14.13	0.63	164.96	29.21	6	1
Set 4 (HAH20)	k1	k2	k3	k4	k5	C	kC	$k_C^{\prime }$	${\xi }_C^{\prime }$	L
	78.84	85.51	0.49	0.79	14.13	1	0	0	6	1
Initial [Ranges]	75 [50:200]	75 [50:100]	0.40 [0.1:0.8]	0.80 [0.7:1]	15 [5:20]	0.75 [0.5:1]	105 [10:200]	105 [10:200]

For the HAH₂₀ model, only the reverse loading coefficients are considered at this stage. Three sets are calibrated to assess the influence of the type of input curve in the identification process. One, called Set 3-1, is calibrated from data extracted from the tension-compression cycles in Fig. 3(a) with only one single reversal. Two tests with pre-strains of 3 and 6% are considered leading to approximately 9% and 18% accumulated strains, respectively. This type of input data with only one single reversal is quite commonly used as surveyed from previous researches.^{38,39,40,48,49,50,51,52)} Another set, called 3-2, is calibrated from the tension-compression-tension cycles with the same strain intervals as shown in Fig. 3(a), that is, two reversals leading to approximately 15% and 30% accumulated strains, respectively. The last Set 3-3 is calibrated, as for Set 2, using the cycles of Fig. 3(b) with six reversals leading to approximately 54% accumulated strains. All the reverse loading coefficients of Sets 3-1, 3-2 and 3-3 are listed in Table 5. Some differences between these three sets of k₁, k₂, k₃, k₄ and k₅ are noticeable, which means that the nature of the input data affect the calibration results. This aspect should be considered further as it might affect forming simulation results.

Figures 9(a) to 9(e) compares the approximated and experimental tension-compression stress-strain curves determined with the different coefficient sets. Figure 9(a) indicates that the predicted stress-strain curves calculated assuming isotropic hardening are not in very good agreement with the experimental curves. However, for anisotropic hardening and the different coefficients sets, the agreement is good in all cases, in particular, for HAH₁₁ (Set 2) in Fig. 9(b) and HAH₂₀ (Set 3-3) in Fig. 9(e), for which six reversals were carried out. Only before the 4th and 6th reversals, the strain hardening is slightly underpredicted. However, the HAH formulations (HAH₁₁ and HAH₂₀) are not likely the reason for this discrepancy because an examination of Fig. 9(a) indicates that the isotropic strain hardening is also underestimated before the 4th and 6th reversals. Since the main focus of this work pertains to the plastic behavior during strain path change, the elastic model was described with a simple chord modulus approach. Figure 9 indicates that this description of the elastic behavior is acceptable at first order.

Fig. 9.

Experimental and calculated true strain-stress curves in tension-compression (T-C) with different histories for DP780 sheet sample; (a) Cycles with six reversals and increasing strain intervals; Isotropic hardening and coefficient Set 1; (b) Cycles with six reversals and increasing strain intervals; HAH₁₁ model and coefficient Set 2; (c) T-C with one reversals and two pre-strains (3 and 6%); HAH₂₀ model and coefficient Set 3-1; (d) T-C-T with two reversals at ±3 and ±6%; HAH₂₀ model and coefficient Set 3-1; (e) Cycles with six reversals and increasing strain intervals; HAH₂₀ model and coefficient Set 3-3.

The calibration of the cross-loading coefficients was performed for the HAH₂₀ model to complete the three Sets 3-1, 3-2 and 3-3. In general, it is necessary to employ the reverse loading coefficients for the calibration of the cross-loading behavior because the second loading segment is usually not pure cross-loading. However, in many cross-loading scenarios, this influence is inexistent. This aspect is discussed further in the next paragraph. In the two-step tension test of the DP780 steel sheet in Fig. 4, only cross-loading contraction occurs, no latent hardening, because no flow stress overshooting is observed. Therefore, in the absence of latent hardening, L is set to 1, which renders k_L irrelevant. In addition, the coefficient ${\xi }_C^{\prime }$ is set to the value of 6, as recommend in the original HAH₂₀ article (see also Appendix). Thus, only three cross-loading coefficients are involved in this calibration, namely, C, k_C and $k_C^{\prime }$. These coefficients are optimized by minimizing Eq. (16) with the suitable flow stresses and strains as variables y_i and x_i. The complete coefficient Sets 3-1, 3-2 and 3-3 are listed in Table 5. Note that, in addition to the HAH₂₀ coefficients, the initial ranges and search regions needed for the simplex optimization method are listed in this table.

Figure 10 shows the experimental flow curves of the two-step tension test with those calculated with HAH₂₀ and the different calibrated coefficient sets (3-1, 3-2 and 3-3). These curves are all in excellent agreement with the experimental curve, irrespective of the coefficient set. Although some differences in reverse loading coefficients were observed, the cross-loading coefficients were almost identical. This indicates that the reverse loading coefficients did not affect the calibration of the cross-loading coefficients significantly. This is because the two-step tension test in Fig. 10 corresponds to cosχ = 0.25 which is positive. This means that the strain path change is a combination of monotonic and cross-loading and does not involve the reverse loading coefficients in HAH₂₀. This observation means that the calculation of the reverse loading and cross-loading coefficients with the procedure described in this article is, as expected, completely independent.

Fig. 10.

True stress-effective strain curves in DD uniaxial tension after 6% TD uniaxial tension pre-strain, measured and predicted with HAH₂₀ with three coefficient sets (3-1, 3-2 and 3-3), which lead to the same excellent prediction.

In Section 4.2, the springback predictions in the U-draw bend test are performed for three constitutive models, namely, isotropic hardening (Set 1), HAH₁₁ (Set 2) and HAH₂₀ (Set 3-3). In addition, the springback predictions are conducted with the same model (HAH₂₀) but with the three sets of coefficients (3-1, 3-2 and 3-3) corresponding to the three types of input curves for reverse loading. Finally, a new set of coefficients call Set 4 was generated to exclude cross-loading effects from Set 3-3 by imposing C = 1, which renders k_C and $k_C^{\prime }$ irrelevant. Since it is expected that only reverse loading occurs in U-draw bending, this trial is carried out to further validate the HAH₂₀ model and its implementation.

4.2 U-draw Bending Simulations

The U-draw bending simulations were first carried out to investigate the influence of the constitutive model, isotropic hardening, HAH₁₁ and HAH₂₀ with the corresponding coefficient Sets 1, 2 and 3-3. The experimental and predicted punch load-displacement curves are compared in Fig. 11 to verify the validity of the finite element model. This comparison shows that for all the predicted cases the load saturates to about the same level in good agreement with the experimental value. This indicates that the selection of the friction coefficient is reasonable although the sensitivity of this parameter was not investigated in detail. In addition, the load computed using the isotropic hardening model is overall slightly larger than the loads predicted with the distortional plasticity models. This was expected since the Bauschinger effect and permanent softening, captured by HAH₁₁ and HAH₂₀, tend to decrease the strength of the material. Therefore, these results indicate that the finite element model of the U-draw bend test is reliable.

Fig. 11.

Punch load-displacement curves in U-draw bending for DP780 steel sheet, measured and predicted with isotropic hardening (Set 1), HAH₁₁ (Set 2) and HAH₂₀ (Set 3-3).

Figure 12(a) shows the experimental and predicted U-shape profiles after springback and Fig. 12(b) the associated parameters P, namely, θ₁, θ₂ and ρ, defined in Fig. 5(b). The values of these parameters, as well as the relative errors calculated as   

$$\text{Relative}\mathrm{\hspace{0.17em} } \text{error}\left(\text{\%}\right)=100*\left|\frac{P_{experiment}-P_{simulation}}{P_{experiment}}\right|$$(18)

using the experimental and simulated value of P (θ₁, θ₂ or ρ), are also summarized in Table 6. None of the model results lead to an accurate prediction of the experimental springback profile. However, as mentioned earlier, the absolute springback profile is difficult to predict because it is controlled by the values of three parameters θ₁, θ₂ and ρ. Small errors in these parameters are amplified by the long section of the flange, which does not deform plastically during the process. Since many other variables (physical and numerical) are likely to affect the forming and springback predictions, in particular, an accurate description of the unloading elastic modulus, no attempt was made to reduce the difference between the predicted and experimental profiles. In fact, the elastic behavior is one of the most important parameters to consider for the prediction of springback.^37,42,43) A number of studies^{34,45,51,53,54,55,56)} indicate that non-linear elasticity is more suitable than the chord modulus and provides more accurate predictions of springback. However, this has not been investigated further in the present work.

Fig. 12.

(a) 2d springback profiles in U-draw bending for DP780 steel sheet, measured and predicted with isotropic hardening (Set 1), HAH₁₁ (Set 2) and HAH₂₀ (Set 3-3); (b) Corresponding experimental and predicted springback parameters.

Regarding the relative influence of the model, isotropic hardening (Set 1) leads to the highest springback, and HAH₁₁ (Set 2) and HAH₂₀ (Set 3-3) lead to similar lower springback. Therefore, as many researches already concluded,^{31,39,41,42,43,44)} the influence of the hardening effects occurring during strain path changes on the springback predictions of the U-draw bending is significant. In fact, the parameters θ₁, θ₂ and ρ that characterize the plastic deformation are almost identical for these two anisotropic hardening models. The only noticeable difference is the wall curvature, which is not large enough to result in a significant change in the springback profile. Since HAH₁₁ and HAH₂₀ lead to almost identical calculated stress-strain tension-compression cycles, they both produce about the same amount of springback.

Table 6. Experimental and predicted springback parameters.

DP780	Experiment	Set 1 (isotropic)	Set 2 (HAH11)	Set 3-1 (HAH20)	Set 3-2 (HAH20)	Set 3-3 (HAH20)	Set 4 (HAH20)
θ1 (°)	113.32	116.15	115.22	116.02	114.95	114.54	114.63
Relative error (%)	—	2.49	1.68	2.38	1.44	1.08	1.15
θ2 (°)	76.55	68.49	71.20	69.33	70.63	71.34	71.25
Relative error (%)	—	10.53	6.99	9.43	7.73	6.79	6.93
ρ (mm)	98.13	87.07	93.38	87.76	93.16	95.73	95.81
Relative error (%)	—	11.27	4.84	10.57	5.06	2.45	2.36

Then, the U-draw bending profiles and parameters θ₁, θ₂ and ρ (see also Table 6) predicted using HAH₂₀ with Sets 3-1, 3-2 and 3-3 were compared to assess the influence of the input data used to calibrate the model coefficients. Figures 13(a), 13(b) indicate that Set 3-2 and Set 3-3 lead to better springback predictions than Set 3-1. In fact, using Set 3-1, the predicted profile is rather close from the isotropic hardening prediction even though the calculated tension-compression curve is in excellent agreement with the experimental curve. This result points out that the performance of the constitutive model in FE application is strongly affected by the input used for the calibration. In this case, the input data with at least 30% accumulated strain and two reversals leads to a better prediction of springback. In previous researches,^39,40,50) the HAH coefficients were calibrated from only one reversal or with a small amount of accumulated strain. However, this work clearly shows that two or more reversals with a larger accumulated strain is likely to provide better results. In addition, examination of Table 5 indicates that the value of k₅, the coefficient controlling the rate at which permanent softening develops, has probably the largest influence on the amount of springback in this case.

Fig. 13.

(a) 2d springback profiles in U-draw bending for DP780 steel sheet, measured and predicted with HAH₂₀ with coefficients Sets 3-1, 3-2 and 3-3. (b) Corresponding experimental and predicted springback parameters.

In fact, the permanent softening can be explained by the amount of recovery that occurs during reversal by the annihilation of a certain amount of dislocations.⁵⁷⁾ This phenomenon is, in fact, very difficult to model properly because for a given accumulated strain, the amount of recovery might depend on the number of reversals. Nevertheless, with only one reversal and a small amount of accumulated strain, an accurate assessment of permanent softening is likely not easy to achieve because it first requires a reliable reference flow stress in the absence of reversal recovery. This is the reason why the selected monotonic stress-strain curve was calibrated up to an effective strain larger than 0.5 in this work. However, this require the use of a bulge test for calibration at large strains.

In addition, the strength-differential (S-D) effect, which appears to be an important feature of high strength materials, reflects a significant role of the hydrostatic stress on plasticity.⁵⁸⁾ The S-D effect leads to a larger flow stress in compression while the permanent softening leads to the opposite effect. Therefore, with only one reversal, the S-D and permanent softening effects cannot be clearly separated. This is possibly the reason why, as noted above, the value of k₅ is different for the coefficient Set 3-1 compared with Sets 3-2 and 3-3. It takes at least two reversals to clearly deconvolute the S-D and permanent softening effects. This explains why, although not used in this work, the pressure effect was considered in the HAH₂₀ formulation. In fact, this might be one of the factors leading to the gap between the predicted and experimental amounts of springback in the present study. This issue is left up to a future investigation.

Finally, Sets 3-3 and 4 are compared to investigate the influence of the cross-loading effect in the U-draw bend test. All the coefficients of Set 4 are identical to those of Set 3-3 except for C which is set to 1 thus, disabling cross-loading contraction and making k_C and $k_C^{\prime }$ irrelevant. These two sets of coefficients, 3-3 and 4, were designed to extract the separate impacts of reverse loading and cross loading effects on springback. Figure 14 shows the springback profiles and the corresponding parameters θ₁, θ₂ and ρ (see also Table 6) for the experiment and the two simulations. Not only the springback profiles are the same but the three parameters are identical as well.

Fig. 14.

(a) 2d springback profiles in U-draw bending for DP780 steel sheet, measured and predicted with HAH₂₀ with coefficients Sets 3-3 and 4. (b) Corresponding experimental and predicted springback parameters.

These results indicate that the cross-loading effect has no influence on springback in the U-draw bend test. This is likely because, during U-draw bending, the dominant path change is pure reverse loading when a material is flowing from a flat area to a corner radius and vice-versa. Therefore, the small difference between HAH₁₁ and HAH₂₀ predictions observed in Fig. 11 are exclusively due to the change of the reverse loading formulation in the HAH₂₀ model. Consequently, the U-draw bend test is not appropriate to assess the improvement of HAH₂₀ compared with HAH₁₁ because the difference induced from the reverse loading is not distinguishable and the improvement of cross-loading cannot be estimated. In fact, these results were expected because both HAH₁₁ and HAH₂₀ models are accurate for pure reverse loading in Section 4.1.2. These simulations tend to confirm that the implementations of these models in the FE code were successful.

5. Conclusions

In this work, the calibration of a new version of the homogeneous anisotropic hardening model, HAH₂₀, was conducted on a DP780 steel sheet sample and the procedure was critically analyzed. The determination of the coefficients corresponding to elastic modulus degradation, reference flow curve and yield function were performed in three stages as if hardening was isotropic. For distortional effects, the calibration was conducted in two stages, one for the coefficients characterizing the effects of reverse loading and the other for those of cross-loading. This constitutive model was applied to the finite element (FE) simulation of the U-draw bend test. For comparison purpose, the calibrations and the simulations were also conducted using an older version of the homogeneous anisotropic hardening model namely, HAH₁₁,²¹⁾ which does not include the cross-loading effects.

• In spite of the relative complexity of the new HAH₂₀ model and the large number of coefficients, the sequential calibration procedure in five stages presented in this work is well adapted and relatively easy to carry out.

• The nature of the input data for the calibration of the HAH₂₀ model in reverse loading, namely, reversal number and total accumulated strain, has an influence on the value of the coefficients which, in turn, affects the prediction of springback.

• In the U-draw bend test, the amounts of springback predicted using HAH₂₀ and HAH₁₁ is very close because the description of the reverse loading effects is similar in both models and the influence of the cross-loading effects in this test is irrelevant.

• Based on this work, a calibration with a minimum of two load reversals and 30% accumulated strains is recommended for springback predictions and, more generally, for applications to sheet metal forming simulations.

Acknowledgments

The authors gratefully acknowledge the generous financial support of POSCO.

References

• 1)   A.  Makinouchi: J. Mater. Process. Technol., 60 (1996), 19. https://doi.org/10.1016/0924-0136(96)02303-5
• 2)   D.  Banabic and  M.  Sester: Mater. Manuf. Process., 27 (2012), 273. https://doi.org/10.1080/10426914.2011.578005
• 3)   E.  Billur and  T.  Altan: Proc. New Developments in Sheet Metal Forming, MAT INFO Werkstoff-Informationsgesellschaft, Frankfurt am Main, (2012), 285.
• 4)   P.  Fu,  L.  Peng and  W.  Ding: Strateg. Study Chin. Acad. Eng., 20 (2018), No. 1, 84.
• 5)   F.  Barlat,  J. M.  Ferreira Duarte,  J. J.  Gracio,  A. B.  Lopes and  E. F.  Rauch: Int. J. Plast., 19 (2003), 1215. https://doi.org/10.1016/S0749-6419(02)00020-7
• 6)   L.  Sun and  R. H.  Wagoner: Int. J. Plast., 45 (2013), 174. https://doi.org/10.1016/j.ijplas.2013.01.018
• 7)   T.  Mánik,  B.  Holmedal and  O. S.  Hopperstad: Int. J. Plast., 69 (2015), 1. https://doi.org/10.1016/j.ijplas.2015.01.004
• 8)   J.  Liao,  J. A.  Sousa,  A. B.  Lopes,  X.  Xue,  F.  Barlat and  A. B.  Pereira: Int. J. Plast., 93 (2017), 64. https://doi.org/10.1016/j.ijplas.2016.03.010
• 9)   F.  Yoshida,  T.  Uemori and  K.  Fujiwara: Int. J. Plast., 18 (2002), 633. https://doi.org/10.1016/S0749-6419(01)00049-3
• 10)   J.  Ha,  M. G.  Lee and  F.  Barlat: Mech. Mater., 64 (2013), 11. https://doi.org/10.1016/j.mechmat.2013.04.004
• 11)   P. J.  Armstrong and  C. O.  Frederick: A Mathematical Representation of the Multiaxial Bauschinger Effect, Report RD/B/N731, Central Electricity Generating Board, Berkeley, UK, (1966), 1.
• 12)   J.-L.  Chaboche: Int. J. Plast., 2 (1986), 149. https://doi.org/10.1016/0749-6419(86)90010-0
• 13)   F.  Yoshida and  T.  Uemori: Int. J. Plast., 18 (2002), 661. https://doi.org/10.1016/S0749-6419(01)00050-X
• 14)   R. K.  Verma,  K.  Chung and  T.  Kuwabara: ISIJ Int., 51 (2011), 482. https://doi.org/10.2355/isijinternational.51.482
• 15)   R. H.  Wagoner,  H.  Lim and  M. G.  Lee: Int. J. Plast., 45 (2013), 3. https://doi.org/10.1016/j.ijplas.2012.08.006
• 16)   U.  Kowalsky,  H.  Ahrens and  D.  Dinkler: Comput. Mater. Sci., 16 (1999), 81. https://doi.org/10.1016/S0927-0256(99)00049-X
• 17)   M.  Yeganeh: Int. J. Plast., 23 (2007), 2029. https://doi.org/10.1016/J.IJPLAS.2007.02.004
• 18)   M. P.  Pietryga,  I. N.  Vladimirov and  S.  Reese: Mech. Mater., 44 (2012), 163. https://doi.org/10.1016/J.MECHMAT.2011.07.014
• 19)   H. P.  Feigenbaum,  J.  Dugdale,  Y. F.  Dafalias,  K. I.  Kourousis and  J.  Plesek: Int. J. Solids Struct., 49 (2012), 3063. https://doi.org/10.1016/J.IJSOLSTR.2012.06.006
• 20)   J.  Qin,  B.  Holmedal and  O. S.  Hopperstad: Int. J. Plast., 101 (2018), 156. https://doi.org/10.1016/J.IJPLAS.2017.10.013
• 21)   F.  Barlat,  J. J.  Gracio,  M. G.  Lee,  E. F.  Rauch and  G.  Vincze: Int. J. Plast., 27 (2011), 1309. https://doi.org/10.1016/j.ijplas.2011.03.003
• 22)   F.  Barlat,  J.  Ha,  J. J.  Grácio,  M. G.  Lee,  E. F.  Rauch and  G.  Vincze: Int. J. Plast., 46 (2013), 130. https://doi.org/10.1016/j.ijplas.2012.07.002
• 23)   F.  Barlat,  G.  Vincze,  J. J.  Grácio,  M. G.  Lee,  E. F.  Rauch and  C. N.  Tomé: Int. J. Plast., 58 (2014), 201. https://doi.org/10.1016/j.ijplas.2013.11.002
• 24)   F.  Barlat,  S.-Y.  Yoon,  S.-Y.  Lee,  M.-S.  Wi and  J.-H.  Kim: Int. J. Solids Struct., 202 (2020), 947. https://doi.org/10.1016/j.ijsolstr.2020.05.014
• 25)   S.  Bin Zaman,  F.  Barlat and  J. H.  Kim: Int. J. Solids Struct., 134 (2018), 20. https://doi.org/10.1016/j.ijsolstr.2017.10.029
• 26)   S.-Y.  Lee,  J.-M.  Kim,  J.-H.  Kim and  F.  Barlat: Int. J. Mech. Sci., 183 (2020), 105769. https://doi.org/10.1016/j.ijmecsci.2020.105769
• 27)   M. S.  Wi,  S. Y.  Lee,  J. H.  Kim,  J. M.  Kim and  F.  Barlat: Int. J. Mech. Sci., 182 (2020), 105770. https://doi.org/10.1016/j.ijmecsci.2020.105770
• 28)   R. K.  Boger,  R. H.  Wagoner,  F.  Barlat,  M. G.  Lee and  K.  Chung: Int. J. Plast., 21 (2005), 2319. https://doi.org/10.1016/j.ijplas.2004.12.002
• 29)   F.  Barlat,  J. C.  Brem,  J. W.  Yoon,  K.  Chung,  R. E.  Dick,  D. J.  Lege,  F.  Pourboghrat,  S. H.  Choi and  E.  Chu: Int. J. Plast., 19 (2003), 1297. https://doi.org/10.1016/S0749-6419(02)00019-0
• 30)   A. V.  Hershey: J. Appl. Mech., 21 (1954), 241.
• 31)   J. Y.  Lee,  J. W.  Lee,  M. G.  Lee and  F.  Barlat: Int. J. Solids Struct., 49 (2012), 3562. https://doi.org/10.1016/j.ijsolstr.2012.03.042
• 32)   S.-Y.  Yoon,  S.-Y.  Lee and  F.  Barlat: Comput. Methods Appl. Mech. Eng., 372 (2020), 113449. https://doi.org/10.1016/j.cma.2020.113449
• 33)   J.  Lee,  H. J.  Bong,  J.  Ha,  J.  Choi,  F.  Barlat and  M.-G.  Lee: JOM, 70 (2018), 1560. https://doi.org/10.1007/s11837-018-2910-4
• 34)   J. Y.  Lee,  M. G.  Lee,  F.  Barlat,  K. H.  Chung and  D. J.  Kim: Int. J. Solids Struct., 87 (2016), 254. https://doi.org/10.1016/j.ijsolstr.2016.01.020
• 35)   X.  Xue,  J.  Liao,  G.  Vincze,  A. B.  Pereira and  F.  Barlat: Int. J. Mech. Sci., 117 (2016), 1. https://doi.org/10.1016/j.ijmecsci.2016.08.003
• 36)   S.  Sumikawa,  A.  Ishiwatari,  J.  Hiramoto,  F.  Yoshida,  T.  Clausmeyer and  A. E.  Tekkaya: Procedia Eng., 207 (2017), 179. https://doi.org/10.1016/j.proeng.2017.10.758
• 37)   S.  Sumikawa,  A.  Ishiwatari,  J.  Hiramoto and  T.  Urabe: J. Mater. Process. Technol., 230 (2016), 1. https://doi.org/10.1016/j.jmatprotec.2015.11.004
• 38)   K. Y.  Seo,  J. H.  Kim,  H. S.  Lee,  J. H.  Kim and  B. M.  Kim: Metals, 8 (2018), 18. https://doi.org/10.3390/met8010018
• 39)   J.  Choi,  J.  Lee,  H. J.  Bong,  M. G.  Lee and  F.  Barlat: Int. J. Solids Struct., 151 (2018), 152. https://doi.org/10.1016/j.ijsolstr.2017.09.030
• 40)   Y.  Choi,  J.  Lee,  S. S.  Panicker,  H. K.  Jin,  S. K.  Panda and  M. G.  Lee: Int. J. Mech. Sci., 170 (2020), 105344. https://doi.org/10.1016/j.ijmecsci.2019.105344
• 41)   F.  Yoshida and  T.  Uemori: Int. J. Mech. Sci., 45 (2003), 1687. https://doi.org/10.1016/j.ijmecsci.2003.10.013
• 42)   S. L.  Zang,  M. G.  Lee and  J. H.  Kim: Int. J. Mech. Sci., 77 (2013), 194. https://doi.org/10.1016/j.ijmecsci.2013.09.033
• 43)   J. Y.  Lee,  F.  Barlat and  M. G.  Lee: Int. J. Plast., 71 (2015), 113. https://doi.org/10.1016/j.ijplas.2015.04.005
• 44)   J. H.  Choi,  S. L.  Zang,  M. G.  Lee,  J. H.  Kim and  F.  Barlat: Int. J. Mech. Sci., 157–158 (2019), 428. https://doi.org/10.1016/j.ijmecsci.2019.04.050
• 45)   F.  Yoshida and  T.  Amaishi: Int. J. Plast., 130 (2020), 102708. https://doi.org/10.1016/j.ijplas.2020.102708
• 46)   R. W.  Logan and  W. F.  Hosford: Int. J. Mech. Sci., 22 (1980), 419.
• 47)   T.  Kuwabara,  T.  Mori,  M.  Asano and  T.  Hakoyama: Int. J. Plast., 93 (2017), 164. https://doi.org/10.1016/J.IJPLAS.2016.10.002
• 48)   M. G.  Lee,  D.  Kim,  C.  Kim,  M. L.  Wenner,  R. H.  Wagoner and  K.  Chung: Int. J. Plast., 23 (2007), 1189. https://doi.org/10.1016/j.ijplas.2006.10.011
• 49)   K.  Ahn,  D.  Yoo,  M. H.  Seo,  S. H.  Park and  K.  Chung: Met. Mater. Int., 15 (2009), 637. https://doi.org/10.1007/s12540-009-0637-z
• 50)   J.  Choi,  J.  Lee,  G.  Bae,  F.  Barlat and  M. G.  Lee: JOM, 68 (2016), 1850. https://doi.org/10.1007/s11837-016-1924-z
• 51)   J.  Jung,  S.  Jun,  H. S.  Lee,  B. M.  Kim,  M. G.  Lee and  J. H.  Kim: Metals, 7 (2017), 480. https://doi.org/10.3390/met7110480
• 52)   V. C.  Tong and  D. T.  Nguyen: Adv. Mech. Eng., 10 (2018), 1. https://doi.org/10.1177/1687814018797436
• 53)   J.  Lee,  J. Y.  Lee,  F.  Barlat,  R. H.  Wagoner,  K.  Chung and  M. G.  Lee: Int. J. Plast., 45 (2013), 140. https://doi.org/10.1016/j.ijplas.2013.01.011
• 54)   E. H.  Lee,  T. B.  Stoughton and  J. W.  Yoon: Int. J. Plast., 98 (2017), 217. https://doi.org/10.1016/j.ijplas.2017.08.003
• 55)   J. Y.  Lee,  M. G.  Lee,  F.  Barlat and  G.  Bae: Int. J. Plast., 93 (2017), 112. https://doi.org/10.1016/j.ijplas.2017.02.004
• 56)   S.  Sumikawa,  A.  Ishiwatari and  J.  Hiramoto: J. Mater. Process. Technol., 241 (2017), 46. https://doi.org/10.1016/j.jmatprotec.2016.11.005
• 57)   T.  Hasegawa and  T.  Yakou: Scr. Metall., 14 (1980), 1083. https://doi.org/10.1016/0036-9748(80)90210-0
• 58)   W. A.  Spitzig,  R. J.  Sober and  O.  Richmond: Metall. Trans. A, 7 (1976), 1703.

Appendix

State Variable Evolution Equations

All the state variables defined in the HAH₂₀ model as well as their initial values are listed in Table A1 and the associated coefficients in Table A2. The state variable evolution equations depend on $\lambda =sgn(\hat{\mathbf{s}}:\hat{\mathbf{h}})=sgn(cos\chi )$ and ${\lambda }^{\prime }=sgn(\hat{\mathbf{s}}:{\hat{\mathbf{h}}}^{\prime })=sgn(cos{\chi }^{\prime })$.

Table A1. State variables and initial values.

Variable	g−, g+	g3 , gP , gS	gL	gC	$\hat{\mathbf{h}}$, ${\hat{\mathbf{h}}}^{\prime }$‡
Initial value	1.	1.	1.	1.	${\hat{\mathbf{h}}}^{\prime }=\hat{\mathbf{h}}={\hat{\mathbf{h}}}_0$†
Note	Bauschinger effect	Permanent softening	Latent hardening	Cross- loading contraction	Microstructure deviator

^‡  Normalized tensor, e.g. $\hat{\mathbf{h}}$:$\hat{\mathbf{h}}$ = 1

^†  Stress deviator associated with first plastic strain increment

Table A2. HAH₂₀ model coefficients with suggested values and approximate ranges. Coefficients marked with arrow have to be determined. Coefficients without arrow may be set to suggested values to simplify calibration.

Coefficient		Suggested	Range‡	Remark
q		3	2–4	Yield condition
p		3	2–4
ξC		4	~4
ξL		0.5	0–1
k		k1	~k1	Microstructure deviator
k′		k1/2	~k1/2
ξR		8	4–12
k1	←	100	50–200	Bauschinger effect
k2	←	75	50–200
k3	←	0.4	0.1–0.8
ξB		4	1–8
${\xi }_B^{\prime }$		1.5	1–3
k4	←	0.9	0.6–1	Permanent softening
k5	←	15	5–20
kS		k1	k1/2−2k1
ξS	←	2	0.5–5
C	←	0.7	0.5–1	Cross-loading contraction
kC	←	60	10–200
$k_C^{\prime }$	←	80	10–200
${\xi }_C^{\prime }$		6	2–8
L	←	1.5	1–5	Latent hardening
kL	←	30	10–300
${\xi }_L^{\prime }$		0.25	0.05–2.0

^‡  Range is approximate

Microstructure Deviators

  

$$\frac{d{\hat{\mathbf{h}}}^{\prime }}{d\overline{\epsilon }}={\lambda }^{\prime }{k}^{\prime }\left(\hat{\mathbf{s}}-cos{\chi }^{\prime }{\hat{\mathbf{h}}}^{\prime }\right)$$(19)

  

$$\frac{d\hat{\mathbf{h}}}{d\overline{\epsilon }}=\lambda k\left(cos{{\chi }^{\prime }}^{2{\xi }_R}+cos{\chi }^{2{\xi }_R}\right)\left(\hat{\mathbf{s}}-cos\chi \hat{\mathbf{h}}\right)$$(20)

Bauschinger Effect

  

$$\begin{array}{l}{\displaystyle \frac{d{g}_{-}}{d\overline{\epsilon }}}=\\ {\displaystyle \frac{1+\lambda }{2}}\left\{k_1{\displaystyle \frac{1-g_{-}}{g_{-}^{{{\xi }^{\prime }}_B}}}\left[1-cos{\chi }^{2{\xi }_B}\right]+k_2\left(k_3{\displaystyle \frac{{\sigma }_{\text{y}}}{\sigma \left(\overline{\epsilon }\right)}}-g_{-}\right)cos{\chi }^{2{\xi }_B}\right\}\\ +{\displaystyle \frac{1-\lambda }{2}}\left\{k_1{\displaystyle \frac{1-g_{-}}{g_{-}^{{{\xi }^{\prime }}_B}}}\right\}\end{array}$$(21)

  

$$\begin{array}{l}{\displaystyle \frac{d{g}_{+}}{d\overline{\epsilon }}}=\\ {\displaystyle \frac{1-\lambda }{2}}\left\{k_1{\displaystyle \frac{1-g_{+}}{g_{+}^{{{\xi }^{\prime }}_B}}}\left[1-cos{\chi }^{2{\xi }_B}\right]+k_2\left(k_3{\displaystyle \frac{{\sigma }_{\text{y}}}{\sigma \left(\overline{\epsilon }\right)}}-g_{+}\right)cos{\chi }^{2{\xi }_B}\right\}\\ +{\displaystyle \frac{1+\lambda }{2}}\left\{k_1{\displaystyle \frac{1-g_{+}}{g_{+}^{{{\xi }^{\prime }}_B}}}\right\}\end{array}$$(22)

Permanent Softening

At the beginning of the deformation, and thereafter, at the exact moment the sign of cosχ changes, the state variables g_P, g₃ and $\hat{\mathbf{h}}$ are memorized as $g_P^{*}$, $g_3^{*}$ and ${\hat{\mathbf{h}}}^{*}$, respectively. An intermediate variable denoted $g_P^{\prime }$ is calculated   

$$g_P^{\prime }=g_P^{*}-\left(g_P^{*}-g_3^{*}\right){\left|\mathbf{s}:{\hat{\mathbf{h}}}^{*}\right|}^{2{\xi }_S}$$(23)

which becomes the new value of g_P if $\mathrm{\Delta }{g}_P=g_P^{\prime }-g_P<0$. The evolution of the state variable g₃ is   

$$\frac{d{g}_3}{d\overline{\epsilon }}=k_5{g}_S\left(k_4-g_3\right)$$(24)

The Intermediate Variable $g_3^{\prime }$

  

$$g_3^{\prime }=g_3+d{g}_3$$(25)

becomes the new value of g₃ if $g_3^{\prime }-g_P\le 0$. Otherwise, g₃ is set equal to g_P. Finally, the variable g_S is introduced to decrease the amount of permanent softening evolution when the material is not subjected to pure reversal.   

$$\frac{d{g}_S}{d\overline{\epsilon }}=-k_S\left[1-\left|\mathbf{s}:{\hat{\mathbf{h}}}^{*}\right|\right]g_S$$(26)

Cross-loading Contraction

  

$$\frac{d{g}_C}{d\overline{\epsilon }}=k_C\left(C-g_C\right)$$(27)

or   

$$\frac{d{g}_C}{d\overline{\epsilon }}=k_C^{\prime }\frac{1-g_C}{g_C^{{\xi }_C^{\prime }}}$$(28)

Equation (27) is used only when |cosχ| = 1.0, but practically when |cosχ| ≥ 0.996, and Eq. (28) in all the other cases.

Latent Hardening

  

$$\begin{array}{l}{\displaystyle \frac{d{g}_L}{d\overline{\epsilon }}}=\\ k_L\left[{\displaystyle \frac{{\sigma }_r\left(\overline{\epsilon }\right)-{\sigma }_r\left(0\right)}{{\sigma }_r\left(\overline{\epsilon }\right)}}\left(\sqrt{L+\left(1-L\right)cos{{\chi }^{\prime }}^{2{\xi }_L^{\prime }}}-1\right)+1-g_L\right]\end{array}$$(29)

Application

Figure A1(a) shows the deviatoric plane with s_xx, s_yy and s_zz the corresponding normal components in the RD, TD and ND, respectively. This figure describes how the yield surface distorts during a pre-strain in uniaxial tension in the TD. The dash and solid lines represent the initial and distorted anisotropic yield surfaces, respectively. The distortion of the yield surface by flattening of the side opposite to the current loading (s) occurs as g₋ decreases from 1 to a lower value dictated by its evolution Eq. (21), leading to the Bauschinger effect (BE). The flow stress is σ_t in uniaxial tension but changes to σ_c in compression if the load is suddenly reversed. The ratio of these two stresses is given by   

$$\frac{{\sigma }_c}{{\sigma }_t}={\left(\frac{1}{f_{-}^q+1}\right)}^{1/q}$$(30)

which can be calculated with Eq. (11) and the current value of g₋. In addition to the BE effect that develops during TD uniaxial tension, the yield surface contracts in all the directions orthogonal to $\hat{\mathbf{h}}$ as g_C decreases from 1 to a lower value as dictated by its evolution Eq. (27), leading to cross-loading contraction (CC). Figure A.1(b) describes the first plastic increment after a change from TD uniaxial tension to the loading corresponding to the new loading direction s. In a first case, g_C = 1 leading to s_X = s_X1 = 0, and s_L = s_L1, which is the stress state if no distortion occurs. In a second case, g_C < 1 leading to s_X = s_X2 ≠ 0 and s_L = s_L2 = s, which is the new stress state on the distorted yield surface. Note that s_X2 does influence the value of the yield function, Eq. (9) combined with Eq. (13), that is, the length of s.

Fig. A1.

Normalized yield surfaces in deviatoric plane assuming isotropic hardening (dash line) and distortional plasticity (solid line) (a). During TD uniaxial tension pre-strain (b). During the first strain increment after pre-strain and reloading in direction s. (Online version in color.)