2024 Volume 65 Issue 10 Pages 1193-1217
Forming simulation is an indispensable analysis tool in industry, the most important objective of which being the reproduction of the material deformation behavior during the process as accurately as possible to predict forming defects and determine optimum forming conditions precisely. This paper reviews the material models suitable for metal forming simulations and the material test methods that can reproduce the various types of stress states occurring in real forming operations. These test methods are essential to verify the validity of material models. Examples of forming simulation results for aluminum alloys are presented, illustrating the significant influence of the material models on the accuracy of the process predictions. Finally, the emerging research trends in the field of material modeling and testing are discussed.
Simulation is now an indispensable analysis tool for predicting defects in forming processes, such as fracture and springback, and determining the optimal manufacturing conditions. The most important objective of a simulation is to reproduce the deformation behavior of the material during a forming operation as accurately as possible. Aluminum alloys (AAs) are prone to fracture and exhibit large springback when the forming loads are removed. Therefore, it is of prime importance to predict forming defects in advance, and determine the optimal processing conditions through forming simulations.
Many factors affect the deformation behavior of materials during forming [1–3]. If the influencing factors in the real world were perfectly reproduced in the virtual space, i.e., forming simulations, the calculated and experimental results should be identical. However, one of the factors in which there is a large discrepancy between the real world and forming simulations is the material model because it determines the strain and strain rate occurring in response to a given set of stresses and boundary conditions. The development of a perfect material model would be achieved only if the discrepancy between the actual elastoplastic deformation characteristics of the material and those predicted by the model was nil.
In order to guarantee the accuracy of material models, it is also necessary to develop material testing methods to confirm their validity. The simplest material test method is uniaxial tension. However, in real sheet forming processes, the materials are subjected to multiaxial stress states, very often combined with stress reversals or other non-linear loading modes. Therefore, it is important that material test methods reproduce such complex stress states.
This review paper consists mostly of four sections. Section 2 introduces several basic concepts necessary for material testing and modeling. Section 3 reviews the material models suitable for metal forming simulations, i.e., yield conditions and their evolution as a function of the accumulated plastic strain, with more emphasis for sheet product applications. Section 4 reviews the material test methods suitable at reproducing the various types of stress states occurring in real forming operations. These test methods are essential to verify the validity of material models. In addition, examples of test data and material modeling results for AAs are presented. Section 5 shows examples of finite element analysis (FEA) results for AAs, illustrating the significant influence of the material model on the accuracy of the forming simulations. Finally, the emerging research trends in the field of material modeling and testing are discussed.
This section reviews several basic concepts that are essential to a good understanding of material testing and modeling. Experimental plasticity is too broad of a subject for an exhaustive review in a manuscript of this size. The reader is also referred to the excellent reviews of earlier work on experimental plasticity [4–13].
2.1 Yield surface and contour of equal plastic workThe yield surface is a closed surface in stress space consisting of a set of multiaxial stress states. The material behavior is purely elastic or elasto-plastic depending whether the stress state remains strictly inside or on the yield surface, respectively. The yield surface can be measured by applying as many linear and non-linear loading modes as possible to a material in different directions, recording the stress marking the transition between elastic and plastic domains, i.e., the yield stress, and plotting them in stress space.
In the following, the x-, y-, and z-axes denote the rolling (RD), rolling transverse (TD), and thickness normal (ND) directions of a sheet metal, respectively. Assuming that the sheet material is subjected to a plane stress state with components σxx, σyy, and σxy, while σzz = σzx = σzy = 0, the yield surface can be schematically represented as shown in Fig. 1. This figure also depicts the schematic diagrams of material test methods suitable for reproducing various stress states on the yield surface. An overview of each test method and related literature are provided in Section 4.
Schematic diagram of a yield surface in σxx − σyy − σxy stress space and schematic representation of material test methods associated with specific stress states. The angles along Line F represent θTC, the angle of the maximum principal stress direction from the RD of the test material.
It is an arduous task to properly measure a yield surface experimentally because it is difficult to clearly identify the starting point of plastic deformation [5], except for low carbon steel, in which the yield point is clearly defined. In addition, when determining the yield surface of a material that was work-hardened by a given prestrain along a certain strain path/stress history (referred to as the subsequent yield surface), it is necessary to conduct experiments to measure the entire yield surface for this prestrain [6, 8]. Furthermore, if the shape of the subsequent yield surface changes depending on the magnitude of the prestrain [14, 15], the number of experiments for determining subsequent yield surfaces becomes astronomical. Yet, if the shape of the subsequent yield surface changes with an increase of the prestrain [6], it may be doubtful to justify the use of a material model created from the initial yield surface data to analyze metal forming problems involving large plastic strains.
Therefore, Hill et al. [16, 17] assumed that the evolution of subsequent yield surfaces with work-hardening is equivalent to that of contours of equal plastic work per unit volume (hereafter referred to as work contour) in stress space. A work contour is experimentally determined as follows: The true uniaxial stress − logarithmic plastic strain curve ($\sigma_{0} - \varepsilon_{0}^{\text{p}}$ curve) measurements in the RD is assumed to be a reference data for work hardening. Namely, the plastic work per unit volume, $w_{0}^{\text{p}}$, associated with particular values of $\varepsilon_{0}^{\text{p}}$ is determined. The stress points that give the same plastic work as $w_{0}^{\text{p}}$ for many linear strain paths form a work contour associated with $\varepsilon_{0}^{\text{p}}$ in stress space.
The present authors have performed biaxial stress tests on various metal sheets and tubes, and found that, at least for linear stress paths, the work contours can be regarded as plastic potentials; the outward normal vectors to the work contours are almost parallel to the direction of the plastic strain rate, Dp. Furthermore, it has been confirmed that the method of constructing a material model based on the measured work contours and the directions of Dp is effective in improving the accuracy of forming simulations [12].
In the following, unless otherwise specified, the terms “work contour” and “yield surface” are used interchangeably.
2.2 Isotropic hardening, differential hardening, and anisotropic hardeningAs $w_{0}^{\text{p}}$, or equivalently $\varepsilon_{0}^{\text{p}}$, increases, the material work-hardens and the size of the work contour becomes larger. When the work contour expands while maintaining the same shape, it is referred to as isotropic hardening (IH) while, when the shape changes, it is referred to as differential hardening (DH) [16, 17]. Many metals, including aluminum alloys, exhibit DH to varying degrees, but some materials almost exhibit IH [18].
Generally, the yield surface evolution is quite complex. More accurate models consist in combinations of expansion, translation, distortion and rotation of the yield surface [19]. This non-isotropic hardening behavior is referred to as anisotropic hardening. The anisotropic hardening effects are particularly important for the simulations of forming involving non-linear strain paths as it occurs in many processes. On the other hand, DH is more reserved for proportional loading.
Beyond the elastic range, metallic materials exhibit plastic, irreversible deformations. At each point of a continuum, the yield surface in stress space separates pure elastic and elasto-plastic states. A stress state is always inside (elastic state) or on (elasto-plastic state) the yield surface but never outside. The equation defining this surface is called the yield condition and writes
| \begin{equation} f(\boldsymbol{\sigma} ,\mathbf{v}) \leq 0 \end{equation} | (1) |
where f is the yield function, σ the Cauchy stress tensor, and v symbolically represents a number of internal variables, scalar or tensorial, that characterize the microstructural state of the material. The stress state is elastic when f < 0 and elasto-plastic when f = 0. Because of the plastic deformation mechanisms of polycrystals, the yield surface should be convex in stress space and be a potential for the plastic strain increment, i.e., dεp = Dpdt, where dt is an increment of time [20]. This indicates that the strain increment should be orthogonal to the yield surface at the loading point which, mathematically, translates in
| \begin{equation} d\varepsilon_{mn}^{p} = d\lambda \frac{\partial f(\boldsymbol{\sigma} ,\mathbf{v})}{\partial \sigma_{mn}} \end{equation} | (2) |
The plastic multiplier dλ controls the size of the deformation increment. Equation (2) represents what is called the associated flow rule (AFR) but non-associated (NAFR) flow rules, where this orthogonality does not exist, have also been adopted by authors, e.g. Stoughton [21]. In this case, a plastic potential, different from the yield function must be introduced. Experimental results [4–6, 12] indicate that the AFR is quite suitable for metals. The plastic part of the strain increment dεp is a deviatoric tensor as a consequence of plastic incompressibility observed in deforming metals, a condition expressed as
| \begin{equation} d\varepsilon_{11}^{p} + d\varepsilon_{22}^{p} + d\varepsilon_{33}^{p} = 0 \end{equation} | (3) |
Because of the multiple mechanisms involved in plastic deformation, the evolution of the yield condition as a function of the accumulated plastic strain may be very complex. Examples of this complexity will be introduced in a later section with succinct descriptions of kinematic and distortional hardening evolutions, necessary when strain path changes occur. However, isotropic hardening, for which the yield surface expands only, is a reasonable assumption for forming problems and will be first examined in detail. Assuming isotropic hardening, the yield condition may write as
| \begin{equation} f(\boldsymbol{\sigma} ,\bar{\varepsilon}) = \bar{\sigma}(\boldsymbol{\sigma} ) - \sigma_{r}(\bar{\varepsilon}) \leq 0 \end{equation} | (4) |
where $\bar{\sigma }(\boldsymbol{\sigma})$ is an effective stress and $\sigma_{r}(\bar{\varepsilon })$ a reference flow curve. The effective stress lumps all the components of the stress tensor into a single function $\bar{\sigma }$ while the effective strain, defined incrementally, lumps all the components of the plastic strain increment tensor into a variable $d\bar{\varepsilon }$. The function $\bar{\sigma }$ and the incremental variable $d\bar{\varepsilon }$ may be those defined as the von Mises equivalent stress and equivalent strain increment for instance. The flow stress may be chosen as a uniaxial tension curve with strain ε and modeled by, for instance, a Hollomon hardening model, σr(ε) = εn. However, these are only examples, as the choice of each element of the set $\bar{\sigma }$, $d\bar{\varepsilon }$ and σr might be totally different. For instance, a more general definition of $d\bar{\varepsilon }$ is the plastic work-based effective strain, i.e.,
| \begin{equation} d\bar{\varepsilon} = \frac{\sigma_{mn}d\varepsilon_{mn}^{p}}{\bar{\sigma}} = \frac{\sigma_{mn}d\varepsilon_{mn}^{p}}{\sigma_{r}} \end{equation} | (5) |
that reduces to von Mises equivalent strain increment when $\bar{\sigma }$ is the von Mises effective stress. Note that, in general, an effective stress has the unit of a pressure and, therefore, is a mathematically homogeneous function of first degree, that is, $\bar{\sigma }(k\boldsymbol{\sigma}) = k\bar{\sigma }(\boldsymbol{\sigma})$ with an arbitrary scalar k. Thus, combining the associated flow rule with the plastic-work effective strain leads to
| \begin{equation} d\bar{\varepsilon} = \frac{\sigma_{mn}d\varepsilon_{mn}^{p}}{\bar{\sigma}} = \frac{\sigma_{mn}}{\bar{\sigma}}d\lambda \frac{\partial \bar{\sigma}(\boldsymbol{\sigma} )}{\partial \sigma_{mn}} = d\lambda \end{equation} | (6) |
which, through the use of Euler theorem on homogeneous function, i.e.,
| \begin{equation} \sigma_{mn}\frac{\partial \bar{\sigma}(\boldsymbol{\sigma} )}{\partial \sigma_{mn}} = \bar{\sigma} \end{equation} | (7) |
reduces $d\bar{\varepsilon }$ to dλ. Therefore, based on the AFR, the definition of the theory of plasticity for isotropic hardening is completely defined when the effective stress $\bar{\sigma }$ and the reference flow stress σr are selected. The former gives the shape of the yield surface in stress space (e.g., von Mises, Tresca, etc.). The latter determines the rate at which it expands and corresponds to the best approximation of an experimental flow curve for a given stress state, usually uniaxial tension (one normal stress only) or balanced biaxial tension (two equal normal stresses only). The following discussion focused on the description of $\bar{\sigma }$ because the expressions and applications of σr are more straightforward.
3.2 Isotropic hardeningSince it is easier to model isotropic hardening for isotropic materials, this case is examined first, and the extension of the theory to plastic anisotropy second.
3.2.1 Stress invariantsFor an isotropic material, $\bar{\sigma }$ is independent of the chosen reference set of unit and mutually orthogonal base vectors, say, e1, e1 and e3 in which the tensor components are expressed. This indicates that no matter what the set of base vectors at each point is, the isotropic yield condition yields to the same plastic state. Hence, the yield condition for an isotropic material should be a function of stress tensor invariants only. The Cauchy stress tensor σ has three independent invariants but various combinations of such invariants are possible to express the yield condition (e.g., Życzkowski [22]). Solving the cubic equation
| \begin{equation} \mathcal{P}(\mu) = \det (\boldsymbol{\sigma} - \mu \mathbf{I}) = \begin{vmatrix} \sigma_{11} - \mu & \sigma_{12} & \sigma_{31}\\ \sigma_{12} & \sigma_{22} - \mu & \sigma_{23}\\ \sigma_{31} & \sigma_{23} & \sigma_{33} - \mu \end{vmatrix} = 0 \end{equation} | (8) |
where the characteristic cubic polynomial equation writes
| \begin{equation} \mathcal{P}(\mu) = - \mu^{3} + 3\tilde{I}_{1}\mu^{2} + 3\tilde{I}_{2}\mu + 2\tilde{I}_{3} = 0 \end{equation} | (9) |
leads to the principal values of the stress tensor σ1, σ2 and σ3. Once the principal values are found, the principal directions are obtained through a standard analysis that is not further discussed here. The coefficients of eq. (9) are called the principal invariants
| \begin{align} I_{1} & = 3\tilde{I}_{1} = \sigma_{11} + \sigma_{22} + \sigma_{33}\\ I_{2} & = 3\tilde{I}_{2} = \sigma_{23}^{2} + \sigma_{31}^{2} + \sigma_{12}^{2} - \sigma_{22}\sigma_{33} - \sigma_{33}\sigma_{11} - \sigma_{11}\sigma_{22}\\ I_{3} & = 2\tilde{I}_{3} = 2\sigma_{23}\sigma_{31}\sigma_{12} + \sigma_{11}\sigma_{22}\sigma_{33} - \sigma_{11}\sigma_{23}^{2} \\ &\quad - \sigma_{22}\sigma_{31}^{2} - \sigma_{33}\sigma_{12}^{2} \end{align} | (10) |
$\tilde{I}_{1} = \text{tr}(\boldsymbol{\sigma})/3 = \sigma_{M}$ is also called the mean stress. Moreover, the stress deviator
| \begin{equation} \sigma'_{pq} = \sigma_{pq} - \delta_{pq}\sigma_{M} \end{equation} | (11) |
where δpq is the Kronecker symbol, should have its own invariants called J1, J2 and J3, with a different notation to indicate that they specifically correspond to a deviator. These quantities are calculated with eq. (10) as well but with the components of the stress deviator, leading to $J_{1} = 3\tilde{J}_{1} = 0$, $J_{2} = 3\tilde{J}_{2}$ and $J_{3} = 2\tilde{J}_{3}$. Since J1 = 0, the mean stress σM is added to J2 and J3 to form a system of three independent invariants to characterize the stress tensor. It can be shown that J2 is always non-negative. Based on the principal invariants, the so-called “angular invariant” is also introduced, i.e.,
| \begin{equation} \theta = \arccos \left[\frac{2\tilde{I}_{1}^{3} + 3\tilde{I}_{1}\tilde{I}_{2} + 2\tilde{I}_{3}}{2(\tilde{I}_{1}^{2} + \tilde{I}_{2})^{3/2}} \right] \end{equation} | (12) |
For a deviatoric stress tensor, the angular invariant reduces to
| \begin{equation} \theta' = \arccos\left(\frac{\tilde{J}_{3}}{\tilde{J}_{2}^{3/2}} \right) \end{equation} | (13) |
from which the ordered principal stresses σ1 ≥ σ2 ≥ σ3 and corresponding deviatoric values can be calculated (see Barlat et al. [23, 24]). These expressions are as follows
| \begin{align} &\left\{ \begin{array}{l} \sigma_{1} = 2\sqrt{\tilde{I}_{1}^{2} + \tilde{I}_{2}} \cos \left(\dfrac{\theta}{3} \right) + \tilde{I}_{1}\\ \sigma_{2} = 2\sqrt{\tilde{I}_{1}^{2} + \tilde{I}_{2}} \cos \left(\dfrac{\theta}{3} - \dfrac{2\pi}{3} \right) + \tilde{I}_{1}\\ \sigma_{3} = 2\sqrt{\tilde{I}_{1}^{2} + \tilde{I}_{2}} \cos \left(\dfrac{\theta}{3} - \dfrac{4\pi}{3} \right) + \tilde{I}_{1} \end{array} \right.\\ & \left\{ \begin{array}{l} \sigma_{1}' = 2\sqrt{\tilde{J}_{2}} \cos \left(\dfrac{\theta'}{3} \right)\\ \sigma_{2}' = 2\sqrt{\tilde{J}_{2}} \cos \left(\dfrac{\theta'}{3} - \dfrac{2\pi}{3} \right)\\ \sigma_{3}' = 2\sqrt{\tilde{J}_{2}} \cos \left(\dfrac{\theta'}{3} - \dfrac{4\pi}{3} \right) \end{array} \right. \end{align} | (14) |
Note that eq. (14) demonstrates that the principal stresses or deviatoric stresses consist of two sets of invariants. Although many invariant combinations can be used to express isotropic yield functions (Życzkowski [22]), it is the author’s opinion that the use of principal stresses is the most effective to extend yield condition from plastic isotropy to plastic anisotropy. This is explained in more details below.
3.2.2 Isotropic yield conditionsTresca [25] proposed the oldest yield condition by assuming that plastic yielding occurs when the maximum shear stress reaches a critical value, the simple shear flow stress. Simple considerations, based on the Mohr circle, lead to
| \begin{equation} \bar{\sigma} = \frac{\sigma_{1} - \sigma_{3}}{2} = \frac{\sigma_{1}' - \sigma_{3}'}{2} = \sigma_{s} = \frac{1}{2}\sigma_{r}(\bar{\varepsilon}) \end{equation} | (15) |
where, in this case, σs is the shear flow stress. If $\sigma_{r}(\bar{\varepsilon })$ is the uniaxial tension flow curve, then, σs = σr/2. Using the principal deviatoric stresses in eq. (14), Tresca yield condition rewrites as
| \begin{equation} \bar{\sigma} = \sqrt{J_{2}} \sin \left(\frac{2\pi}{3} - \frac{\theta'}{3} \right) = \frac{1}{2}\sigma_{r} \end{equation} | (16) |
The Tresca yield surface is a cylinder of axis defined by σ1 = σ2 = σ3 in principal stress space, with a hexagonal cross-section viewed in the deviatoric plane orthogonal to the cylinder axis. The most widely used yield condition, proposed by von Mises [26], is expressed with a quadratic function reducing to the corresponding effective stress written, in an expanded form, as
| \begin{equation} \bar{\sigma} = \left\{\frac{(\sigma_{1} - \sigma_{2})^{2} + (\sigma_{2} - \sigma_{3})^{2} + (\sigma_{3} - \sigma_{1})^{2}}{2} \right\}^{\frac{1}{2}} {}= \sigma_{u} \end{equation} | (17) |
where $\sigma_{r}(\bar{\varepsilon }) = \sigma_{u}(\bar{\varepsilon })$ in this case denotes the uniaxial flow stress. The corresponding yield surface is a cylinder with circular cross-section orthogonal to the axis defined by σ1 = σ2 = σ3, or an elliptical cross-section in a plane stress state, that is, when one of the three principal stresses is nil. Note that, since $\sigma '_{p} = \sigma_{p} - \tilde{I}_{1}$, with p = 1 to 3, von Mises yield condition writes as a function of the principal deviatoric stresses, as follows
| \begin{equation} \bar{\sigma} = \left\{\frac{(\sigma_{1}' - \sigma_{2}')^{2} + (\sigma_{2}' - \sigma_{3}')^{2} + (\sigma_{3}' - \sigma_{1}')^{2}}{2} \right\}^{\frac{1}{2}}{} = \sigma_{u} \end{equation} | (18) |
and as a function of the second stress deviator invariant as
| \begin{equation} \bar{\sigma} = \sqrt{9\tilde{J}_{2}} = \sqrt{3J_{2}} = \sigma_{u} \end{equation} | (19) |
The von Mises yield condition was further extended by Drucker [27] to
| \begin{equation} \bar{\sigma} = \sqrt{3} [J_{2}^{3} - cJ_{3}^{2}]^{1/6} = \sigma_{u} \end{equation} | (20) |
where c is a coefficient. Another description of a cylindrical yield surface, which includes Tresca and von Mises as particular cases, was proposed by Hershey [28] to represent yield surfaces calculated with a self-consistent crystal plasticity model
| \begin{equation} \bar{\sigma} = \left\{\frac{| \sigma_{1} - \sigma_{2} |^{a} + | \sigma_{2} - \sigma_{3} |^{a} + | \sigma_{3} - \sigma_{1} |^{a}}{2} \right\}^{\frac{1}{a}} {}= \sigma_{u} \end{equation} | (21) |
For a = 2 or 4, eq. (21) reduces to the Mises yield condition, whereas for a = 1 or +∞, it leads to Tresca. Hershey’s yield condition may be expressed with principal deviatoric stresses only
| \begin{equation} \bar{\sigma} = \left\{\frac{| \sigma_{1}' - \sigma_{2}' |^{a} + | \sigma_{2}' - \sigma_{3}' |^{a} + | \sigma_{3}' - \sigma_{1}' |^{a}}{2} \right\}^{\frac{1}{a}}{} = \sigma_{u} \end{equation} | (22) |
and, using $J_{2} = 3J'_{2}$ and θ′, as
| \begin{equation} \bar{\sigma} = 2\sqrt{J_{2}} \left\{\cfrac{\biggl| \sin \biggl(\cfrac{\pi}{3} - \cfrac{\theta'}{3} \biggr) \biggr|^{a} + \biggl| \sin \biggl(- \cfrac{\theta'}{3} \biggr) \biggr|^{a} + \biggl| \sin \biggl(\cfrac{2\pi}{3} - \cfrac{\theta'}{3} \biggr) \biggr|^{a}}{2} \right\}^{\frac{1}{a}}{} = \sigma_{u} \end{equation} | (23) |
A yield condition should be convex because it is an important property obtained from the physics of plastic deformation, and an essential requirement for stable numerical simulation of processes. It is relatively simple to show that the Hershey yield condition is convex with eq. (21) or (22). However, it is much more difficult with eq. (23). The point is that, although isotropic yield conditions based on invariants such as J2, J3 or θ′ are correct, their expressions based on principal stresses is more convenient.
According to Hosford [29], a particularly good agreement with many crystal plasticity results can be obtained with eq. (21). Karafillis and Boyce [30] proposed a generalization of Hershey’s model as
| \begin{align} \bar{\sigma} &= \bigg\{(1 - c)\frac{| \sigma_{1}' - \sigma_{2}' |^{a} + | \sigma_{2}' - \sigma_{3}' |^{a} + | \sigma_{3}' - \sigma_{1}' |^{a}}{2}\\ &\quad + 3^{a}c\frac{| \sigma_{1}' |^{a} + | \sigma_{2}' |^{a} + | \sigma_{3}' |^{a}}{2^{a} + 2} \bigg\}^{\frac{1}{a}} = \sigma_{u} \end{align} | (24) |
where c, different from that in eq. (20), is a coefficient between 0 and 1. Cazacu and Barlat [31] introduced an isotropic yield function that accounts for the strength differential (SD) effect in pressure-independent materials such as HCP polycrystals
| \begin{align} \bar{\sigma} &= \left\{\frac{|| \sigma_{1}' | - k\sigma_{1}' |^{a} + || \sigma_{2}' | - k\sigma_{2}' |^{a} + || \sigma_{3}' | - k\sigma_{3}' |^{a}}{2} \right\}^{\frac{1}{a}}{} \\ &= \sigma_{u} \end{align} | (25) |
where k is a coefficient that characterizes the tension-compression asymmetry. Some authors have proposed isotropic, pressure-dependent, yield conditions to explain the SD effect in cubic metals. Based on experimental results obtained in a pressure chamber (Spitzig et al. [32]), Richmond and Spitzig [33] proposed the following yield condition
| \begin{equation} \bar{\sigma}(\boldsymbol{\sigma} ) + \alpha cI_{1}(\boldsymbol{\sigma} ) = c \end{equation} | (26) |
where c is a flow stress and α a pressure coefficient equal to 20 TPa−1 for steel and 55 TPa−1 for aluminum. The literature on yield conditions is wide as, for instance, Życzkowski [22] and Yu [34] listed several hundreds.
3.2.3 Anisotropic yield conditions a. QuadraticAnisotropy is the dependence of material properties on the testing direction. Unlike isotropic theories for which the choice of a reference frame is arbitrary, anisotropic formulations must be expressed in a system of axes attached to the material. For instance, due to prior thermo-mechanical processing history, sheet metals exhibit orthotropic symmetry, i.e., they possess three mutually orthogonal planes of symmetry. In this paper, the unit vectors associated with these planes are denoted ex, ey and ez, corresponding to the sheet rolling (RD), transverse (TD) and normal (ND) directions, respectively. This means that the yield condition should always be evaluated with stress components expressed in the suitable unit vectors, i.e., σxx, σyz, etc. The first anisotropic yield condition was proposed by von Mises [35] as a generalization of his isotropic version. When it includes the additional requirement of independence of the mean stress σM, it contains 15 anisotropy coefficients. As a particular case, Hill [36] proposed the following quadratic form for orthotropic materials with six coefficients, F, G, H, L, M and N
| \begin{align} \bar{\sigma}& = \{F(\sigma_{yy} - \sigma_{zz})^{2} + G(\sigma_{zz} - \sigma_{xx})^{2} + H(\sigma_{xx} - \sigma_{yy})^{2}\\ &\quad + 2L\sigma_{yz}^{2} + 2M\sigma_{zx}^{2} + 2N\sigma_{xy}^{2}\}^{\frac{1}{2}} = \sigma_{U} \end{align} | (27) |
In the case of a sheet exhibiting planar isotropy that is, the properties measured in any direction of the sheet plane are identical, it is always possible to adopt the principal stress axes as references. Thus, in plane stress, Hill’s yield function reduces to either
| \begin{align} i) &\ \bar{\sigma} = \left\{\sigma_{1}^{2} - \left[2 - \left(\frac{\sigma_{u}}{\sigma_{b}} \right)^{2} \right]\sigma_{1}\sigma_{2} + \sigma_{2}^{2} \right\}^{\frac{1}{2}} {}= \sigma_{u}\quad \text{or}\\ \textit{ii}) &\ \bar{\sigma} = \left\{\sigma_{1}^{2} - \frac{2r}{1 + r}\sigma_{1}\sigma_{2} + \sigma_{2}^{2} \right\}^{\frac{1}{2}} {}= \sigma_{u} \end{align} | (28) |
where the r-value, also called Lankford coefficient, is the width-to-thickness strain ratio in uniaxial tension. Note that these two yield conditions do not lead to the same results, depending on which parameter, σb or r, is considered. From eq. (28), the balanced biaxial stress σb is of magnitude $\sigma_{u}\sqrt{(1 + r)/2} $. Woodthorpe and Pearce [37] found that the yield stress in balanced biaxial tension for aluminum alloy sheet having a r-value between 0.5 and 0.65 is significantly higher than the uniaxial yield stress in the plane of the sheet. However, Hill’s criterion cannot describe this behavior, i.e., materials with r < 1 and σb > σu. To capture this so-called “anomalous” behavior, non-quadratic yield formulations were considered.
b. Non-quadraticHill [38] proposed a non-quadratic form of the yield condition from which four special cases were extracted. The most useful expression of this yield criterion is the so-called “Special Case IV”, which applies to materials exhibiting planar isotropy
| \begin{equation} \bar{\sigma} = \left\{\frac{| \sigma_{1} + \sigma_{2} |^{a} + (1 + 2r) | \sigma_{1} - \sigma_{2} |^{a}}{2(1 + r)}\right\}^{\frac{1}{a}} {}= \sigma_{u} \end{equation} | (29) |
and reduces to eq. (28) when a = 2. Although the “anomalous” behavior is captured with this function for a < 2, the predicted yield surface shapes are sometimes different from those either determined experimentally or predicted with crystal plasticity models. Bassani [39] developed a family of yield functions, which depend on four parameters, in order to approximate a relatively broad range of transversely (planar) isotropic yield surfaces obtained with the Bishop-Hill [20] crystal plasticity model. Logan and Hosford [40] proposed a particular expression of Hill’s general yield criterion [36] to describe crystal plasticity results
| \begin{align} \bar{\sigma} &= \{F|\sigma_{yy} - \sigma_{zz}|^{a} + G|\sigma_{zz} - \sigma_{xx}|^{a} + H|\sigma_{xx} - \sigma_{yy}|^{a}\}^{\frac{1}{a}}\\ & = \sigma_{u} \end{align} | (30) |
where F, G and H are material coefficients. It leads to a good approximation of yield loci computed using the Bishop-Hill crystal plasticity model by setting a = 6 for BCC materials and a = 8 for FCC materials (Logan and Hosford [40]). Hosford [41] included a shear stress component in eq. (30) for plane stress cases but found some inconsistencies in the formulation. Barlat and Lian [42] successfully extended Hershey’s criterion with the influence of the shear term for plane stress. Hill [43] proposed an alternate way to include a shear stress component in a yield function, which was slightly modified by Lin and Ding [44]. However, none of the yield functions presented captured the experimentally observed behavior of certain aluminum alloys and brass (Stout and Hecker [45]) for which the flow stresses and r-values in the RD and TD are such that σRD = σTD but rRD ≠ rTD, and reciprocally. To account for this so-called “anomalous behavior of second order”, Hill [46] proposed a criterion expressed based on a third order polynomial of the stresses. Other non-quadratic yield criteria were proposed by Gotoh [47], Zhou [48], Yoshida et al. [49] and others. Additional details can be found in Barlat et al. [50] and Banabic et al. [51, 52].
c. Linear transformationsThe concept of linear transformation initiated in Barlat et al. [23] and Karafillis and Boyce [30] allows a systematic transformation of convex isotropic yield functions, described using the principal (deviatoric) stresses such as those given in the range of eqs. (15) to (25), to convex anisotropic yield functions. In this case, the equations of the previous Sections and, notably the isotropic yield conditions, are fully part of the anisotropic theory. For instance, for a pressure independent material, the stress deviator can be transformed linearly as
| \begin{equation} \begin{split} & i)\ \tilde{\boldsymbol{\sigma} }' = \mathbf{C}{:}\boldsymbol{\sigma} ' = \mathbf{C}{:}\mathbf{T}{:}\boldsymbol{\sigma} = \mathbf{L}{:}\boldsymbol{\sigma} \quad \text{or, equivalently}\\ & \textit{ii})\ \tilde{\sigma}_{pq}' = C_{pqmn}\sigma_{mn}' = C_{pqmn}T_{mnrs}\sigma_{rs} = L_{pqrs}\sigma_{rs} \end{split} \end{equation} | (31) |
where C and L are fourth order tensors, expressed in the material symmetry axes, containing the anisotropy coefficients. The structure of such tensors is similar to those of the elastic moduli of crystals or polycrystals with the suitable symmetry. T transforms the Cauchy stress tensor to its deviator. Then, the principal values of $\tilde{\boldsymbol{\sigma}}'$ substitute the principal deviatoric stresses of the associated isotropic yield function as eq. (22) for instance, leading to
| \begin{equation} \bar{\sigma} = \left\{\frac{| \tilde{\sigma}_{1}' - \tilde{\sigma}_{2}' |^{a} + | \tilde{\sigma}_{2}' - \tilde{\sigma}_{3}' |^{a} + | \tilde{\sigma}_{3}' - \tilde{\sigma}_{1}' |^{a}}{2} \right\}^{\frac{1}{a}}{} = \sigma_{u} \end{equation} | (32) |
This equation, combined with a proper orthotropic linear transformation tensor, was proposed by Barlat et al. [23] as a yield condition called Yld91. Note that the reference stress is not necessarily the uniaxial stress but a proper identification must be performed. Moreover, if C is the fourth order identity tensor, eq. (32) reduces to the isotropic function in eq. (22). It is worth noting that linear transformations of the stresses do not change the convexity of a function. Therefore, the convexity of the yield condition expressed in eq. (32) is guaranteed. In a similar way, the more general isotropic yield function in eq. (24) with different types of anisotropy and material symmetries were considered in Karafillis and Boyce [30]. The isotropic approach that accounts for a SD effect in eq. (25) was extended to orthotropic anisotropy in Cazacu et al. [53].
In general, any isotropic yield function expressed with the principal stresses can be applied to anisotropic behavior by defining a number of linear transformations such as
| \begin{equation} \tilde{\boldsymbol{\sigma} }^{(k)} = \mathbf{C}^{(k)}{:}\boldsymbol{\sigma} ' \end{equation} | (33) |
with suitable anisotropy coefficient tensors. With two linear transformations, i.e., k = 1, 2, Bron and Besson [54] extended eq. (24) to orthotropic anisotropy, and Barlat et al. [24] proposed the so-called Yld2004-18p yield function
| \begin{equation} \bar{\sigma} = \left\{\frac{1}{4} \sum_{p,q = 1}^{3} | \tilde{\sigma}_{p}^{(1)} - \tilde{\sigma}_{q}^{(2)} |^{a} \right\}^{\frac{1}{a}} {}= \sigma_{u} \end{equation} | (34) |
which generalizes eq. (17). Although 18 coefficients are available in this yield condition, van den Boogaard et al. [55] demonstrated that only 16 are independent. Note that the transformed stresses are not necessarily deviatoric, even though the formulation is independent of the mean stress. Thus, the principal stresses σ1, σ2 and σ3 in eq. (14) must be used to calculate the principal values of the transformed stresses. Aretz and Barlat [56] proposed different isotropic forms for extension to anisotropy. Based on the same principle, Plunkett et al. [57] extended eq. (25) further with two linear transformations.
Specific formulations were also developed for plane stress cases, which are useful in finite elements (FE) for shell formulation in sheet forming. For instance, the so-called Yld2000-2d yield function (Barlat et al. [58]) is also a generalization of Hershey’s isotropic yield function
| \begin{align} \bar{\sigma} &= \{(| \tilde{\sigma}_{1}^{(1)} - \tilde{\sigma}_{2}^{(1)} |^{a}) + (| 2\tilde{\sigma}_{1}^{(2)} + \tilde{\sigma}_{2}^{(2)} |^{a} + | 2\tilde{\sigma}_{2}^{(2)} + \tilde{\sigma}_{1}^{(2)} |^{a})\}^{\frac{1}{a}} \\ &= \sigma_{U} \end{align} | (35) |
The yield function itself is isotropic because, for each transformation, it is possible to switch the two indices ‘1’ and ‘2’. Note that this yield condition was proposed under different forms in the literature (Aretz [59, 60]; Banabic et al. [61]) as pointed out by Barlat et al. [62]. Comsa and Banabic [63] introduced a yield condition expressed as a series in which the terms can be obtained by linear transformations (see also Vrh et al. [64]). Figure 2 shows the application of Yld2000-2d in eq. (35) to the case of a dual-phase steel.
Differential hardening behaviour of a 590 MPa steel sheet subjected to linear stress paths. (a) Measured stress points comprising contours of plastic work. $\varepsilon_{0}^{\text{p}}$ indicates a uniaxial tensile plastic strain in the RD representing a level of plastic work. The solid lines indicate the differential hardening model based on the Yld2000-2d yield function (Barlat et al. [58]). (b) The directions of plastic strain rates β measured at different levels of plastic work corresponding to $\varepsilon_{0}^{\text{p}}$ and those calculated using the associated Yld2000-2d yield function shown in (a) [84].
Since convexity is not affected by a change of reference frame (Lippman [65]) and by linear transformations of the variables such as in eqs. (31) and (33) (Rockafellar [66]), the anisotropic formulations presented in this section preserve the convexity of the associated isotropic function. The yield surface shapes obtained by this type of formulations are consistent with experimental results for steel sheets (Andar et al. [67], Deng et al. [68], Kuwabara et al. [69–71], Kuwabara and Nakajima [72], Kuwabara and Sugawara [73], Verma et al. [74], aluminum alloy sheets (Banabic [75], Yamanaka et al. [76], Kuwabara et al. [77], Yanaga et al. [78]), pure titanium sheets (Ishiki et al. [79], Nixon et al. [80], Sumita and Kuwabara [81]), magnesium alloy sheets (Andar et al. [82]) and a pure zirconium sheet (Plunket et al. [83]). However, the recommended exponents of 8 and 6 for FCC and BCC materials, respectively, are only estimates in the absence of better information. When the exponent is measured as in Hakoyama and Kuwabara [84], Kawaguchi et al. [85], and Yanaga et al. [86], the yield condition is in better agreement with experimental results (see Fig. 2).
d. Other formulationsOther approaches to anisotropy are available as presented in Yu et al. [87], and Skrypeck and Granczarski [88], but only a few are discussed in this section. For planar isotropy, any two-dimensional yield condition writes as a function of two variables
| \begin{equation} x = \frac{\sigma_{1} + \sigma_{2}}{2\sigma_{b}}\quad \text{and}\quad y = \frac{\sigma_{1} - \sigma_{2}}{2\sigma_{s}} \end{equation} | (36) |
Based on this idea, Budiansky [89] proposed a general yield criterion expressed in polar coordinates
| \begin{equation} x = g(\zeta)\cos \zeta\quad \text{and}\quad y = g(\zeta)\sin \zeta \end{equation} | (37) |
where ζ and g(ζ) > 0 are the polar angle and radial coordinate of a point located on the yield surface, respectively. Tourki et al. [90] proposed a yield criterion using a variable substitution similar with that introduced by Budiansky. This formulation is more general because it allows for planar anisotropy cases due to the dependence of g on two parameters ζ and χ. In Tourki et al. [90], the isotropic Drucker [27] criterion was extended to orthotropic symmetry such as a rolled sheet. Tong [91] introduced Fourier series to describe plane stress yield conditions for orthotropic sheets. Vegter et al. [92] proposed to use a Bézier interpolation of experimental data points to define the yield locus. In order to describe planar anisotropy, Vegter’s criterion contains 17 parameters. The most important advantage of this criterion is the flexibility ensured by the large number of parameters. More recently, Peng et al. [93] proposed another type of interpolation with Hermite polynomials. However, these approaches require a large number of experimental data and are restricted for the description of particular stress states only.
Based on the framework of the theory of representation of tensor functions (e.g., Boehler [94], Liu [95]), Cazacu and Barlat [96, 97] proposed generalizations of the stress deviator invariants. These expressions are required to be homogeneous functions of degree two and three in stresses, respectively, reduce to J2 and J3 for isotropic conditions, be insensitive to pressure, and be invariants to any transformation belonging to the symmetry group of the material. In Cazacu and Barlat [96, 97], this approach was used to extend Drucker’s 1949 [27] isotropic yield criterion.
Recent formulations have been proposed mostly to describe the S-D effect either using a pressure-independent approach such as Cazacu and Barlat [96, 97] or a pressure-dependent. These formulations are based on the use of the first or third invariant of the stress tensor, or both, for isotropic and anisotropic materials, and combine with a non-associated flow rule as, for instance, in Lou et al. [98–100] and Wang et al. [101]. Yoon et al. [102] extended Drucker yield condition, eq. (20), for a general stress case, using three tensor invariants with suitable linear transformations. Lou and Yoon [103] extended this model further, without the first invariant to differentiate the behavior of BCC and FCC materials. Chen et al. [104] developed a user-friendly formulation for the same purpose. Based on a plane stress formulation, Hou et al. [105, 106] proposed a model coupling additively two yield conditions, polynomial and invariant-based. A similar coupling is performed for the plastic potential. Other interesting researches were conducted towards the numerical investigation of convexity by Lou et al. [107], or the simplification of non-quadratic yield functions in the case of integer exponents by Cazacu [108, 109]. Additional approaches have been proposed in the last decade but are out the scope of the present review due to its limited length. The reader is referred to the books of Cazacu et al. [110], and Cazacu and Revil-Baudard [111] for further insight on the topic.
e. Strain rate potentialsHill [112] proved the existence of the dual conjugate of a stress potential for rate independent perfect plasticity, i.e., its reciprocal potential expressed in terms of the dual strain rate variables. For instance, Fortunier [113] proposed a formulation for the dual conjugate associated with a single crystal stress potential obeying the Schmid law and Van Houtte [114] with a polycrystal. In general, it is difficult to find the dual conjugate of stress potentials analytically. However, it is possible to define strain rate potentials that describe the material anisotropy as an independent definition of the material behavior. The strain rate potential, proposed by Barlat et al. [115], operates on a linearly transformed plastic rate of deformation tensor $\dot{\tilde{\boldsymbol{\varepsilon}}} = \mathbf{L}\text{:}\dot{\boldsymbol{\varepsilon}}$, and writes as follows
| \begin{align} &\{|\dot{\tilde{\varepsilon}}_{1} - \dot{\tilde{\varepsilon}}_{2} - \dot{\tilde{\varepsilon}}_{3}|^{b} + |\dot{\tilde{\varepsilon}}_{2} - \dot{\tilde{\varepsilon}}_{3} - \dot{\tilde{\varepsilon}}_{1}|^{b} + |\dot{\tilde{\varepsilon}}_{3} - \dot{\tilde{\varepsilon}}_{1} - \dot{\tilde{\varepsilon}}_{2}|^{b}\}^{\frac{1}{b}} \\ &\quad = 3^{a}k\dot{\tilde{\varepsilon}}_{r} \end{align} | (38) |
This potential is viewed as a pseudo-dual of Yld91 and contains six independent coefficients. Numerical applications showed that this potential leads to good approximations of polycrystal potentials using the exponent b = 4/3 and b = 3/2 for FCC and BCC materials, respectively. Kim et al. [116] extended this potential further with two linear transformations. Among more recent developments, it is worth to mention the work of Rabahallah et al. [117] and Cazacu et al. [118]. Stress and strain rate potentials can be used in numerical simulations with equal degree of success (e.g., Li et al. [119]) but stress potentials have received more attention in the literature because they are more convenient for implementation in commercial FE codes. Strain rate potentials are useful for rigid plasticity FE analysis and design codes.
3.3 Anisotropic hardeningThe non-exhaustive review presented above is an attempt to assess the progress made in the field of anisotropic yield conditions over the last few decades. The approach based on linear stress transformation allows a generalization of any isotropic yield function written with principal stresses to anisotropic behavior while preserving the convexity. The above review focuses only on yield conditions that can be applied directly to the flow theory of plasticity assuming isotropic hardening.
However, as already mentioned above, the yield surface evolution is quite complex. For instance, differential hardening effects can be captured using evolving yield function coefficients (Aretz [120], Hakoyama and Kuwabara [84], Peters et al. [121], Yanaga et al. [86]). Generally, though, more accurate models consist in combinations of expansion, translation, distortion and rotation of the yield surface (Kurtyka and Życzkowski [19]). Non-isotropic hardening effects are particularly important for the simulations of forming involving non-linear strain paths as it occurs in many processes [122]. However, it is also important to keep the constitutive model as simple as possible, i.e., optimized, for application in industrial forming problems.
Approaches based on isotropic-kinematic hardening, in which the yield surface translates and expands (e.g., Chaboche [123], Yoshida et al. [124, 125], Lee et al. [126]), are among such optimized models. Another possible approach is a distortional response (Barlat et al. [127, 128]) with no kinematic hardening variable where the effective stress reduces to any isotropic or anisotropic yield function such as described in the sub-sections above, in the absence of distortion.
3.3.1 Kinematic hardeningTypically, kinematic hardening, with a back-stress as a tensorial state variable, has been used to capture the Bauschinger effect. The so-called back-stress β defines the center of the yield surface that translates, thus breaking the center of symmetry in stress space. In the earlier models, the evolution of the back-stress led to linear hardening (Prager [129]; Ziegler [130]) but were hardly applicable to materials such as metals. Later, Frederick and Armstrong [131] and Chaboche [132], proposed non-linear kinematic hardening models in which the yield condition of an initially isotropic material is defined as
| \begin{equation} \bar{\sigma}(\boldsymbol{\sigma} - \boldsymbol{\beta} ) = \sigma_{r}(\bar{\varepsilon}) \end{equation} | (39) |
where $\bar{\sigma }$ and $\bar{\varepsilon }$ are the Mises effective stress and effective strain, respectively and $\sigma_{r}(\bar{\varepsilon })$ is a reference hardening function. In these models, it is customary to write the hardening function as the sum of the yield stress and $R(\bar{\varepsilon })$, the added hardening stress due to plastic deformation, i.e.
| \begin{equation} \sigma_{r}(\bar{\varepsilon}) = \sigma_{y} + R(\bar{\varepsilon}) \end{equation} | (40) |
In Chaboche model [132], the non-linear back-stress evolution equation takes the form
| \begin{equation} d\boldsymbol{\beta} = \frac{2}{3}Cd\boldsymbol{\varepsilon} - \gamma \boldsymbol{\beta} d\bar{\varepsilon} \end{equation} | (41) |
where C and γ are two constant coefficients. Finally, the classical associated flow rule, indicating the normality between the strain rate and the yield surface, and symbolically expressed as
| \begin{equation} d\boldsymbol{\varepsilon} = d\lambda \frac{\partial \bar{\sigma}(\boldsymbol{\sigma} - \boldsymbol{\beta} )}{\partial \boldsymbol{\sigma} } \end{equation} | (42) |
is generally recognized as a good approximation. Equation (42) corresponds to nine equations each for the component of each tensor with a suitable index. As usual, the symmetry of the tensors involved reduces the model to six independent components. The effective strain is calculated in a similar fashion as in isotropic hardening
| \begin{equation} d\bar{\varepsilon} = \frac{(\sigma_{pq} - \beta_{pq})d\varepsilon_{pq}}{\bar{\sigma}} \end{equation} | (43) |
and leads to $d\lambda = d\bar{\varepsilon }$. Kinematic hardening has been mostly developed to capture the material behavior after load reversal, which includes all the phenomena associated with the Bauschinger effect. In particular, this model provides good results for cyclic plasticity. Chaboche [133] provides an excellent review on the topic. Such a kinematic hardening approach was refined by Teodosiu and Hu [134] to account for additional effects that occur during arbitrary strain path changes such as latent hardening. A detailed description of this approach is out of the scope of the present article. This model is particularly well suited to low carbon steels as it was designed based on the dislocation mechanisms observed by transmission electron microscopy in this material (Peeters et al. [135, 136]).
3.3.2 Multi-surface modelsA multi-surface yield condition was proposed by Mroz [137]. Each surface translates inside the surface of immediate larger size and is associated with a given constant hardening modulus. Therefore, any stress-strain curve is approximated by a sequence of linear hardening segments. Two-surface models were proposed by Krieg [138] and Dafalias and Popov [139], and lead to a smooth stress-strain response except, of course, when the strain path changes abruptly. Among the numerous multi-surface kinematic hardening models proposed in the literature, the formulation proposed by Yoshida and Uemori [140], denoted YU, has received much attention in the forming community. Since kinematic hardening has been reviewed in detail in Chaboche [133], the remaining of this section focuses more on distortional plasticity.
3.3.3 Distortional modelsIn order to improve the predictive accuracy, kinematic hardening has been supplemented by yield surface distortion such as in Ortiz and Popov [141], Voyiadjis and Foroozesh [142], François [143] and Feigenbaum and Dafalias [144] to name a few. Kinematic, distortional and other forms of hardening have also been proposed by Kurtyka and Życzkowski [19]. More recently, a distortional plasticity approach was introduced in which the yield surface expands and distorts but does not translate in stress space (Barlat et al. [127]). It is important to note that this model was the first that describes the Bauschinger effect without a back-stress. The most recent version of this model is described in more details in this section. Other distortional models were developed to predict the material behavior for complex loading such as by He et al. [145], Mánik et al. [146] and Holmedal [147]. Qin et al. [148–150] investigated the performance of different constitutive models on the prediction of strain path changes and suggested suitable modifications.
The HAH (Barlat et al. [127, 128, 151]) and HEXAH (Reyne and Barlat [153]) families have no back-stress and are therefore similar to the anisotropic plasticity model with isotropic hardening. The main difference is that, although only one state variable is available for isotropic hardening, $\sigma_{R}(\bar{\varepsilon })$, additional state variables have been introduced in the HAH-HEXAH family. Since the distortion depends on the actual load and its direction, a so-called microstructure deviator denoted $\hat{\mathbf{h}}$ tends to remain aligned with the stress tensor. For monotonic loading, $\hat{\mathbf{h}}$ is always superimposed with the loading direction. $\hat{\mathbf{h}}$ is a deviator with constant size, in particular, $\hat{h}_{pq}\hat{h}_{pq} = 1$ in the later versions of the model. The distortional yield condition for HAH11 (Barlat et al. [127]) is
| \begin{equation} f(\boldsymbol{\sigma} ',f_{+},f_{-},\hat{\mathbf{h}},\bar{\varepsilon}) = \bar{\Sigma} - \sigma_{R} = (\bar{\sigma}^{q} + \phi^{q})^{\frac{1}{q}} - \sigma_{R} = 0 \end{equation} | (44) |
When ϕ = 0, $\bar{\Sigma } = \bar{\sigma }(\boldsymbol{\sigma}') = \sigma_{R}(\bar{\varepsilon })$, the model reduces to the isotropic and anisotropic plasticity formulations presented above for isotropic hardening. ϕ controls the reverse loading distortion and writes as
| \begin{align} &\phi (\boldsymbol{\sigma} ',f_{+},f_{-},\hat{\mathbf{h}}) \\ &\quad = \{f_{+}^{q}| \sigma_{pq}'\hat{h}_{pq} + | \sigma_{pq}'\hat{h}_{pq}||^{q} + f_{-}^{q}| \sigma_{pq}'\hat{h}_{pq} - | \sigma_{pq}'\hat{h}_{pq} ||^{q}\}^{\frac{1}{q}} \end{align} | (45) |
ϕ is nil when f+ = f− = 0. ϕ corresponds to two planes orthogonal to the direction supporting $\hat{\mathbf{h}}$ in stress space that truncate each side of the standard (isotropic hardening) yield surface. The variables f+ and f− are replaced by two other variables, g+ and g−, respectively, that are directly linked to the compression-tension flow stress ratios,
| \begin{equation} f_{\omega} = \left(\frac{1}{g_{\omega}^{q}} - 1 \right)^{\frac{1}{q}},\quad \omega \equiv \text{“+” or “$-$”} \end{equation} | (46) |
ϕq = 0 when g+ = g− = 1. The evolution equations of the variables g+, g−, $\hat{\mathbf{h}}$ and σR are described in Table 1.
;%0A%09%09%09newWindow.document.open();%0A%09%09%09newWindow.document.write('<img src=%22./Graphics/MT-L2024010tab01.jpg%22>');%0A%09%09%09newWindow.document.close();%0A%09%09)
Note that two additional variables g3 and g4 were introduced to capture permanent softening, i.e., when the material remain asymptotically softer after a reversal. Moreover, $d\hat{\mathbf{h}}/d\bar{\varepsilon }$ represents a rotation of $\hat{\mathbf{h}}$ in the deviatoric plane. The initial value of all the variables gk is 1 while the initial value of $\hat{\mathbf{h}}$ is that of the stress deviator when yielding occurs first. It is out of the scope of the present article to describe the other versions of the HAH-HEXAH family in detail. Only short descriptions of the improvement brought by the subsequent models are given. The reader is referred to the original papers for a detailed description of these models.
In a subsequent version of the model, HAH14, Barlat et al. [128] introduced cross-loading effects. It means that the yield surface shape contracts or expands in the cross-loading direction, which is orthogonal to $\hat{\mathbf{h}}$ in the deviatoric stress space. The contraction is similar with a Bauschinger effect in the cross-loading direction and the expansion occurs when the reloading stress overshoot the original stress-strain curve due to latent hardening. This expansion or contraction is produced by a linear transformation of $\bar{\sigma }$ in the appropriate direction. Although this model works for any stress state and any strain path, it exhibits a singularity near the pure cross-loading condition due to the rotation of $\hat{\mathbf{h}}$.
This singularity is almost removed in the next version of the model called HAH20 (Barlat et al. [151]) in which a second microstructure deviator is introduced with the unique purpose of delaying the rotation of the first microstructure deviator when the load changes. Another improvement was the introduction of a pressure effect in spite of a negligible volume change that occurs during the plasticity of metals. This was done for application to advanced high strength steels (AHSS) for which this effect is significant. Although, the goal of HAH20 was achieved, this model contains more variables than its predecessors and a relatively large number of coefficients, although recommended values are provided to the users. As an application example of the HAH model family, Fig. 3 represents uniaxial tension-compression-tension cycles with two different amplitudes for a DP780 steel sheet sample [152]. The HAH20 model leads to an excellent prediction and so does the HAH14 model, although not shown in the figure. Figure 4 shows examples of two-step uniaxial tension tests in different directions for the DP780 steel. In one case, both models lead to excellent predictions while, in the other case, HAH14 underestimate the experiment significantly, which was one of the reasons HAH20 was developed.
Experimental tension-compression-tension cycles with 3% and 6% amplitudes for a DP780 steel sample, and predictions using the HAH20 model [152].
Experimental strain-strain curves for a DP780 steel sheet sample and predictions using HAH14 and HAH20 models. (a) TD uniaxial tension after 6% pre-strain in RD uniaxial tension and; (b) TD uniaxial tension after 6% pre-strain in uniaxial tension at 45° from RD [152].
Although these models, particularly HAH20, have been used successfully in a number of examples, including industrial cases, the introduction of a modified approach called HEXAH was deemed necessary. Since rotations of the microstructure deviator tends to complicate further development of HAH models, the main improvement of HEXAH (Reyne and Barlat [153]) was to introduce a new deviator as soon as the rotation of $\hat{\mathbf{h}}$ is too large. This leads to a suite of deviators $\hat{\mathbf{h}}^{(n)}$, $\hat{\mathbf{h}}^{(n + 1)}$ and $\hat{\mathbf{h}}^{(n + 2)}$ for which only the last is active while all the previous ones become inactive. As a result, the influence of the inactive deviators tends to decrease as the accumulated plastic deformation increases. When this influence becomes negligible, the corresponding deviator is removed. Figure 5 demonstrates the difference in successive yield surface evolution between HEXAH and HAH20. The stress state and the yield surface shape after pre-straining using the HAH20 model corresponds to Shape “1” in Fig. 5. When an abrupt load change follows, the yield surface eventually ends up to Shape “2” after a certain amount of accumulated strain. The angle between Shapes 1 and 2 is slightly larger than π/2. The green lines represent different stages in the evolution of the yield surface between Shape “1” and Shape “2”. When compared with Shape “3” changing to Shape “4” for an angle slightly lower than π/2, it appears that the evolution is somewhat different. This difference would virtually vanish even if Shapes “1” and “2” (or “3” and “4”) were even closer from each other. This is what is called a singularity for this particular loading. However, using HEXAH for the same non-linear strain path, Fig. 5 indicates that this singularity has disappeared because the yield surface evolution is identical.
Yield surface evolution in deviatoric plane at the end of plane strain tension in transverse direction ey (vertical) and plane strain tension in rolling direction ex (horizontal). HAH20 model in top row with yield surface (in position labeled 1) and $\hat{\mathbf{h}}$ rotating yield surface in (to position labeled 2). Same loading change with HEXAH model in bottom row with to microstructure deviators, $\hat{\mathbf{h}}^{(1)}$ vertical and $\hat{\mathbf{h}}^{(2)}$ horizontal [153].
Although accurate finite element implementations are necessary for all the constitutive models, this is all the more important for models such HAH and HEXAH because some state variables may vary drastically upon strain path change and may take slightly inaccurate values that are beyond the range acceptable for these variables. In that case, calculations cannot proceed further and the code generate an error message. Moreover, robustness is an essential aspect because of the complexity of an industrial problem including boundary conditions, contacts, etc. Therefore, Yoon et al. [154–156] developed a number of improvements and analysis tools based on the work by Scherzinger [157]. Moreover, Yoon et al. [156] developed a non-iterative stress integration algorithm allowing a very robust convergence as well as significant gains in computational efficiency. These methods are not described here but are well-documented in the original articles. Moreover, a code including these implementations is available at a website [158]. To demonstrate the robustness of these approaches, the FE simulation of a B-pillar with complex geometry was conducted for a high strength TRIP-1180 steel sheet using the pressure sensitive version of the HAH20 model. In Fig. 6, the contours of the stress as well as those of an index describing the load history in the last few percent strains are represented.
FE forming simulation results of a B-pillar for a TRIP1180 steel blank showing contours of: (a) stress levels achieved in the part and; (b) index characterizing the severity of the strain path change. This example demonstrates the robustness of the HAH20 model and its FE implementation [156].
This section reviews the material test methods that can reproduce the various stress states occurring in real forming operations. These test methods are essential to verify the validity of material models. Examples of test data and material modeling results for AAs are also presented.
4.1 Material test methods using linear loadingsIn this section, the typical plane stress states that form a yield surface in the σxx − σyy − σxy stress space and the common material test methods for measuring every section of the yield surface (Fig. 1) are reviewed.
In Fig. 1, Curve A is a set of yield stresses measured using uniaxial tensile tests in different directions, defined by their angle θUT from the RD. The uniaxial tensile test method is standardized as ISO 6892-1 [159]. As a uniaxial tensile test method for circular tube materials, Dick and Korkolis [160] investigated the mechanics of a Ring Hoop Tension Test using 3D Digital Image Correlation for the strain measurement and applied it to assess the properties of a AA6061-T4 tube in the hoop direction.
Curve B in Fig. 1 defines stress states in the first quadrant of the yield surface (σxx ≥ 0; σyy ≥ 0). It can be measured using cruciform specimens or cylindrical tubular specimens loaded with axial force and internal pressure. Various cruciform specimen shapes have been proposed in the literature [12, 161]. Figure 7(a) shows the cruciform specimen specified in ISO16842 [162]. This specimen was first proposed by Kuwabara et al. [69]. According to finite element analyses, the stress measurement errors associated with the use of the cruciform specimen is less than 2% [163, 164]. The maximum equivalent plastic strain, $\varepsilon_{\text{max}}^{\text{p}}$, that can be measured with the cruciform specimen ranges between 0.002 and 0.05. It is affected by the arm slit width, the work hardening exponent of the test material and the shape of the yield surface [16]. As a way to extend $\varepsilon_{\text{max}}^{\text{p}}$ to higher levels, a method of decreasing the thickness of the gauge area where the stress is measured, see Fig. 7(b) [68, 165–167], and a method of increasing the strength of the arms [168, 169] have been proposed.
Using the biaxial tensile test methods with cruciform specimens, the following studies have identified the proper yield functions of pure aluminum and AA sheets: 1145 [170], 3104-H19 and 5182-O [171], 5754 [167], 5754-O [172], 6016-T4 [78, 86], 5182-O [85], 5XXX [173], 6014-T4 [174], 6016-O and -T4 [175], 6016-O and -T4P [176], and 5052-O [177].
Yamanaka et al. [178] proposed deep neural network (DNN) approaches to efficiently estimate biaxial stress-strain curves of aluminum alloy sheets from a digital image representing the sample’s crystallographic texture. It was observed that the DNNs could estimate biaxial stress-strain curves with high accuracy.
Several cruciform specimens with reduced thickness in the gauge area were proposed to measure the forming limit and/or fracture limit diagrams of 5000-series AA sheets [179–183]. A method for measuring the forming limits of a 6061-O sheet at fracture from shear to plane strain tension region was developed [184]. A new biaxial testing system for forming limit diagram measurements of sheet metals under hot stamping conditions (400–500°C) was developed [185].
Circular tube test pieces with uniform wall thickness have been used for a long time to measure the yield surface by applying an axial force and an internal pressure [14, 15, 186–192]. This biaxial stress test method is hereinafter referred to as the “Biaxial Tube Bulge Test (BTBT) method”.
Using the BTBT, the material models and/or forming limit strains/stresses under proportional and non-proportional loadings were measured for circular AA tubes: 5154-H112 [193, 194], 6260-T4 [195–198], 6061-O [199], and AA2219-T852 [200].
Yoshida et al. [194] experimentally verified that the forming limit stresses of an extruded AA5154-H112 tube are path-independent. Yoshida and Suzuki [201] analyzed the forming limit stresses of sheet metals subjected to linear and combined stress paths using the M-K model in conjunction with two anisotropic work-hardening models: a work-hardening model capable of describing Bauschinger and cross-hardening effect, and another model that cannot predict the cross-hardening effect. It was found that the forming limit stress is path-independent when the stress–strain curves for the linear and combined stress paths agree well with each other. On the other hand, the forming limit stress for the combined stress path depends on the strain path when the prestrain changes the subsequent stress–strain relation.
Tiji et al. [202] developed a servo-controlled experimental setup for BTBT. Using a FEM model with PID control and simple material model such as von Mises yield condition, the boundary conditions to create different deformation modes in BTBT were generated. Considering the equivalency of the plastic work, the evolution of the yield surface and strain-rate potential of an AA7075-O extruded tube were characterized using the Yld2004-18p yield function [24].
The cryogenic deformation mechanism of an AA 6061 tube at −196°C was revealed by microstructure characterization. It was also found that the maximum expansion strain of the tube was 34.0% ± 0.6%, that is, 99.8% higher than that at room temperature [203].
An experimental method for generating combined tension and shear stress-states in thin-walled tubes was proposed to calibrate anisotropic yield functions of an AA 6061-T4 tube [204].
The BTBT can be used as a biaxial stress test method for sheet metals uniformly bent and welded to fabricate a circular tube test piece that is subjected to axial force and internal pressure [73, 79, 84]. Then, an arbitrary biaxial tensile stress state is applied to the mid-section of the tube. The RD of the sheet sample may be either the circumferential or the axial direction of the tubular test piece. Methods for measuring the strain include strain gauges [79], strain measuring devices using displacement meters [73], and digital image correlation (DIC) [84].
Circular tube specimens are work-hardened by bending. For this reason, the stress-strain curves measured by the BTBT are generally higher than those measured by the cruciform specimen. In order to correct the effect of the work hardening due to bending, a method has been proposed in which the stress-strain curve obtained from the BTBT is translated in a direction parallel to the strain axis and smoothly connected to the stress-strain curve measured using the cruciform specimen [73].
Employing the cruciform specimen test and the BTBT, the DH of AA sheets 5182-O [85], and 6016-T4 and -O [175] were modeled using the Yld2000-2d yield function [58]. Figure 8 [175] shows the stress points that form the contours of plastic work for 6016-O and -T4. The maximum value of $\varepsilon_{0}^{\text{p}}$, for which the work contour has a full set of stress points measured using cruciform specimens, was $\varepsilon_{0}^{\text{p}} = 0.04$ for 6016-O and 0.06 for 6016-T4. The stress points for a larger strain range were obtained using the BTBT. The stress data in balanced biaxial tension (σx:σy = 1:1) were measured using the hydraulic bulge test (HBT) method for $\varepsilon_{0}^{\text{p}} \geq 0.05$ (6016-O) and for $\varepsilon_{0}^{\text{p}} \geq 0.07$ (6016-T4).
The point C in Fig. 1 corresponds to the stress state where the yield locus B protrudes most in the σx direction, where the material is subjected to plane-strain tension in the x-axis: the y-component of Dp is 0. Several authors proposed specimen geometries suitable for plane-strain tension tests [205–207]. Gille et al. [208] proposed a specimen shape suitable for measuring the forming limit strain in plane strain tension of 5000-series and 6000-series AA sheets and verified its performance by measuring the strain distribution using DIC. It should be noted that the specimen geometry is different from that suitable for the plane-strain tension of high-strength steel sheets [209].
The stress point D in Fig. 1 corresponds to the balanced biaxial tensile yield stress, σb. As a test method for measuring σb, either a balanced biaxial tensile test using a cruciform specimen [162] or a HBT is a suitable choice. ISO standard has been established for the HBT using an optical strain measuring device [210]. Failure is delayed in the HBT, allowing the measurement of the material response at significantly larger strains than in the uniaxial tension test. If the material is anisotropic, the stress state at the apex deviates from equibiaxial tension [211]. A more accurate HBT that takes into account elastic strain, the effect of bending, the thickness ratio to the apex radius and the difference in the radius of curvature between RD and TD, has been proposed [212]. It was shown that the flow curve could be measured with a very high accuracy by fitting the surface coordinates to an ellipsoid shape function and by considering the local strain data to approximate the curvature of the midplane [213]. Experimental data obtained for bulge diameters as small as 50 mm and with die opening diameter to sheet thickness ratio ranging from 42 to 150 were compared with the results calculated using a hardening law identified with an algorithm specific to small bulge diameters [214]. A methodology of determining the anisotropy material parameters of the Yld2004-18p yield function [24] was developed using uniaxial tension, plane strain tension and HBT data, and the deformation and bursting behavior of an AA 2024 sheet was successfully predicted [215].
If the yield stress of a material is independent of the hydrostatic pressure, a classical assumption in the plasticity literature, the balanced biaxial yield stress, σb, is identical to the uniaxial compressive yield stress in the thickness direction [216–218]. In the compression test of a stacked sheet specimen, the elastic deformation and friction coefficient associated with the compression dies are greatly reduced if the latter is composed of a cemented carbide surface coated with diamond-like carbon, thereby improving the measurement accuracy of the compressive yield stress in the thickness direction, i.e., $\sigma_{z}^{\text{C}}$ [219, 220].
Curve E in Fig. 1 corresponds to the yield stresses of the material subjected to biaxial tension with the principal axes inclined by 45° from the RD.
Curve F in Fig. 1 is obtained by performing a combined tension-compression test with a ratio of the principal stresses of σ1:σ3 = 1:−1 [221], or a simple shear test in various directions. Figure 9 shows the specimen geometries proposed for a simple shear test [222–231]. Yin et al. [232] concluded that the shear stress-strain curves are almost identical regardless of specimen shape when an optical strain measuring device is employed.
Material test methods for identifying a yield surface by the combined use of uniaxial, plane-strain, and biaxial tensile tests were proposed [233–237]. Takizawa and Kodama [236] and Saito and Takizawa [237] determined a polygon that passes through the yield point in uniaxial, equibiaxial and plane strain tension tests, and proposed a method to determine the work contour using the curve inscribed in the polygon as shown in Fig. 10. Zhang et al. [238] performed the biaxial tension test of an AA5086 sheet to identify the yield surface by comparing the experimental and numerical principal strains along a diagonal direction of the specimen central area. They concluded that a single biaxial tensile test seems sufficient to obtain all the material parameters of a complex yield criterion for the test sample.
Schematic illustration for determining a yield locus as a polygon circumscribing stress points forming a work contour [237]. (a) Stress-strain curves in uniaxial, plane strain and equibiaxial tensile tests. (b) Identified yield loci using the Yld2000-2d [58] and Yld2004-18p yield functions [24]. Material: A5052-O and -H32.
Curves G in Fig. 1 are obtained from combined plane strain tension-shear tests, see Fig. 11 [239–245]. In this test method, the normal stress on a plane parallel to the tensile axis, i.e., σyy in Fig. 11(c), cannot be measured. Therefore, it is often assumed equal to half of the normal stress from the analogy of plane-strain tension for an isotropic material.
The stress point H in Fig. 1 is an in-plane uniaxial compressive yield stress, which is measured to evaluate the Bausinger effect (BE) during stress reversal or the asymmetry of the yield stress between tension and compression (Tension-Compression Asymmetry (TCA); TCA is often referred to as Strength Differential Effect (SDE)). There are two types of uniaxial compression test methods for sheet metals [12]. One is a method in which a thin sheet specimen is subjected to in-plane compression while being sandwiched between rigid blank-holders that prevent buckling of the specimen. The other one is a method in which a stacked specimen is subjected to uniaxial compression [246].
The following studies have been conducted on material models that reproduce the BE of AA sheets. Qin et al. [149] proposed a new model for predicting stress transients caused by strain path changes. The microstructure stress deviator is used to memorize and model the evolution history of the microstructure. The model captured well earlier published experiments for a commercially pure aluminum, an extra deep drawing quality steel and a dual-phase steel. Barrett et al. [247] investigated the continuous-bending-under-tension (CBT) on an AA6022-T4 sheet. Cyclic tension-compression experiments were performed on this material to determine a non-linear kinematic hardening model. The model replicated the development of strain on the surface during CBT, and compared well with post-test strain measurements. Daroju et al. [248] applied an elasto-plastic self-consistent (EPSC) crystal plasticity model to predict and interpret reverse loading in two precipitation-hardened AAs, 6016-T4 and 7021-T79. Comparison of the experimental and modeling results revealed that the unloading behavior is primarily driven by the backstress, the BE by the backstress and inter-granular stresses, and the hardening rates upon load reversals primarily by the strain-path sensitive evolution of the dislocation density. Figure 12 shows examples of the measured and calculated stress-strain curves of an AA 7021-T79 sheet under loads reversals.
Comparison of simulated and measured true uniaxial stress-true strain response in strain path reversal deformation for AA 7021-T79 sheet along TD [248].
Richmond and Spitzig [33, 249] performed uniaxial tension and compression tests on various high-strength steels, Fe-based single crystals, and 1100 commercially pure aluminum under the superposed action of a range of hydrostatic pressures. They found that the flow stress σ increases proportionally with the hydrostatic pressure p and expressed this relationship by the following equation
| \begin{equation} \sigma = \sigma_{0}(1 + 3\alpha p) \end{equation} | (48) |
where σ0 is the yield stress at p = 0. These authors measured the so-called pressure coefficient α and found it constant in a range of 13 to 23 TPa−1 for ferrous materials, and 56 TPa−1 for the 1100 pure aluminum.
Several authors measured the TCA in AAs. Barlat et al. [250] conducted tension and compression tests on an AA sheet, 2090-T3, at increments of 15 degrees from the RD to measure the respective 0.2% proof stresses, σT and σC. They found that σT was 12% larger than σC in the RD and that σC was several percent larger than σT in the directions at 45°, 60°, 75°, and 90° from the RD, confirming the TCA. Kuwabara et al. [251] developed comb-shaped dies for performing the in-plane compression test of a sheet metal without causing buckling. They found that two mild steel sheets exhibited the TCA and that an AA sheet 5182-O did not. Fourmeau et al. [252, 253] observed the TCA in a 20 mm-thick AA 7079-T651 plate. They found that the compressive stress was higher than the tensile stress in the diagonal direction (45° from RD) but did not observe the TCA in the RD for $\varepsilon_{x}^{\text{p}} \leq 0.05$ [252]. The compressive stress was higher than the tensile stress in the normal direction (ND) whereas the tensile stress was higher than the compressive stress in the RD for a strain range of $\varepsilon_{x}^{\text{p}} >0.05$ [253]. Kabirian et al. [254] found no TCA in an AA5182-O sheet. Seidt and Gilat [255] measured the TCA in a 12.7-mm-thick 2024-T351 aluminum plate at a strain rate of 1 s−1 and found that the flow stresses in compression were slightly higher (∼4%) than those in tension. Holmen et al. [256] investigated the TCA in four 6xxx series AAs in several tempers with yield strengths varying from 27 to 373 MPa. The TCA was as large as 6.3% but it was not observed in the annealed materials. Ku et al. [257] measured the thermo-mechanical behavior and texture evolution of two overaged Al 7056 alloy plates, in T761 and T721 tempers, over a wide range of strain rates (10−4–3 × 10−3 s−1) and temperatures (22–300) °C under uniaxial tension and compression along the ND. They concluded that the difference in the Schmid factor between tension and compression contributed to the TCA. Barros et al. [258] performed the numerical simulation of the cylindrical cup drawing of a 2090-T3 AA considering the TCA. They recommended to use a yield criterion that is flexible enough to describe the anisotropy in both the yield stresses and r-values, including those in compression.
Figure 13(a) [259] shows the TCA observed in an AA 5083-O sheet with a higher flow stress in compression than in tension. Figure 13(b) compares the $\sigma_{\text{b}} - |\varepsilon_{z}^{\text{p}}|$ curve measured using the balanced biaxial tensile test with a cruciform specimen [162] and the uniaxial compression curve in the thickness direction, $|\sigma_{z}^{\text{C}}|$-$|\varepsilon_{z}^{\text{p}}|$. It is clearly observed that $|\sigma_{z}^{\text{C}}| \approx \sigma_{\text{b}}$ for $|\varepsilon_{z}^{\text{p}}| < 0.08$ indicating the absence of a pressure effect on the flow stress. Thus, the cause of the TCA is thought to be due to factors other than pressure.
Comparison of true stress-true strain curves of AA5083-O sheet between (a) uniaxial tension and compression in RD and (b) equibiaxial tension, σb-$|\varepsilon_{z}^{\text{p}}|$, measured using a cruciform specimen [162], and uniaxial compression in the thickness direction, $|\sigma_{z}^{\text{C}}|$-$|\varepsilon_{z}^{\text{p}}|$ [259].
The yield locus I in Fig. 1 is located in the third quadrant of stress space (σxx ≤ 0 and σyy ≥ 0) and can be measured by biaxially compressing a stacked specimen [260], a piece of square sheet [261], or a cruciform specimen with short arms [262]. The biaxial compression test using a stacked specimen is effective to measure a yield surface in the π-plane [260]. A cruciform test piece having short arms was proposed to measure the BE of sheet metals subjected to biaxial tension-compression stress reversals [262]. The test method is effective in reproducing the stress history of automotive outer panels in a dent resistance test; it is subjected to biaxial tension during press forming followed by biaxial compression during a dent test.
4.2 Material test methods using nonlinear loadings 4.2.1 Loading-unloading-reloading methodLoading-unloading-reloading (LUR) test methods have long been used to measure the subsequent yield surface of prestrained materials [4–12].
Barlat et al. [263] performed non-linear deformation path experiments using uniaxial tension followed by simple shear tests for an AA1050-O sheet sample in different specimen orientations with respect to the material symmetry axes. The relative contributions of the crystallographic texture and dislocation microstructure evolution to the anisotropic hardening behavior of the material were discussed. Figure 14 [263] shows the measured stress-strain curves and TEM micrographs of typical dislocation structures after deformation. The authors concluded that for monotonic loading, or at strains far from a strain path change, the differences observed in the hardening rate in different testing orientations appear to be controlled mainly by crystallographic texture evolution, but that for sequences of two linear strain paths, the anisotropic hardening behavior just after reloading appears to be strongly dependent on the dislocation microstructure evolution.
(a1), (b1), and (c1) Stress vs. plastic work curves obtained in TD uniaxial tension at strains of 0.0, 0.07 and 0.14 followed by 45° simple shear. (a2), (b2), and (c2) The TEM micrographs of typical dislocation structures developed after plastic deformation in TD uniaxial tension up to a strain of 0.14, followed by simple shear up to a strain of 0.15 in different directions; (a2) 45° shear illustrating cell reorganization; (b2) 90° shear showing cell superimposition; (c2) 135° shear showing cell dissolution [263].
The shape of the yield surface changes as plastic deformation progresses. The yield surface measured after the application of a plastic prestrain is called subsequent yield surface. Regarding the possibility that there is a corner on the subsequent yield surface at the loading point, as predicted by crystal plasticity [264], it has been argued by Hecker [5] that if a corner exists, it will be erased by the elastic unloading needed to probe the yield surface. But whether or not a vertex has formed on the yield surface is very important for the predictions of plastic instabilities (Stören and Rice [265], Hutchinson and Tvergaard [266], and Needleman and Tvergaard [267]).
Kuroda and Tvergaard [268] proposed a method for determining the subsequent yield surface in the vicinity of a current loading point without unloading. The procedure, (i) to (iii), is as follows.
This test method does not include unloading. Therefore, if a vertex exists at the point of loading, the stress trajectory before and after the abrupt change in the strain path should naturally draw a vertex. Figure 15 [269] shows the stress trajectories measured by applying the abrupt strain path change method to a 6000-series AA sheet. The total strain rate ratio (RD:TD) during the first loading was set to D11:D22 = 1:1 and, at the moment the nominal strain component reached e11 = e22 = 0.01, the total strain rate ratio D22/D11 was abruptly changed to the values shown in the figure. Just after the abrupt strain path change, the stress trajectory depicted a clear vertex. The stress path calculated using the generalized Hooke’s law from the total strain rate ratio after the path change is completely different from the measured stress path; therefore, it is clear that the stress path after the abrupt strain path change is not elastic unloading. The measured Dp is inclined in the direction of the stress path, and this phenomenon is also predicted by polycrystalline plasticity analysis [268].
Stress paths for abrupt strain path change following (a) equibiaxial tension and (b) plane strain tension. Material: 6000-series AA sheet. The first and second markers on the stress paths (◇) indicate the stress point at which the accumulated equivalent plastic strain Δεp is 0.001 and 0.002, respectively (Δεp is measured from the strain path change point). Directions of the Dp for the second step of straining are shown by lines on the stress paths after strain path change. Open markers in (a) indicate the subsequent yield surfaces defined as contours of equal plastic work measured by LUR method: the levels of plastic deformation are $\varepsilon_{0}^{\text{p}}$ = □ 0.0001 ○ 0.001, △ 0.005, where $\varepsilon_{0}^{\text{p}}$ is the logarithmic plastic strain attained in the uniaxial tension in the RD of the material. -●-: Work contour measured using linear radial stress paths [269].
Yoshida and Tsuchimoto [242] investigated the elastoplastic deformation behavior of thin-walled tubes (30 mm diameter and 1.2 mm wall thickness) made of pure aluminum and steel subjected to linear loading (LL) and nonlinear loading (NL), Fig. 16. In the LL experiments, the ratio between the displacement and rotation of the grip was held constant during loading (Fig. 16(a), left). In the NL experiments, a specimen was subjected to uniaxial tension followed by the simultaneous application of tension and torsion (Fig. 16(a), right). For the LLs, the associated flow rule predicted the plastic flow behavior observed in the experiments with sufficient accuracy, see Fig. 16(b). For the NLs, the plastic flow behavior was much different from that for the LLs; the direction of Dp was rotated toward the direction of the stress/strain rate, see Fig. 16(b). Moreover, the authors found that the difference in the direction of Dp between the LL and NL under the same stress state is roughly proportional to the difference in the total strain rate directions. They proposed a phenomenological model that reproduces the non-normality.
Stress paths and plastic flow of pure aluminum specimen: (a) schematic illustration of displacement of a grip for linear and nonlinear loadings, (b) stress paths for the linear and nonlinear loadings with φNL = 90°. The results for the nonlinear loadings are represented by red lines. “L” and “NL” denote the quantities measured in the linear and nonlinear loadings, respectively [242].
Hama et al. [270] investigated the plastic deformation behavior under various strain path changes in an A6022-T4 Al alloy sheet, focusing on the evolution of the direction of the plastic strain rate, β, and the stress ratio, φ, after abrupt strain path changes. In the cruciform specimen experiments, before the strain path change, the deformation was well represented by the associated flow rule with the Yld2000-2d yield function [58] and the isotropic hardening assumption. After the abrupt strain path change, the deformation temporarily deviated from that predicted by the isotropic hardening model with associated flow rule. However, the isotropic hardening model recovered after the plastic work increased by roughly 4.0 MJ·m−3, which corresponds to a strain increment of approximately 0.024 under uniaxial tension in the rolling direction.
This section introduces the identification of complex material models for AA sheets using biaxial stress test methods, and the verification of the prediction accuracy of the forming simulations using finite element analyses (FEA) based on these models.
5.1 HBTFigure 17 [85] shows the work contours of an AA5182-O measured using the cruciform test [162] for a strain range of $\varepsilon_{0}^{\text{p}} \leq 0.05$ and the BTBT for the strain range of $0.05 < \varepsilon_{0}^{\text{p}} \leq 0.13$, showing that the test sample clearly exhibits DH. The Yld2000-2d yield function [58] accurately reproduces the change in the shape of the work contours (Fig. 17(a)), while those of von Mises and Hill’s quadratic lead to large deviations from the measurement (Fig. 17(b)). Figure 18 [85] compares the experimental and predicted HBT result, $|\varepsilon_{z}^{\text{p}}|$ at apex vs. pressure diagram. The Yld2000-2d (IH) with isotropic hardening is applied with an exponent of 8.86 to accurately approximate the work contour for $\varepsilon_{0}^{\text{p}} = 0.05$. Moreover, two types of differential hardening (DH) models based on the Yld2000-2d yield function were used in the FEA. In the DH-A model, the stress-strain curve in the RD was approximated using Swift’s hardening law in the strain range, $0.002 \leq \varepsilon_{x}^{p} \leq 0.18$. In the DH-B model, the $\sigma_{\text{b}} - \varepsilon_{z}^{\text{p}}$ stress-strain curve was measured using the HBT, converted into a $\sigma_{x} - \varepsilon_{x}^{\text{p}}$ and fitted with a work-hardening law in the strain range, $0.002 \leq \varepsilon_{x}^{p} \leq 0.43$. The calculated curve using the DH-B model is almost identical to the measurement in the whole strain range of $\varepsilon_{z}^{\text{p}}$. The application of DH model is also effective in improving the FEA accuracy of the HBT [86] and the hole expansion test [175] for 6000-series AA sheets.
Biaxial stress test results of AA sheet 5182-O. (a) Stress points forming contours of plastic work and the yield surfaces based on the Yld2000-2d yield function [58]. (b) The stress values are normalized by the σ0 associated with each work contour and the yield surfaces based on the von Mises and Hill’s quadratic yield functions [85].
Hydraulic bulge test result: $|\varepsilon_{z}^{\text{p}}|$ at apex vs. pressure diagram [85].
Extensive benchmark analyses on the cup drawing of an AA 6016-T4 sheet were performed by 11 research teams relying on phenomenological or crystal plasticity approaches, using commercial or self-developed finite element (FE) codes [271]. The ability of the constitutive law and of the FE model to predict r-value and yield stress in different directions was verified. Then, the simulations of the earing profile (cup peak number, average height and amplitude), the punch force evolution and thickness in the cup wall were evaluated and analyzed.
Saito and Takizawa [237] performed a cup drawing FEA for AA5052-O and -H32 sheets using the Yld2000-2d [58] and Yld2004-18p [24] yield functions determined as shown in Fig. 10(b) and compared the calculated results with the experiment, see Fig. 19. The overall tendency of the earing behavior calculated using the Yld2008-18p were consistent with the experiment, however, the calculated cup heights were 3–5% lower than the experiment.
Comparison between experimental and calculated wall height distributions in circumferential direction [237].
Kuwabara et al. [175] performed the hole expansion experiment and FEA using AA6016-O and -T4 and investigates the influence of heat treatment on the anisotropic plastic deformation behavior of the test samples. Although the materials exhibited the DH effect (see Fig. 8), the Yld2000-2d yield function provided proper material representations of the plastic behavior of both materials because it correctly predicted the fracture or localized neck locations that occurred in the hole edge vicinity.
Ha et al. [272] performed a FEA of the hole expansion of AA6022-T4. The orthotropy of the material was modeled with the Yld2000-2d yield criterion [58] and calibrated using uniaxial tension, plane-strain tension and disk compression experimental results. The greatest thinning occurred at ±45° from the RD and the Yld2000-2d yield criterion [58] led to a well predicted thinning variation at the hole edge.
Iizuka et al. [273] performed a hole expansion experiment using a AA6116-T4 sheet and found that the fracture in the hole expansion always occurred, regardless of the hole diameter, when the circumferential strain at the hole edge reached the fracture limit strain (FrLS) identified in the uniaxial tensile test. Therefore, the authors concluded that the FrLS criterion is suitable for predicting the hole expansion limit.
With the emergence of data-driven approaches or artificial intelligence (AI), it is legitimate to ask whether constitutive models will still play a role in the future. Data-driven or AI methods are extremely useful for the optimization of known manufacturing operations because a huge amount of data is available. However, in the design stage of a process, none of the process parameters are imposed and no solid database is available. An alternative option is to conduct forming simulations using data-driven approaches to describe the constitutive relationships [274]. However, experimental data for material models can be generated only for a limited number of stress states unless full field methods [275], discussed in more details below, can extend the number of probed stress states significantly enough. Nevertheless, as constitutive models are becoming more complex, for instance by accommodating strain path change effects, the generation of experimental data to completely describe a full and robust constitutive response is likely to remain out of reach. Accurate constitutive models, which are more powerful for interpolation and extrapolation purposes, are likely to remain indispensable in direct numerical simulations of forming processes, and in the generation of data necessary for an acceptable numerical constitutive description. For the latter, it is yet to demonstrate that the data-driven approach will be efficient for application in forming process modeling of complex components. This, of course, justifies the need of improving the material elasto-plasticity models further.
Regarding the formulation of new generations of constitutive models, the dilemma will be to develop accurate descriptions that are not too complex for application in industry. This includes a robust identification method with a limited number of experimental datasets. This is discussed with further details below. This will also require accurate, robust and time efficient stress integration algorithms. Another related issue is the understanding of the accuracy required for a given forming process. Very sophisticated models might require large costs of development and utilization but lead to insignificant enhanced accuracy in forming process predictions. Therefore, a cost-benefit comparison will be necessary.
As stated in this paper, uniaxial tensile tests alone are not sufficient to verify the accuracy of material models; multiaxial stress tests are essential. However, compared with the uniaxial tensile test method, the multiaxial stress test approach is not as commonly used. The multiaxial stress test methods for which international standards have already been established, include the HBT method [210] and the biaxial tensile test method using a cruciform test piece [162]. In the future, international standards should be established for other mechanical tests, in particular for in-plane reversal test under uniaxial and biaxial stress states, simple shear test, combined plane strain tension-shear stress test and biaxial tube bulge test.
An alternative test approach to the conventional method, which is often referred to as Material testing 2.0 (MT2.0) (Pierron and Grédiac [275], Rossi et al. [276], and Andrade-Campos et al. [277]), consists in identifying material model coefficients based on heterogeneous strain field data collected on a non-uniform deforming specimen. Specifically, the deformation of the specimen is analyzed using a FEA, and the coefficients of the material model assumed in the FEA are optimized so that the errors between the strain field measured using DIC and those calculated using the FEA are minimized. Although, a model must be assumed a priori in the FEA analysis which, for instance, makes the identification of DH model [278] difficult, it has the advantage of requiring fewer material tests than those required with conventional mechanical tests. Further development of the research on MT2.0 is expected.
The authors are thankful for the financial supports from the Light Metal Educational Foundation Inc. of Japan and the Global Innovation Research Organization in TUAT.
