Phenomenological Modeling of Deformation-Induced Anisotropic Hardening Behaviors: A Review

Abstract

Numerous studies indicate that the hardening behaviors of materials are closely related to their deformation history. In the forming processes with loading path changes, such as sheet metal forming, anisotropic hardening behaviors are universally observed. In this situation, selecting or constructing a suitable anisotropic hardening model is essential. This paper presents a review of the phenomenological modeling of the deformation-induced anisotropic hardening behaviors. At the beginning, the deformation-induced hardening behaviors are introduced together with the relevant experiments. Different from other published review works, this paper is not laid out according to the description of a series of models. Instead, the modeling is emphasized by generalizing the main mathematical modeling ideas among various hardening models and sorting out the description methods for the decomposed anisotropic hardening behaviors. Some prospective development directions for the modeling of anisotropic hardening behaviors are suggested at the end of this work. This review work tries to provide the researchers with an instruction on how modeling for the anisotropic hardening behaviors according to the materials and forming processes.

1. Introduction

Finite element simulation of the sheet metal forming has been widely used in industry. To achieve a realistic simulation, the use of improved constitutive models for the reproduction of material’s plastic behaviors, is essential. In constitutive models, the plastic behaviors are reflected by the description of the yield surface, including its initial shape and its evolution during the plastic deformation.

The initial shape of the yield surface is influenced by the crystal structure of the material. The most widely used Mises yield function provides a concise framework under the assumption of an ideal isotropic yield surface. Whereas, as for most actual metal materials, especially the FCC, HCP, or porous materials, more potential features of the yield surface, such as non-quadratic, tension-compression asymmetry, hydrostatic pressure sensitivity, are necessary to be reproduced. Accordingly, other two basic modeling frameworks are brought forward by Hershey [1] and Drucker [2].

Plastic deformation induces the evolution of the yield surface. In the metal forming industry, the rolling acts as a universal process to produce the sheet metal. The pre-rolling of the sheet metal might induce the orthotropic yield or anisotropic yield, which is generally incorporated into the description of the initial yield surface. Among the modeling for the pre-rolling induced initial yield behavior, Hill 48 [3], Barlat 89 [4], and Cazacu 01 [5] could be viewed as the representatives. In regards of the modeling of the initial yield surface, focusing on the influences of crystal structure and pre-rolling, an excellent review has been carried out by Barlat et al. [6] and will not be involved much in this paper. For another, in the forming processes, the yield surface evolves with the plastic deformation, including the changes in its size, central position, and shape. These changes could be associated with various hardening behaviors. The hardening models are often specialized to reproduce the hardening behaviors by describing the evolution of the yield surface.

Hardening behaviors are strain path dependent. If the strain path is monotonic, a material will undergo continual hardening. The hardening behaviors in all loading directions can be approximated by using a unified stress–strain curve, which is defined as isotropic hardening. A typical consequential behavior of the anisotropic hardening is that the subsequent stress–strain curves deviate from the one measured in the proportional loading.

In most of the sheet metal forming processes, the materials are subject to multiaxial loading, where the loading path changes are very likely to occur. In the deep-drawing forming process, for instance, the stress states of the materials at the flange edge gradually evolves from hoop compression to radial tension [7]. Furthermore, in the multi-step forming processes, the loading path changes are almost inevitable. In these conditions, it is no more appropriate to use the isotropic hardening models. For example, in the twist springback simulation of a dual-phase steel conducted by Liao [8], the isotropic hardening model always over-predicts the twist angle, because it omits the anisotropic hardening behaviors after the loading path changes. By adopting the anisotropic hardening models, the accuracy is improved. Not only the springback [8,9], but the prediction of residual stress [10,11] and the failure [12,13,14] is also sensitive to the description of the hardening behaviors. Thus, to ensure the reliability of the sheet metal forming simulation, the modeling of deformation-induced anisotropic hardening behaviors is of vital importance.

To describe the anisotropic hardening behaviors, the hardening models have been constructed either physically or phenomenologically. Since the hardening behaviors are basically dislocation related, physical models are constructed to predict the materials’ macroscopic hardening behaviors based on the evolution of dislocation density. Among them, the Taylor model [15] is one of the classic physical hardening models and has been universally applied in the crystal plasticity simulation. Rauch [16,17] proposed a series of physical hardening models specialized for the deformation-induced hardening behaviors. For more details about the physical mechanisms of the plastic hardening behaviors and the crystallographic models, readers can refer to the relevant chapters in the handbook edited by Lemaitre [18]. Despite these physical modeling methods are instrumental for the users to thoroughly understand the plastic deformation, they are hardly used to deal with the engineering problems due to their complexity and high computational cost. Phenomenological hardening models do not contain the materials’ physical quantities (e.g., the dislocation density) in their expressions. Instead, they aim to reproduce the hardening behaviors in a pure mathematical manner. Owing to their simplicity, time-efficiency, and the accuracy, the phenomenological hardening models cater to the needs of industry [19].

The state-of-the-art of the phenomenological modeling for the deformation-induced anisotropic hardening behaviors is reviewed in this paper. It tries to generalize the main ideas of the modeling among various hardening models and sort out the different description methods of each decomposed anisotropic hardening behavior. For more details of each hardening model, readers are advised to refer to the original text. This paper aims at providing an instruction for the readers to select or construct a proper hardening model according to their materials and forming processes. For format uniformity, the model formulations cited in this paper are not strictly consistent with their original versions. Some symbols are replaced, some coefficients with less importance are omitted, and some formats of the expressions are transformed into equivalent expressions.

2. Anisotropic Hardening Behaviors

Before introducing the modeling for the anisotropic hardening behaviors, an overall understanding of the deformation-induced anisotropic hardening behaviors is essential. At the beginning, a detailed decomposition of the hardening behaviors is introduced. After that, decomposed hardening behaviors possibly to be observed in specific materials and forming processes are investigated. Finally, a brief review work of typical experiments with loading path changes to detect the materials’ anisotropic hardening behaviors is conducted.

2.1. Decomposition of the Anisotropic Hardening Behaviors

To get a comprehensive description of the anisotropic hardening behaviors, a phenomenological decomposition needs to be made based on the modes loading of loading path changes. Different modes of loading path changes may activate different slip systems [20], and macroscopically result in different hardening behaviors.

The description of the loading mode relies on the loading directions before and after the loading path changes. The current plastic strain increment $d\mathit{\epsilon }$ or stress deviator $\mathit{s}$ is used to define the loading direction. The terms ${\mathit{N}}_{\mathit{\epsilon }}$ and ${\mathit{N}}_{\mathit{s}}$ (or ${\mathit{N}}_{\mathit{s}-\mathit{\alpha }}$), which represent the normalized directions of $d\mathit{\epsilon }$ and $\mathit{s}$ (or $\left(\mathit{s}-\mathit{\alpha }\right)$), respectively, are expressed as follows:

$${\mathit{N}}_{\mathit{\epsilon }}=\frac{d\mathit{\epsilon }}{\Vert d\mathit{\epsilon }\Vert }$$

(1)

$${\mathit{N}}_{\mathit{s}}=\frac{\mathit{s}}{\Vert \mathit{s}\Vert } (\mathrm{or} {\mathit{N}}_{\mathit{s}-\mathit{\alpha }}=\frac{\mathit{s}-\mathit{\alpha }}{\Vert \mathit{s}-\mathit{\alpha }\Vert })$$

(2)

where ‖ ‖ denotes the norm operation of the inner tensor; $\mathit{\alpha }$ is the so-called back stress denoting the center of the yield surface and this term is explained in Section 3. If the directions of ${\mathit{N}}_{\mathit{\epsilon }}$ and ${\mathit{N}}_{\mathit{s}}$ need not be distinguished, the term $\mathit{N}$ is used as a general representation of the loading direction. For simplicity, we first consider a case of two-stage loading with an abrupt loading path change. Schmitt et al. [21] proposed a method for measuring the abrupt loading path changes, which is widely used. The Schmitt factor $\cos \mathsf{\chi }$ is expressed as a scalar product of the normalized direction tensors before and after loading path changes:

$$\cos \mathsf{\chi }={\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}:\mathit{N}$$

(3)

where ${\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}$ and $\mathit{N}$ respectively denotes the loading direction before and after the loading path changes and $\mathsf{\chi }$ refers to the Schmitt angle. In the origin version of Schmitt et al. [21], $\mathit{N}$ is given by Equation (1), while in Barlat’s work [22], $\mathit{N}$ is given by Equation (2) to simplify their model formulation. Different loading modes can be easily specified by the Schmitt factor: $\cos \mathsf{\chi }=1$ represents no path changes (i.e., proportional loading); $\cos \mathsf{\chi }=-1$ or $0$ represents reverse (e.g., compression after tension loading) or orthogonal path changes (e.g., torsion loading after axial tension loading of the round tube), respectively; and $\cos \mathsf{\chi }\ne 1,-1, 0$ represents the remaining intermediate modes of loading path changes. For the more universal case of continuous strain path changes, the incremental descriptions of Equation (3), which is often applied for the finite element (FE) method, will be used. In this situation, $\mathit{N}$ denotes the current loading direction, while ${\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}$ is expressed by the deformation history related state variables. More details of such descriptions are presented in the beginning of Section 3.2.2.

To better introduce the hardening behaviors, the concept of the hardening modulus, $h,$ is adopted, which is expressed as $h=d\sigma /d\overline{\epsilon }$. In which, $\overline{\epsilon }$ is the equivalent plastic strain calculated according to the principle of plastic work equivalence the hardening modulus. According to the sign of the hardening modulus, the hardening behaviors are described by hardening ($h>0$), stagnation ($h\approx 0$), and softening ($h<0$). The hardening behaviors are decomposed according to three typical modes of loading path changes. As is mentioned before, isotropic hardening is assumed in proportional loading ($\cos \mathsf{\chi }=1$). The other two typical loading modes are reverse loading ($\cos \mathsf{\chi }=-1$) and orthogonal loading ($\cos \mathsf{\chi }=0$). The typical anisotropic hardening behaviors corresponding to these two loading modes are respectively displayed in Figure 1a,b. The anisotropic hardening behaviors under reverse loading are listed as follows:

(i)

The Bauschinger effect: The lower re-loading yield stress than the unloading yield stress, and the higher hardening rate than that of the monotonic loading.

(ii)

Hardening stagnation: A temporary reduction in the hardening rate, which introduces a plateau into the stress–strain curve near the end of the Bauschinger effect.

(iii)

Permanent softening: The lower stress level than monotonic loading under a large strain, even when the hardening rate recovers to the level of monotonic loading. Permanent softening is usually supposed to be caused by the two aforementioned transient hardening behaviors.

The anisotropic hardening behaviors under orthogonal loading are listed as follows:

(i)

Orthogonal hardening or orthogonal softening: When the load path changes orthogonally, the stress level increases or decreases compared with the one under monotonic loading.

(ii)

Permanent softening: Permanent softening can also occur in orthogonal loading, with the phenomenon similar to that observed in reverse loading.

The differences in the hardening behaviors triggered by reverse loading and orthogonal loading are usually physically explained through the dislocation evolution. In reverse loading, dislocation motion persists in identical slip systems despite the reversal of the slip direction. In orthogonal loading, the original slip system is blocked, and dislocation motion transfers to new activated slip systems.

In the remaining modes of loading path changes, the anisotropic hardening behaviors is observed as the transition state between the proportional, orthogonal, and reverse loading. The reproduction of the hardening behaviors in the intermediate state of loading path changes is fulfilled by the description of the yield surface during the deformation process. The approximate geometric shape of the yield surface after deformation is depicted in Figure 2b. This surface forms a corner with high curvature along the direction of pre-strain and a flat zone at its rear. The distorted yield surface is often represented as an “egg-shaped” surface.

In the preceding part of the text, the anisotropic hardening behaviors are sorted according to the features of stress evolution; however, the modeling of the deformation direction and extent is also crucial. In the engineering field, the Lankford coefficient (also called r-value) is a macroscopic reflection of the plastic deformation direction. And the evolution of r-value is counted as a kind of anisotropic hardening behaviors which is supposed to be modeled. Although the r-value can be accordingly predicted by the yield function of every moment under the hypotheses of associated flow, whether the prediction result can faithfully reproduce the evolution of r-value remains ambiguous. Since the modeling for the r-value evolution is another topic which is not involved much in this paper, interested readers can refer to the work of Mánik et al. [24], Lee et al. [25], and Zaman et al. [26].

2.2. Influences from Materials and Forming Processes

Section 2.2 summarizes the typical anisotropic hardening behaviors, but not all these behaviors are pronounced in the specific materials and forming processes. Moreover, a hardening model is usually proposed according to the specific materials and forming processes. This section provides an overview of the hardening behaviors possibly observed in the different materials and forming processes.

Currently, most anisotropic hardening models are specific for steels and aluminum alloys. Typical hardening behaviors existing in these materials are summarized here. Silvestre [21] compared the cyclic hardening behavior of different steel families and found that the Bauschinger effect is more pronounced in multiphase steels with higher strength. Yoshida [27] conducted the forward-reverse loading test of different steels and concluded that the permanent softening behavior after loading reverses is more pronounced in the mild steel than that observed in the advanced high strength steel. Moreover, the hardening stagnation after loading reverses is hardly observed in high strength steels and aluminum with high alloy components [28]. Khan et al. [29,30] systematically studied the deformation-induced hardening behaviors of aluminum alloys and found that the orthogonal softening tends to take place in low work hardening aluminum alloy while the orthogonal hardening is often observed in high work hardening aluminum alloy. In general, the commercial pure aluminum exhibits the similar anisotropic hardening effect after loading path changes with mild steel [31], but the degree of anisotropic hardening of commercial pure aluminum is weaker [24]. Among the materials investigated, AA6XXX (six series aluminum alloy) shows mild anisotropic hardening behaviors [32]. The anisotropic hardening behaviors observed in different materials are summarized in

Table 1. Anisotropic hardening behaviors observed in typical steels and aluminum alloys.

Loading Modes	Hardening Behaviors	Materials
		Steel	Aluminum Alloy (AA)
Reverse loading cos χ = −1	The Bauschinger effect	pronounced in multiphase steels with high strength	moderate in AA6XXX
	Hardening stagnation	hardly observed in high strength steels Mild, IF, DP (moderate)	hardly observed in AA with high alloy components AA1XXX, AA2XXX, AA3XXX
	Permanent softening	pronounced in mild steel	hardly observed in AA6XXX
Orthogonal loading cos χ = 0	Orthogonal hardening	observed in mild steel Mild, IF	occur in low work hardening Al alloy AA1XXX, AA2XXX, AA3XXX
	Orthogonal softening	observed in high strength steels DP, TRIP	observed in high work hardening Al alloy AA5XXX, AA6XXX
	Permanent softening	pronounced in mild steel Mild, IF	hardly observed in AA6XXX

Note: Materials involved in the table: Steel: mild steel [16,17,29,31,32,36], IF (interstitial free) steel [32], DP (dual-phase) steel [23,32,34], TRIP (transformation induced plasticity) steel [23,33]; Aluminum alloy: AA1XXX [16,17,29,30], AA2XXX [37], AA3XXX [20], AA5XXX [32], AA6XXX [16,17,28,31].

. Moreover, some representative materials with each decomposed hardening behavior are listed according to the published research (Vincze et al. [31], Haddadi et al. [32], Rauch et al. [16,17], Khan et al. [29,30], Ha et al. [33], Mánik et al. [24], Zhou et al. [34], Qin et al. [20], Lee et al. [35]).

As a footnote, the loading modes defined by the Schmitt factor may be not enough. The materials’ anisotropic hardening behaviors are also influenced by some other potential process factors. In other words, some material’s hardening behaviors may not strictly correspond to those listed in Table 1. For example, the experiments of Qin et al. [20] found that the orthogonal hardening behaviors of AA3103 after pre-tension and pre-rolling are not fully consistent, despite the similar Schmitt angles. The influence of other factors related to forming processes, for example, the hydrostatic pressure, may not be ignored. Anyhow, at present, the Schmitt angle still acts as the dominant factor to define the loading modes in the modeling of anisotropic hardening behaviors. By referring to Table 1, the modeling objectives can be preliminarily specified according to the materials and the characteristics of forming processes, while the final determination of the hardening model relies on the further experimental result.

2.3. Experiments with Loading Path Changes

A series of experiments with loading path changes are designed to study the materials’ hardening behaviors. According to the loading condition, the experiments are classified into proportional loading, reverse loading (or cyclic loading), orthogonal loading or the experiments in the remaining loading modes; the experiments with abrupt or continuous loading path changes; the experiments with uniaxial loading or multiaxial loading. According to the specimens, the experiments are classified into sheet specimen and thin-walled circular tube specimen. Moreover, some practical forming experiments are applied for the model calibration and validation.

Representative experiments used by different researchers to study the materials’ anisotropic hardening behaviors are sorted out in this paragraph. The uniaxial experiments are widely used because of their simplicity. The basic method is uniaxial tension-compression cyclic test [24,33]. To measure the stress–strain in an adequate deformation range, an anti-buckling device must be specially designed. Some anti-buckling devices with different geometries are summarized in the review work of Bruschi et al. [38]. Alternatively, the reverse shear tests are adopted to provide the homogenous deformation in a larger strain range, which include the forward and reverse simple shear test [16] and in-plane torsion test [39]. To investigate the materials’ hardening behaviors under orthogonal or other abrupt loading path changes, multi-step loading tests are designed. In the first loading step, a pre-strain is carried out in the specimen with large geometry size, and then the specimens with small size are fabricated from the homogenous deformation region of the pre-strained large specimens to conduct the subsequent deformation in other loading paths. Common pre-strain methods include the rolling [24], uniaxial tension [40], plane strain tension [41], simple shear and bulging [32]. Furthermore, in some newly published works, three-step loading tests with two loading path changes are carried out to study the hardening behaviors in more complex loading conditions [22,38]. For the convenience of reproducing the continuous loading path changes between more varied loading modes, the multiaxial testing methods with special setups are developed. Shear-tension test on plane stress specimens [42,43,44] and biaxial tension test on cruciform specimens [45,46] are two typical multiaxial testing methods for sheet metals, and these two methods closely reflect the practical loading condition in forming processes. Furthermore, as for the tube metals, tension–torsion test on thin-walled cylindrical tubes [29,30,47,48] is specially designed.

To some extent, certain experimental factors also have influences on the materials’ hardening behaviors, such as the definition of the yield point, the method of probing the yield surface, the specimen types, and their corresponding experiment methods [49]. Common probing methods of yield surface include the equivalent plastic work method [50], post unloading proportional limit or strain offset method [29,30], and abrupt strain path changes method [51,52,53]. There may be deviation in yield surface shape measured by different test methods. For example, the yield always encloses the origin of the σ-τ space in the total unloading method, whereas the origin may be outside the yield surface in the partial unloading method [49] (Figure 3). Different potential yield surfaces can be observed depending on the material properties of the specimen and testing method used. A more detailed review of the experimental methods to probe the distorted yield together with the modeling methods for this feature is recently carried out by Marek et al. [54].

Based on the result of anisotropic hardening experiments, hardening models are calibrated. Optimal fitting is the universal calibration method for anisotropic hardening models. Since model parameters are not completely decoupled in most cases, the optimal procedure is often conducted by fitting several stress–strain curves measured in different loading paths simultaneously. As for some distortion hardening models where the shape prediction of the subsequent yield surface is emphasized, a set of subsequent yield points of certain pre-strains is also used for the parameter identification [56,57]. In some models, the constrain relation of partial parameters can be determined or initialized in advance according to some features of the experimental data [24,54,58]. Moreover, the measurable shape and size features of simple forming parts can be applied for the inverse determination of model parameters [9].

Some other experiment results are applied to validate the models’ effectiveness to reproduce the anisotropic hardening behaviors. In most cases, the redundant data acquired from the aforementioned calibration experiments are used for model validation. Commonly, the r-value versus strain curves [25,26] and the stress–strain curves acquired from the similar tests but under different loading paths or different level of pre-strains [22,23] are two resources of data applied for the model validation. Whereas, the use of forming experiments which reproduce the practical working conditions for model validation seems to be more convincing. The FE simulation results are compared with the corresponding data measured in the forming processes. Yoon et al. [59] verify the accuracy of their hardening model by comparing the simulated and practical measured loading-stroke curve in the loading process and the springback parameters of U-draw bending parts after unloading. Moreover, a complex forming simulation of a B-pillar is also conducted for comparison, and based on the comparison of geometrical features, the reliability of the hardening model is further verified.

In general, the deformation-induced anisotropic hardening behavior is rather complex in the forming processes of different materials. Thus, the features of materials and processes must be thoroughly considered before designing the experiments and selecting or constructing the hardening models. The universal regulations summarized in this section are for reference only, while the reliable modeling of hardening behaviors remains an in-depth study of the target materials and processes.

3. Mathematical Fundamentals of Anisotropic Hardening Models

Although the specific expressions for different anisotropic hardening models are quite varied, the core mathematic ideas behind these models are similar. In this section, typical modeling ideas and strategies are refined and summarized, to help readers better understand the modeling process.

Anisotropic hardening models are developed from isotropic hardening models. The following expression is a universal framework for the temperature and rate independent isotropic hardening model:

$$\psi \left(\mathit{s}\right)={\mathsf{\sigma }}_f\left(\overline{\epsilon }\right)$$

(4)

where $\overline{\epsilon }$ is the equivalent plastic strain calculated according to the principle of plastic work equivalence; $\psi \left(\mathit{s}\right)$ denotes the yield function, which is unified to a homogeneous function of degree one in this formula; and ${\sigma }_f\left(\overline{\epsilon }\right)$ is a reference monotonic flow stress curve usually measured in the rolling direction, which reflects the basic hardening trend of monotonic loading and indicates the homogenous expansion of the yield surface.

In general, the extension of Equation (4) from isotropic hardening to anisotropic hardening involves the following two steps:

(i)

In the first step, state variables are introduced into the original hardening function (Equation (4)) to describe the shape of the distorted yield surface. In this paper, $\mathit{X}$ represents the set of internal state variables.

(ii)

In the second step, the evolution functions are constructed for the newly introduced state variables, to memorize the loading history and reproduce the corresponding hardening behaviors throughout the deformation process. The evolution function is usually expressed as $d\mathit{X}\left(\mathit{X},\mathit{s},d\mathit{\epsilon }\right)$ which reflects the influences of the loading history ($\mathit{X}$) and the current loading condition ($\mathit{s},d\mathit{\epsilon }$).

The remaining part of Section 3 presents the mathematical description of the aforementioned two steps.

3.1. Strategies of Introducing the State Variables

The state variables are introduced into a hardening model to incorporate the influences of the loading history on the hardening behavior. Depending on the object to be modified, the methods to introduce state variable are classified into two categories. In the first class of methods, the variables are incorporated into the left-hand side of Equation (4) to adjust the expression of yield function ($\psi \left(\mathit{s}\right))$, so that the evolution of the yield surface during the plastic deformation is reproduced. The second class of methods is relatively simple in which the reference stress–strain curve (${\mathsf{\sigma }}_f\left(\overline{\epsilon }\right)$), namely the right-hand side of Equation (4), is modified after the loading direction changes while the description of the yield surface remains unchanged.

In the modeling of plasticity, the associative flow rule is commonly used to determine the direction of plastic strain increment ($d\mathit{\epsilon }$) as follows:

$$d\mathit{\epsilon }=d\lambda \frac{\partial \psi }{\partial \mathit{\sigma }}$$

(5)

where $d\lambda$ represents a non-negative scalar determined according to the consistency condition. Equation (5) reveals that the direction of plastic strain increment ($d\mathit{\epsilon }$) is consistent with that of the outer normal of the current yield surface. As the yield surface distorts with plastic deformation in the aforementioned first method, the direction of plastic strain increment is accordingly reproduced so that the evolution of the r-value can be recorded. Moreover, the associated flow indicates that ${\mathit{N}}_{\mathit{\epsilon }}$ and ${\mathit{N}}_{\mathit{s}}$ (expressed in Equations (1) and (2), respectively) are completely consistent for the von Mises yield function, but the deviations will occur between these two definitions when other yield functions are used. The following section describes the specific implementation of the two methods to introduce the state variables.

3.1.1. Modification of the Yield Function

Depending on the model capability to describe different hardening behaviors, the modeling methods can be further classified into two categories: kinematic and distortional hardening models. The schematics of the yield surface evolutions subject to these two hardening models are displayed in Figure 2. The yield surface for isotropic hardening is also depicted as a reference. Kinematic hardening models are specially constructed to reproduce the hardening behaviors occurring under forward and reverse loading, which correspond to the stress states represented by two red points in Figure 2a. Whereas, the hardening behaviors occurring in other loading modes between forward and reverse loading (Figure 2a) are not specially captured by this group of models. Distortional hardening models contain the description of the hardening behaviors under both reverse and orthogonal loading (two blue points in Figure 2b). Some distortional hardening models provide a more comprehensive description of all the possible anisotropic hardening behaviors mentioned in Section 2 by tailoring the overall shape of the yield surface.

Kinematic hardening is a classic method for modeling the anisotropic hardening behaviors, and this method has been used in industry. The pure kinematic hardening model assumes that the yield surface translates but not expands homogeneously during the plastic deformation. In practice, a combined description of the translation and expansion is adopted in most cases (the red circle in Figure 2a); note that this description is still classified as kinematic hardening in this paper for convenience. The tensor variable $\mathit{\alpha }$, which denotes the back stress, is introduced into the yield function to trace the motion trail of the yield surface center. The uniform expression of kinematic hardening models is as follows:

$$\psi \left(\mathit{s}-\mathit{\alpha }\right)={\sigma }_{fi\mathit{s}o}\left(\overline{\epsilon }\right)$$

(6)

Note that the notation in Equation (6) is different from that in Equation (4), since when the back stress $\mathit{\alpha }$ is inserted into the hardening model, the reference stress–strain relationship is different from the original one used in isotropic hardening models because the back stress affects the yield even under monotonic loading. In this paper, this newly defined reference curve is denoted as ${\sigma }_{fi\mathit{s}o}\left(\overline{\epsilon }\right)$. According to the kinematic hardening behavior of different materials, many evolution laws have been proposed for $\mathit{\alpha }$ [60,61], and an introduction of these laws is presented in Section 4.1. Since the kinematic hardening models are specialized for forward and reverse loading, they cannot accurately reproduce the anisotropic hardening behaviors occurring for the remaining loading modes. For example, even though the description of yield surface translation indirectly introduces the effect of orthogonal softening (the yield surface shrinks in the direction orthogonal to the current stress; marked by the green arrow in Figure 2a) into the model, the accuracy of the description cannot be guaranteed without any specialized variables. Thus, the kinematic hardening models are mainly used for modeling hardening in processes dominated by reverse loading, such as cyclic loading and springback. In the anisotropic hardening models, the back stresses often act as the basic variables to reproduce the material behaviors under reverse loading.

Distortional hardening models can flexibly adjust the shape of the yield surface in the deformation process. Thus, they can reproduce the hardening behaviors occurring along different loading paths simultaneously. Early proposed hardening modeling methods tend to incorporate the deformation-induced anisotropy by modifying the yield function. Baltov and Sawczuk [62] introduce a fourth-order anisotropic tensor $\mathit{A}$ into the quadratic Mises isotropic yield criterion to reproduce the distortion features:

$$\psi \left(\mathit{s}-\mathit{\alpha }, \mathit{A}\right)={\left[\left(\mathit{s}-\mathit{\alpha }\right):\mathit{A}\left(\mathit{s}-\mathit{\alpha }\right)\right]}^{\frac{1}{2}}={\mathsf{\sigma }}_{fi\mathit{s}o}\left(\overline{\epsilon }\right)$$

(7)

In this formulation, the anisotropic fourth-order tensor $\mathit{A}$ acts as a state variable, so that the anisotropy yield surface is distorted in the deformation process with the evolution of this variable. This modification method originates from a universal modeling strategy for initial orthotropy yield of pre-rolled blanks [63,64]. The difference is that $\mathit{A}$ keeps constant in yield models, while evolves with deformation in hardening model. The back stress $\mathit{\alpha }$ from kinematic hardening models is also incorporated in Equation (7). First, $\mathit{\alpha }$ reproduces the hardening behaviors under reverse loading. Second, this variable memorizes the loading history and controls the evolution of $\mathit{A}$. Under most circumstances, $\mathit{A}$ is orthogonally decomposed to describe the distortion of the yield surface in the directions parallel (${\mathit{A}}_p$) and orthogonal (${\mathit{A}}_o$) to the current loading direction $\mathit{N}$:

$${\mathit{A}}_p=\left(\mathit{N}:\mathit{A}\mathit{N}\right)\mathit{N}\otimes \mathit{N}$$

(8)

$${\mathit{A}}_o=\mathit{A}-{\mathit{A}}_p$$

(9)

The symbol “$\otimes$” denotes the tensor product operation. As explained by Norman [65], the parallel and orthogonal components are supposed to reflect the influences of the current dynamic and latent slip systems on hardening behaviors, respectively. Moreover, the norms of the fourth-order tensors can be used to adjust the reference flow stress curve to reproduce some transient hardening behaviors, which are introduced in Section 4.

Although tensor $\mathit{A}$ can comprehensively describe the anisotropy of the yield surface, the format of this fourth-order tensor is a bit complicated. Moreover, this tensor has weak compatibility with different yield functions. Some alternative tensors with more intelligible forms have been proposed. The hardening models based on these tensors conduct intuitionistic geometric adjustments of the yield surface. Because of their simplicity and compatibility, these models are useful in practical applications.

In the first approach, a distorted stress tensor ${\mathit{s}}_d$ is used to replace the original stress tensor $\mathit{s}$. The modified model is expressed as follows:

$$\psi \left({\mathit{s}}_d-\mathit{\alpha }\right)={\sigma }_{fi\mathit{s}o}\left(\overline{\epsilon }\right)$$

(10)

The distorted stress tensor ${\mathit{s}}_d$ is typically constructed through the procedure described in the following text. First, the current stress $\mathit{s}$ is subjected to orthogonal decomposition according to the loading direction before loading path changes (${\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}$) as follows:

$$\mathit{s}={\mathit{s}}_p+{\mathit{s}}_o$$

(11)

$${\mathit{s}}_p={\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}\Vert \mathit{s}\Vert \cos \mathsf{\chi }$$

(12)

$${\mathit{s}}_o={\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}\Vert \mathit{s}\Vert \sin \mathsf{\chi }$$

(13)

where ${\mathit{s}}_p$ represents the stress component parallel to ${\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}$ and ${\mathit{s}}_o$ denotes the orthogonal component. The parameter ${\mathit{s}}_p$ or ${\mathit{s}}_o$ is then adjusted in different manners to incorporate the description of corresponding hardening behaviors. Finally, a new expression of stress reconstructed by these two parameters is defined as the distorted stress tensor ${\mathit{s}}_d$. The distortion of the yield surface in specific directions can be reproduced by the aforementioned parameters. For example, Barlat et al. [66] adjusted ${\mathit{s}}_o$ to reproduce the orthogonal hardening behavior (details will be introduced in Section 4.4), and François [36] adjusted $\mathit{s}$ to reproduce the “egg-shaped” distortion of the yield surface (details will be introduced in Section 4.6).

The second approach introduces the incorporating plane equations into the yield function to modify the shape of the yield surface in specific directions. Geometrically, the yield locus in the π-plane is locally compressed by an invisible asymptotic line (the blue line in Figure 4). In this paper, the original yield function is denoted as ${\psi }_0$, and the modified yield function is denoted as $\psi$, which is expressed as follows:

$$\psi ={\left({\psi }_0^q+{\psi }_h^q\right)}^{\frac{1}{q}}={\mathsf{\sigma }}_f\left(\overline{\epsilon }\right)$$

(14)

where $q$ is a coefficient and ${\psi }_h$ denotes the function of the asymptotic line. When $q=2$, the expression of $\psi$, which controls the distorted region of the yield surface, is equivalent to an elliptic function. This expression often acts as an alternative for the back stress in pure distortional hardening models; therefore, the reference stress–strain curve can be represented by ${\sigma }_f\left(\overline{\epsilon }\right)$ in Equation (14) without adjustments. The aforementioned elliptic function conducts a localized adjustment of the yield surface shape but not a monotonic translation; thus, it can independently capture the flattening of the yield surface in the direction opposite or orthogonal to ${\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}$ and avoid the coupling of the Bauschinger effect and orthogonal softening.

3.1.2. Modification of the Reference Stress–Strain Curve

Another method for modeling the anisotropic hardening involves modifying the reference flow stress curve during the deformation process. The general expression of this method is as follows:

$$\mathsf{\psi }\left(\mathit{s}\right)={\mathsf{\sigma }}_f\left(\overline{\epsilon },\mathit{X}\right)$$

(15)

where $\mathit{X}$ denotes a set of state variables influenced by the loading history. In common anisotropic hardening models, $\mathit{X}$ is multiplied or added to the original flow stress in the form of a norm, that is, $\mathsf{\psi }\left(\mathit{s}\right)=\Vert \mathit{X}\Vert {\mathsf{\sigma }}_f\left(\overline{\epsilon }\right)$ or $\mathsf{\psi }\left(\mathit{s}\right)={\mathsf{\sigma }}_f\left(\overline{\epsilon }\right)\pm \Vert \mathit{X}\Vert$. In such expressions, the movement of state variables to the left-hand side is nominally equivalent to modifying the yield function. However, the actual shape of the yield surface is not changed through this modification; thus, the direction of plastic strain increment cannot be well captured. In practical applications, the aforementioned modification method often acts as an alternative for or a supplement to the methods mentioned in Section 3.1.1. In the models proposed by Teodosiu and Hu [67] and Wang et al. [68], the norm of a fourth-order tensor $\Vert \mathit{S}\Vert$ plays the role of $\Vert \mathit{X}\Vert$ to incorporate the influences of transient hardening on the flow stress curve or control the saturation of the back stress. The term $\mathit{S} \left(={\mathit{S}}_p+{\mathit{S}}_o\right)$ is orthogonally decomposed according to the decomposition of $\mathit{A}$ expressed in Equations (8) and (9). The models of Mánik et al. [24] and Qin et al. [69] use some straightforward scalar variables as substitutes for the aforementioned tensor and their methods exhibit favorable portability.

3.2. Methods of Constructing Evolution Functions

After internal state variables ($\mathit{X}$) are introduced into the yield function, their evolution rules must be specified to memorize the deformation history and reproduce the relevant transient hardening behaviors. This section describes a typical evolution mode and a general classification of these variables.

Both the directionality and the severity of the yield surface translation or distortion ought to be captured by the state variable set $\mathit{X}$. According to their capability to describe the direction, the state variables might as well be classified into the directivity state variables ${\mathit{X}}_d$ and the magnitude state variables $X_i$. Various combinations of these two classes of state variables are adopted by different models. For example, the Armstrong and Frederick [70] kinematic hardening model uses only one deformation history-related variable, namely the back stress $\mathit{\alpha }$, which can simultaneously represent the direction and magnitude of the Bauschinger effect. Barlat’s distortional hardening model [22] uses a microstructure deviator $\widehat{\mathit{h}}$ and a series of scalar variables $g_i$ to respectively record the direction and magnitude of different anisotropic hardening behaviors. Among the anisotropic models, the second-order tensors variables often take the role of the directivity state variables; while the scalar variables and the fourth-order tensor variables often act as the magnitude state variables. As the distortion of the yield surface is oriented (see in Figure 2b), after the loading direction switches, the distortion of the yield surface may accumulate in another mode. Thus, the evolution of the magnitude state variables is always based on the directivity state variables.

3.2.1. A Universal Evolution Function for the State Variables

Generally, in two-stage loading with an abrupt strain path change, the two-phase features must be modeled. In the first phase, the potential of anisotropic hardening is accumulated, but the phenomenon of anisotropic hardening will not exhibit until strain path changes. Usually, higher pre-strain results in greater anisotropic hardening at the onset of loading path changes, while the increase rate of anisotropic hardening tends to slow as pre-strain increases. In the second phase, after the strain path switches, the degree of anisotropic hardening increases considerably and then decreases gradually as deformation continues, with the hardening rate gradually recovering to the level of proportional loading [71]. The general evolution laws of the accumulation and decay of anisotropic hardening are rather similar. The nonlinear saturated incremental evolution function adopted by Voce [72] and Armstrong and Frederick [70] is used in most anisotropic hardening models as the evolution function for different variables. The universal incremental equation has the following scalar and tensor forms, respectively:

$$dX=k\left(X_{sat}-X\right)d\overline{\epsilon }$$

(16)

$$d\mathit{X}=k\left(X_{sat}{\mathit{N}}_{\epsilon }-\mathit{X}\right)d\overline{\epsilon }$$

(17)

In Equation (16), $X_{sat}$ denotes the saturation value of $X$, and $k$ controls the evolution rate of $X$. In Equation (17), $X_{sat}$ and ${\mathit{N}}_{\epsilon }$ denote the saturation magnitude and direction of $\mathit{X}$, respectively. By integrating the scalar function under proportional deformation, the following explicit evolution function is obtained:

$$X=X_0+\left(X_{sat}-X_0\right)\left[1-\exp \left(-k\overline{\epsilon }\right)\right]$$

(18)

where $X_0$ denotes the initial value of $X$. The influences of $X_{sat}$ and $k$ on the evolution function are illustrated in Figure 5. Most anisotropic hardening behaviors can be roughly modeled by using the simple evolution function (Equation (17)) with two parameters. However, for the more precise reproduction of some hardening behaviors, especially the Bauschinger effect (details will later be illustrated in Section 4.1), a combination of several functions in the form of Equation (17) is used [58]. Other special forms of evolution functions are introduced in Section 4, when they occur in the corresponding hardening models.

3.2.2. Evolutions of the Directivity Variables

This section presents a detailed discussion on the directivity variables. First, a more general case of loading path changes, namely continuous loading path changes, is described. For continuous loading path changes, the original definition of ${\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}$ expressed in Equation (3) is no more sufficient, and the Schmitt angle must be redefined. Physically, the shape of the current yield surface is mainly affected by the materials’ microstructures, which gradually develops during the deformation process. As the second-order directivity variable ${\mathit{X}}_{\mathit{d}}$ acts as the memory of loading history, its direction can be used to redefine ${\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}$ as follows:

$${\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}=\frac{{\mathit{X}}_{\mathit{d}}}{\Vert {\mathit{X}}_{\mathit{d}}\Vert }$$

(19)

Because the microstructure evolution also exhibits transient and saturated features, Equation (17) still acts as the basic framework for the evolution of ${\mathit{X}}_{\mathit{d}}$.

Depending on whether their magnitudes are varied, the directivity variables can be further classified into two categories: variables whose magnitudes evolve during the plastic deformation and variables that only represent the directional information. The first category of variables includes the nonlinear back stress tensor $\mathit{\alpha }$ proposed by Armstrong and Frederick [70] and the dimensionless delayed pointer $\mathit{P}$ used in the models of Teodosiu and Hu [67], Mánik et al. [24], and Qin et al. [23]. The evolution function of such variables is completely consistent with Equation (17). Figure 6a illustrates the evolution mode of the aforementioned variables in the π-plane under two-stage loading. In the first loading stage, $\mathit{P}$ gradually increases from 0 along the loading direction. After the loading direction switches, $\mathit{P}$ evolves along the connecting line of $\mathit{P}$ at the end of the first loading stage and corresponding saturation limit of the second loading stage. Because the magnitudes of variables change simultaneously during this process, the influence of microstructures on anisotropic hardening can be simultaneously represented. Thus, in the basic kinematic hardening models, orthogonal hardening behaviors can be reproduced using only one variable, namely $\mathit{\alpha }$.

Whereas the effect of microstructure evolution on different hardening behaviors is different. Thus, even though the magnitudes of the aforementioned directivity variables are variable, to capture different anisotropic hardening behaviors, corresponding magnitude variables are still necessary. In the other group of directivity variables, the amplitude information is completely omitted (which is independently represented by different magnitude variables) and only the directional information is preserved. The use of such variables avoids the unnecessary coupling, and makes the descriptions of different hardening behaviors relatively independent. Barlat et al. [22] proposed an evolution function for this category of tensor variables with unit magnitude. The following function describes the evolution rule of this tensor:

$$d\widehat{\mathit{h}}=k\left(\mathit{N}-\widehat{\mathit{h}}\cos \mathsf{\chi }\right)d\overline{\epsilon }$$

(20)

where $\widehat{\mathit{h}}$ denotes the directivity variable with a unit magnitude and ${\mathit{N}}_{\mathit{p}\mathit{r}\mathit{e}}=\widehat{\mathit{h}}$ in this framework. An evolution defined by Equation (20) under two-stage loading is illustrated in Figure 6b. In the first incremental step of plastic deformation, $\widehat{\mathit{h}}$ is immediately directed along the deformation direction. This variable keeps constant as long as no loading path changes occur. After the loading path switches, the direction of $\widehat{\mathit{h}}$ changes toward the current loading direction $\mathit{N}$ along an arc in the plane determined by $\widehat{\mathit{h}}$ and $\mathit{N}$. However, as is pointed out by Manopulo et al. [73] and Holmedal [74], some drawbacks exist in this evolution mode. Any slight fluctuation of loading at the beginning of a simulation can result in an unreasonable initialization of $\widehat{\mathit{h}}$, which will influence the subsequent simulation of the hardening behaviors to some extent.

4. Modeling of the Decomposed Anisotropic Hardening Behaviors

This section focuses on how the anisotropic hardening behaviors are reproduced in different models. At first, different hardening models’ description capabilities of each anisotropic hardening behavior are compared. Then, the specific modeling of each decomposed hardening behavior is reviewed in order. Models mainly involved in this section are respectively proposed by Armstrong and Frederick [70] (AF 66), Chaboche [75] (Chaboche 86), Yoshida and Uemori [61] (YU 02); Baltov and Sawczuk [62] (Baltov 65), Rokhgireh et al. [76,77] (Rokhgireh 13, 17), Voyiadjis and Foroozesh [78] (Voyiadjis 90), Feigenbaum and Dafalias [56,79] (Feigenbaum 07, 08), Noman et al. [65] (Noman 10), Pietryga et al. [57] (Pietryga 12), François [36] (François 01), Zhang et al. [80] (Zhang 20); Teodosiu and Hu [67] (TH 98), Wang et al. [68] (Wang 08), Mánik-Holmedal-Hopperstad [24] (MHH 15), Qin-Holmedal-Hopperstad [23] (QHH 18); homogeneous anisotropic hardening [22,66,81,82] (HAH 11, 13, 14, 20), Qin et al. [69](Qin 17), He et al. [83] (He 18) and Holmedal [74] (Holmedal 19) models.

Table 2. Capability of different models to describe each anisotropic hardening behavior.

Hardening Model		Type	Anisotropic Hardening Behaviors							Number of Model Parameters
			Reverse Loading			Orthogonal Loading			Others
			The Bauschinger Effect (Reproduced by Multi Variables)	Hardening Stagnation	Permanent Softening	Orthogonal Softening	Orthogonal Hardening	Permanent Softening	“Egg-Shaped” Distortion
Isotropic (i)	Voce 48 [72]	i								2
Kinematic (k)	AF 66 [70]	k				○				4
	Chaboche 86 [75]	k	√		○	○				7
	YU 02 [61]	k	√	√	○	○				6
Distortional (d)	Baltov 64 [62]	k+d	Rokhgireh 17 [77]						Rokhgireh 13 [76]	5/6
	Voyiadjis 90 [78]	k+d				√	√		√	6/7
	Feigenbaum 07 [56], 08 [79]	k+d				○			√	6
	Noman 10 [65]	k+d				○	√		Pietryga 12 [57]	8/9
	François 01 [36] & Zhang 20 [80]	k+d				○			√	5/6
Varied flow stress (vfs)	TH 98 [67] & Wang 08 [68]	k+vfs		√		√	√			12
	MHH 15 [24]	k+vfs	○	√		○	√	√		11
	QHH 18 [23]	k+vfs +d	○	√	√	○	√	√		13
Pure distortional	HAH 11 [22], 13 [81], 14 [66]	d+vfs		Qin 17 [69]	√	√	√	Qin 17 [69]	He 18 [83]	9/11/13/16/19
	Holmedal 19 [74]	d+vfs			√	○			○	12
	HAH 20 [82]	d+vfs			√	√	√		√	25

Note: “√” indicates that a specialized variable is included in the hardening model to reproduce the corresponding anisotropic hardening behavior; “○” indicates that the hardening behavior is indirectly described by other variables, which means the description of a behavior is coupled with the descriptions of another one; labels such as “Rokhgireh 17” in the table indicate that the description of corresponding hardening behavior is supplemented in “Rokhgireh 17” model; for comparison purpose, each model in this table is supposed to contain the isotropic hardening component, which is here unified as the form of Voce 48 (with 2 parameters, where initial yield strength is not viewed as a hardening parameter); the number of back stresses in MHH 15 and QHH 18 here are unified as 1.

presents the capabilities of the aforementioned models to describe different anisotropic hardening behaviors. Since the Bauschinger effect is a hardening behavior considered in all these models, it is evaluated whether the models use more than one nonlinear variable to more precisely reproduce this behavior. Due to that in some forming processes dominated by cyclic loading, higher modeling accuracy is required for the nonlinear hardening behavior after loading reverse, it is no more adequate to reproduce the Bauschinger effect by using a simple nonlinear function. According to their main mathematical implementation methods summarized in Section 3.1, these models can be roughly classified into four categories: kinematic models (denoted by “k”), distortional models (denoted by “d”), varied flow stress models (denoted by “vfs”), and pure distortional models without kinematic components. Ideally, to obtain a comprehensive and flexible description of anisotropic hardening, each decomposed hardening behavior should be reproduced by the specialized variables in the model, which inevitably results in the introduction of numerous parameters. In fact, some decomposed hardening behaviors may not be completely decoupled physically. For example, as is supposed by Peeters [84], permanent hardening behaviors such as permanent softening are to some extent triggered by the previous transient hardening behaviors such as the Bauschinger effect and hardening stagnation. Therefore, a trade-off between the accuracy and simplicity of the hardening model is often made in practical use. In Table 2, the number of parameters required in each hardening model is also summarized for comparison. By referring to Table 1 and Table 2 successively, the appropriate hardening models can be preliminary selected according to the target material and the process related loading modes. Details of these modeling methods will be introduced in the remaining part of this section.

4.1. Reverse Loading: The Bauschinger Effect

The Bauschinger effect is one of the most universal anisotropic hardening behaviors. Under this effect, hardening asymmetry is observed between forward and reverse loading. The monolithic translation (kinematic hardening) or local contraction (distortional hardening) of the yield surface is used to describe this effect.

The framework of kinematic hardening is presented in Equation (4). In the 1950s, researchers proposed the concept of yield surface translation [85]. Since then, different formulations have been proposed for the back stress $\mathit{\alpha }$ to represent the Bauschinger effect. Prager [60] and Ziegler [86] proposed two simplified linear evolution functions, which can be uniformly expressed as follows:

$$d\mathit{\alpha }=k\mathit{N}d\overline{\epsilon }$$

(21)

In Prager’s formulation, the translation direction of the yield surface is collinear with the normal of the current yield point (${\mathit{N}}_{\mathit{\epsilon }}$), while in Ziegler’s formulation, this direction is along the vector connecting the yield surface center to the current stress point (${\mathit{N}}_{\mathit{s}-\mathit{\alpha }}$). Ziegler’s expression is convenient for solving the plane stress problem, which is a common assumption in sheet forming. However, in linear formulations, the evolution rate remains constant; therefore, they are unable to reflect the transiency of the Bauschinger effect and overestimate the proportion of kinematic hardening, especially over a large strain range. Multilinear and nonlinear models have been proposed for overcoming this problem. Mróz [87] proposed a multi-surface model in which each surface represents a constant hardening rate. In the loading process, these yield surfaces work in sequence so that the nonlinearity of the Bauschinger effect can be captured. Armstrong and Frederick [70] proposed a nonlinear saturated evolution function, which is expressed in Equation (17). By adjusting the values of $X_{sat}$ and $k$ in Equation (17), the transient behavior can be well reproduced. Because of its concise form and relatively accurate description capability, AF 66 model has been widely used and is treated as a basic component in different anisotropic hardening models (such as the model with the kinematic components as listed in Table 2).

For reverse loading-dominated deformation behaviors, especially those under cyclic loading, the high nonlinearity of the Bauschinger effect must be carefully described. In the hardening models used to study the aforementioned behaviors, multiple back stresses in the form of Equation (17) (denoted as “AF back stress”) are combined. Chaboche proposed a series of kinematic hardening models for analyzing cyclic plasticity, especially for reproducing the ratcheting effect [75,88,89,90]. In the Chaboche 86 model, several decomposed AF back stresses are superposed as follows to reproduce the high nonlinearity (corresponding to three-phased nonlinearity) of the Bauschinger effect:

$$d\mathit{\alpha }={\displaystyle \sum }_{i=1}^{n}d{\mathit{\alpha }}_i$$

(22)

where $n$ denotes the number of decomposed back stresses. Usually, $n$ is equal to 3 or $4$. When $n=3$ (Figure 7), three back stresses respectively reproduce (i) the initial high hardening rate at the onset of re-yielding, (ii) the nonlinearity of the transient segment, and (iii) the near-linear or constant hardening rate in larger strain range. Therefore, expression of the third back stress often reduces to the linear Prager model (Equation (21)). For a more accurate simulation of the ratcheting effect, a threshold value or step function is introduced into the expression for the forth decomposed back stresses in the hardening model [89,91]. A detailed discussion on cyclical hardening models is presented in Chaboche’s review [92], which is not discussed in detail in this paper. Lee et al. [93] replaced the fixed material constants in the Chaboche 86 model with condition functions to incorporate the anisotropy of kinematic hardening. In the YU 02 model, two AF back stresses work in sequence, and one of them acts as a varied saturation limit ($X_{sat }$ in Equation (17)) for the actual back stress during the deformation process; namely, a nested constraint relationship of two back stresses exists in YU 02. Figure 8a presents a schematic of the kinematic hardening process described by the YU 02 model. An outer bounding surface restrains the evolution of the inner yield surface. At the onset of yield, the inner yield surface moves at a relatively high speed (red curve in Figure 8b). After it touches the outer bounding surface, the hardening rate is suppressed to be consistent with that of the outer surface; thus, a two-segment hardening feature can be reproduced using the YU 02 model. The description of phased hardening features is more independent in the YU 02 model than in Chaboche’s methods; however, this model does not include a description of the constant rate over a large strain range.

As mentioned in Section 3, unlike isotropic hardening, kinematic hardening inevitably introduces a softening effect in additional loading directions. Barlat et al. [22] solved this problem by introducing a fluctuating term ${\mathsf{\psi }}_h$ to the yield function in the HAH 11 model (Equation (14)). This fluctuating term acts as an asymptotic line in the π-plane to compress the half-facet opposite the microstructure tensor $\widehat{\mathit{h}}$ (Figure 4a) so that the description of the Bauschinger effect does not influence the hardening behaviors occurring in the direction along or orthogonal to $\widehat{\mathit{h}}$. In contrast to the situation under the monolithic translation of the yield surface, under local distortion, two independent expressions are required to reproduce the pre-accumulation and attenuation of the Bauschinger effect, respectively. The term ${\psi }_h$ in Equation (14) is expressed as follows:

$${\mathsf{\psi }}_h^q=f_{-}^q\langle -2{\mathit{s}}_p\rangle {}^q+f_{+}^{q}2{\mathit{s}}_p{}^q$$

(23)

where 〈〉 represents MacCauley brackets. When the internal number is positive, its value will be reserved; otherwise, the value is set as 0. The terms $f_{-}$ and $f_{+}$ denote the intermediate variables controlling the reverse loading’s influence on two facets opposite and coincident to $\widehat{\mathit{h}}$, respectively. The evolutions of $f_{-}$ and $f_{+}$ are controlled by the potential $g_{-}$ (controls the preaccumulation of the Bauschinger effect) and active $g_{+}$ (controls the attenuation of the Bauschinger effect), respectively.

$$d{g}_{-}=k_2\left(k_3\frac{{\mathsf{\sigma }}_f\left(0\right)}{{\mathsf{\sigma }}_f\left(\overline{\epsilon }\right)}-g_{-}\right)d\overline{\epsilon }$$

(24)

$$d{g}_{+}=k_1\frac{1-g_{+}}{g_{+}}d\overline{\epsilon }$$

(25)

Geometrically, $g$ represents the reverse loading-induced distortion ratio, that is, the ratio of the radii of the normalized distorted and original isotropic yield surfaces. Moreover, unlike the back stress $\mathit{\alpha }$ used in kinematic hardening models, the variable ${\psi }_h$ equals to 0 in monotonic loading. Thus, in HAH models, the experimentally measured reference flow stress curve corresponds to the flow stress curve of the isotropic hardening component, which is convenient for parameter identification.

However, as concluded by Holmedal [74], several defects exist in the initial version of the HAH model. For instance, in numerical calculation, a singularity occurs at the conjunction of two distorted facets (i.e., the orthogonal loading stress state with $\cos \mathsf{\chi }=0$). In HAH 20, this drawback is overcome by introducing additional parameters to improve the continuity of the transition region. Moreover, in specific materials with a strong Bauschinger effect, the distorted yield surface does not contain the origin in the space of σ-τ (Figure 3b). This situation cannot be modeled using Equation (14). In Holmedal 19 model, two back stress components are used to represent the role of $\widehat{\mathit{h}}$. The concept of kinematic hardening is introduced into this distortion hardening model, so that the location of the distorted yield surface is controlled by the back stress but not a fixed origin such as the HAH models. In this manner, the Bauschinger effect-induced distortion behaviors can be reproduced with high flexibility.

4.2. Reverse Loading: Hardening Stagnation

Hardening stagnation is represented by a transient plateau in the stress–strain curve, which occurs near the end of the Bauschinger effect. In common hardening models, two strategies are adopted to reproduce this behavior. In the direct method, the hardening rate is temporarily reduced in the segment where hardening stagnation occurs. In the indirect method, a finite stress overshoot is inserted into the model after loading reversal. This overshoot diminishes gradually as deformation continues and forms a plateau in the stress–strain curve.

The YU 02, TH 98, and Wang 08 models adopt the direct method. In the YU 02 model, the hardening level after loading reversal is restricted by a newly defined non-isotropic hardening (non-IH) surface. The role of this surface is illustrated in Figure 9. During plastic deformation, the center of the bounding surface does not cross the boundary of the inner non-IH surface. Mathematically, the evolution of the isotropic hardening component does not begin until it reaches the non-IH surface; thus, the hardening of the bounding surface is restricted in the early period after loading reversal. After the initial high hardening rate of the Bauschinger effect diminishes, the hardening behavior is completely controlled by the bounding surface. Thus, as shown in Figure 9, hardening saturation occurs after the yield surface touches the bounding surface and before the center of the bounding surface reaches the non-IH surface. Moreover, the introduction of a stagnation segment forces the subsequent hardening curve to fall behind the isotropic hardening curve so that permanent softening is also reproduced. The TH 98 and Wang 08 reproduce the anisotropic hardening behaviors by superimposing the norm of a fourth-order tensor on the reference flow stress ${\mathsf{\sigma }}_f\left(\overline{\epsilon }\right)$. The hardening saturation varies with the parallel component ${\mathit{S}}_p$ (discussed in Section 3.1.2).

The indirect method is described in the following text. In the MHH 15 and QHH 18 models, a stress overshoot $S_r$ is superimposed on the reference flow stress curve ${\sigma }_f\left(\overline{\epsilon }\right)$ and takes effect after loading reverses. The evolution function of $S_r$ follows the format of Equation (17) and is expressed as follows:

$$d{S}_r=k_{Sr}\left(S_{rsat}\langle -\Vert \mathit{P}\Vert \cos \mathsf{\chi }\rangle -S_r\right)d\overline{\epsilon }$$

(26)

where $\mathit{P}$ is the delayed pointer (mentioned in Section 3.2.2) and $S_{rsat}\langle -\Vert \mathit{P}\Vert \cos \mathsf{\chi }\rangle$ represents the saturation values of $S_r$. Figure 10 illustrates the evolution trends of $S_r$ and its saturation value. At the onset of loading reversal, $S_{rsat}\langle -\Vert \mathit{P}\Vert \cos \mathsf{\chi }\rangle$ reaches its maximum value ($S_{rsat}$) and then decreases to 0 gradually as deformation continues, causing the curve of the overshoot $S_r$ to first rise and then fall. After appropriate parameters are determined, a plateau is formed in the stress–strain curve (the segment traced in red in Figure 10). Moreover, the modification method used in MHH 15 and QHH 18 models can simultaneously reproduce the high initial hardening rate of the Bauschinger effect (the segment traced in blue in Figure 10). The Qin 17 model is an extension of the original HAH model that can describe hardening stagnation. The Qin 17 model is improved by adjusting the variable that represents the Bauschinger effect. As presented in Equation (25), in the HAH model, $g_{+}$ denotes the influence of the active Bauschinger effect on the yield surface contraction ratio. $g_{+}$ is updated by considering the influence of hardening stagnation. The updated variable is denoted as ${g_{+}}^{\prime }$, which is expressed as follows:

$${g_{+}}^{\prime }=g_{+}+g_{\mathit{s}}\left(1-g_{+}\right)$$

(27)

Because $g_{+}\le 1$, $g_{\mathit{s}}$ introduces the hardening stagnation’s alleviation effect on the Bauschinger effect after the loading reversals. The evolution of $g_{\mathit{s}}$ follows the basic function expressed in Equation (17).

4.3. Reverse or Orthogonal Loading: Permanent Softening

Permanent softening can be observed under reverse or orthogonal loading, depending on material properties. Because their physical mechanism and softening severity are different, in most instances, reverse permanent softening and orthogonal permanent softening are described independently. In some models, permanent softening is assumed to be caused by the influence of former transient hardening behaviors, namely the Bauschinger effect and hardening stagnation; therefore, in these models, permanent softening is indirectly described by the variables which control the aforementioned transient behaviors. Certain specialized variables are used in some other models for independently capturing the permanent softening effect.

In the Chaboche 86 and YU 02 models, the descriptions of reverse permanent softening are included in their description of permanent softening and hardening stagnation. In the Chaboche 86 model, the third decomposed back stress evolves according to a linear law (Equation (21)) without a saturation value; thus, after loading reverses, the global stress level always lags the isotropic hardening level. This persistent lag is mathematically consistent with the effect of permanent softening. However, since Equation (21) merely contains one parameter, loading history’s influence on the permanent softening cannot be well reproduced. For another, the YU 02 model records the permanent softening behavior mainly through the stress hysteresis triggered by hardening stagnation.

The MHH 15 (and QHH 18) model also provides indirect descriptions of permanent softening behavior. This model uses two specialized variables to reproduce the permanent softening behavior induced by reverse and orthogonal loading. The aforementioned model is constructed by modifying the reference flow stress curve (Equation (15)). Hardening is transiently suppressed immediately after the loading path changes so that the corresponding flow stress curve lags that obtained under monotonic loading. The evolution functions of hardening reduction in the flow stress after reverse ($R_{\mathit{s}r}$) and orthogonal ($R_{\mathit{s}o}$) loading path changes are expressed as follows, respectively:

$$\mathrm{d}{R}_{\mathit{s}r}=-k_{\mathit{s}r}\langle -\Vert \mathit{P}\Vert \cos \mathsf{\chi }\rangle d\overline{\epsilon }$$

(28)

$$\mathrm{d}{R}_{\mathit{s}o}=-k_{\mathit{s}o}\Vert \mathit{P}\Vert \sin \mathsf{\chi }d\overline{\epsilon }$$

(29)

HAH models provide relatively independent descriptions of reverse permanent hardening. The original HAH model reproduces this behavior by modifying the variable that represents the Bauschinger effect. Equation (25) presents the attenuation law for the Bauschinger effect. In this equation, the constant “1” in the numerator indicates that the yield surface regains an isotropic nature as long as sufficient proportional deformation occurs. Accordingly, to describe permanent softening, a variable less than 1 is used to substitute the constant “1”. The variable substituting “1” is denoted as $g_c$ in this paper. The evolution formulations for the modified $g_{+}$ and $g_c$ values are expressed as follows:

$$d{g}_{+}=k_1\frac{g_c-g_{+}}{g_{+}}d\overline{\epsilon }$$

(30)

$$d{g}_c=\left(k_c-g_c\right)d\overline{\epsilon }$$

(31)

where the coefficient $k_c\le 1$. This parameter controls the ultimate softening rate. However, an inadequacy exists in this evolution function, because $g_c$ continues increasing once plastic deformation occurs, irrespective of the loading path. Thus, the softening level is overestimated near orthogonal loading path changes. Lee et al. [94] fixed this inadequacy by modifying Equation (30) as follows:

$$d{g}_{+}=k_1\frac{g_c+\left(1-g_c\right){\sin }^2\mathsf{\chi }-g_{+}}{g_{+}}d\overline{\epsilon }$$

(32)

According to the aforementioned formulation, permanent softening does not occur when the orthogonal loading path changes. The HAH 20 model solves the aforementioned inadequacy by modifying the reference flow stress curve. After a series of intermediate variables are introduced into the aforementioned model, the softening level between reverse and orthogonal loading can be freely regulated. Furthermore, Barlat observed that the stress-differential effect, which is often ignored in the yield criteria for a metal, has the same magnitude as the permanent softening effect under reverse loading. Therefore, to acquire a more reliable prediction of the permanent softening-related hardening behaviors, the effect of hydrostatic pressure is incorporated into the HAH 20 model.

4.4. Orthogonal Hardening

In some materials, transient hardening occurs immediately after orthogonal loading path changes. To reproduce this behavior, early studies adjusted the relevant anisotropic fourth-order tensor, whereas later studies adopted more explicit modeling methods.

Fourth-order tensor variables can be used to simulate the distortion of the yield surface (denoted as $\mathit{A}$) or incorporate transient features into scalar evolution equations (such as the flow stress curve) to represent some deformation-induced hardening behaviors (denoted as $\mathit{S}$). The second methodology is adopted in the TH 98 and Wang 08 hardening models. In these models, a variable saturation value of back stress is obtained as follows:

$${\alpha }_{sat}={\alpha }_{sat0}+\left(1-f\right)\sqrt{{\Vert {\mathit{S}}_p\Vert }^2+m\Vert {\mathit{S}}_o\Vert {}^2}$$

(33)

where ${\alpha }_{sat0}$ denotes the initial saturation value of back stress; $f$ and $m$ are two constants; and ${\mathit{S}}_p$ and ${\mathit{S}}_o$ represent the parallel and orthogonal components of the tensor $\mathit{S}$, respectively. These components dominate the effects of reverse and orthogonal loading on increasing the saturation value of the back stress (${\alpha }_{sat}$). The detailed evolution function of $\mathit{S}$ is presented in Wang’ paper [68]. When $m>1$, orthogonal loading increases ${\alpha }_{sat}$ to a more considerable extent; thus, the orthogonal hardening behaviors can be indirectly reproduced. Moreover, in the Norman 10 model, an elongation of the yield surface perpendicular to the current loading direction is modeled by modifying the fourth-order anisotropic tensor $\mathit{A}$. The yield function for such surface is expressed in Equation (7), and the evolution functions of the components of $\mathit{A}$ are expressed in the format of incremental iterations. First, an orthogonal decomposition is conducted according to Equations (8) and (9) to extract the orthogonal component ${\mathit{A}}_o$ and parallel component ${\mathit{A}}_p$ from the initial tensor ${\mathit{A}}_{old}$. Subsequently, for a given strain increment, the corresponding increments of ${\mathit{A}}_o$ and ${\mathit{A}}_p$ are calculated as follows:

$$d{\mathit{A}}_o=k_o\left(A_{osat}\left({\mathit{I}}_{dev}-\mathit{N}\otimes \mathit{N}\right)-{\mathit{A}}_o\right)d\overline{\mathit{\epsilon }}$$

(34)

$$d{\mathit{A}}_p=k_p\left(A_{psat}\mathit{N}\otimes \mathit{N}-{\mathit{A}}_p\right)d\overline{\mathit{\epsilon }}$$

(35)

where $A_{osat}$ and $A_{psat}$ denote the saturation norms of ${\mathit{A}}_o$ and ${\mathit{A}}_p$, respectively. After loading path changes, the orthogonal hardening effect accumulated in the former path gradually weakens; thus, $A_{psat}$ is often set as 0. Finally, the new tensor ${\mathit{A}}_{new}$ is acquired as follows:

$${\mathit{A}}_{new}={\mathit{A}}_{old}+d\mathit{A}={\mathit{A}}_{old}+d{\mathit{A}}_o+d{\mathit{A}}_p$$

(36)

After the aforementioned steps are completed, the value of ${\mathit{A}}_{new}$ is assigned to ${\mathit{A}}_{old}$ to start the next loop. The evolution of the yield surface corresponding to the Norman 10 model is displayed in Figure 11.

The aforementioned high-order tensors are rarely used in newly proposed hardening models. Most of these models reproduce the transiency of orthogonal hardening by employing the Schmitt factor. The Schmitt angle, which initially has a finite value, gradually decreases to 0 when proportional loading occurs after loading path changes. Thus, this variable inherently possesses a transient feature. MHH 15 and QHH 18 models’ description of orthogonal hardening resembles their description of hardening saturation. The additional hardening that occurs after orthogonal path changes is superposed on the flow stress curve. The evolution functions used in the MHH 15 and QHH 18 models are respectively listed as follows:

$$d{S}_o=k_{So}\left(S_{osat}\Vert \mathit{P}\Vert \sin \mathsf{\chi }-{\mathit{s}}_o\right)d\overline{\epsilon }$$

(37)

$$d{\mathit{s}}_o=k_{\mathit{s}o}\left(S_{osat}{\sin }^n\mathsf{\chi }-S_o\right)d\overline{\epsilon }$$

(38)

HAH models provide descriptions of the orthogonal hardening behavior by modifying the stress tensor. The HAH 13 model assumes that the influence of orthogonal softening on the yield surface is isotropic. Thus, the reference flow curve ${\sigma }_f\left(\overline{\epsilon }\right)$ is magnified by multiplying it with a scaling factor $g_l$. For the sake of the comparison, an equivalent form of the method adopted by HAH 13 model is expressed; namely, the stress $\mathit{s}$ is modified as a whole and redefined as ${\mathit{s}}_d$:

$${\mathit{s}}_d=\frac{1}{g_l}\mathit{s}$$

(39)

where $g_l$ represents the homogenous expansion rate of the yield surface after orthogonal loading. In the HAH 14 model, the scaling factor $g_l$ merely acts on the orthogonal component ${\mathit{s}}_o$; thus, expansion only occurs in the orthogonal direction as follows:

$${\mathit{s}}_d={\mathit{s}}_p+\frac{1}{g_l}{\mathit{s}}_o$$

(40)

The evolution of the yield surface in the HAH 14 model under an abrupt loading path change is illustrated in Figure 12. An orthogonal expansion occurs at the onset of loading path changes and gradually diminishes as persistent proportional deformation occurs. The HAH 20 model uses a flexible modification method that can be expressed as follows:

$${\mathit{s}}_d=\frac{1}{{\xi }_l\left(g_l-1\right)+1}{\mathit{s}}_p+\frac{1}{g_l}{\mathit{s}}_o$$

(41)

where ${\xi }_l$ is a coefficient between 0 and 1 for adjusting the expansion in different loading directions. When ${\xi }_l=1,$ Equation (41) is transformed into Equation (39), whereas when ${\xi }_l=0,$ Equation (41) is transformed into Equation (40). The evolution mode of $g_l$ is based on the dislocation-based hardening function proposed by Rauch et al. [17] and follows the universal format presented in Equation (17). This mode can be expressed as follows:

$$d{g}_l=k_l\left(\frac{{\mathsf{\sigma }}_f\left(\overline{\epsilon }\right)-{\mathsf{\sigma }}_f\left(0\right)}{{\mathsf{\sigma }}_f\left(\overline{\epsilon }\right)}\left(\sqrt{l+\left(1-l\right){\cos }^2\mathsf{\chi }}-1\right)+1-g_l\right)d\overline{\epsilon }$$

(42)

where $l$ controls the degree of orthogonal hardening. The evolution of the yield surface in two-step loading is depicted in Figure 12. Orthogonal expansion occurs rapidly at the onset of loading path changes and gradually diminishes after the occurrence of deformation.

4.5. Orthogonal Softening

Orthogonal softening is observed in some materials. This hardening behavior is naturally incorporated into kinematic hardening models by considering the monolithic translation of the yield surface, which represents the Bauschinger effect. Therefore, models that contain back stresses rarely include variables that are directly related to orthogonal softening. However, the TH 98 model includes a specific description of orthogonal softening, which is obtained by constructing a recession function for the orthogonal component of the fourth-order tensor ${\mathit{S}}_o$ that is used in Equation (33).

Nevertheless, pure distortional hardening models (e.g., HAH models) do not specifically model the orthogonal softening behavior. In HAH models, the basic format presented in Equation (14) is adopted for the yield function to consider the orthogonal softening effect. The yield function can be expressed as follows:

$$\psi ={\left({\psi }_0^p\left({\mathit{s}}_d\right)+{\psi }_0^p\left({\mathit{s}}_x\right)\right)}^{\frac{1}{p}}$$

(43)

where ${\mathit{s}}_x$ is the result obtained by linearly transforming the stress orthogonal component ${\mathit{s}}_o$ by using $g_o$ (i.e., ${\mathit{s}}_x={\mathit{s}}_x\left({\mathit{s}}_o,g_o\right)$) and $p$ denotes a constant exponent, different from $q$ used in Equation (14). The term ${\psi }_0\left({\mathit{s}}_x\right)$ represents two lines perpendicular to $\widehat{\mathit{h}}$ in π-plane. These two lines are symmetrically distributed on either side of the origin, and the distance from the origin is controlled by $g_o$. The function ${\psi }_0\left({\mathit{s}}_x\right)$ compresses the yield surface in the orthogonal loading direction so that the softening behavior can be reproduced. For analyzing the transiency of the orthogonal softening behavior, two different evolution functions are used in the HAH 14 and HAH 20 models respectively. Because the evolution modes of these models are similar to the methods describing other transient hardening behaviors in HAH models, these modes are not reintroduced here.

4.6. “Egg-Shaped” Distortion of the Yield Surface

For enhancing the reliability of an anisotropic hardening model, the intermediate hardening behaviors between forward, reverse, and orthogonal loading must be considered (e.g., the “egg-shaped” distortion of the yield surface). Although the aforementioned modeling methods indirectly consider the intermediate distortion to some degree, including specialized variables for representing this anisotropic behavior is still necessary, especially for processes relevant to complex loading conditions.

The HAH 20 model reproduces the intermediate distortion behaviors in a simplified manner. In this model, no additional variables are introduced; however, a series of constant coefficients are included in those variables that describe transient hardening behaviors to reflect their influences on the intermediate distortion. For example, in Equation (41), ${\xi }_l$ is used to describe the intermediate distortion related to orthogonal hardening. Moreover, the exponents $q$ and $p$ are used to adjust the curvature of the yield surface in Equations (14) and (43), respectively. The aforementioned coefficients considerably improve the flexibility of hardening models; however, the use of an excessive number of such coefficients renders these models complicated.

Conventional description methods for intermediate distortion are described in the following text. In numerous hardening models, the distortions of the yield surface are assumed to be collinear to the loading direction variable ${\mathit{X}}_{\mathit{d}}$. The anisotropic tensor $\mathit{A}$ (Equation (7)) or distorted stress tensor ${\mathit{s}}_d$ (Equation (10)) can be used to reflect the directionality of deformation.

Early hardening models, such as the Baltov 65 model, reproduce the torsion of the yield surface by constructing a reasonable fourth-order anisotropic tensor $A^{\prime }$. This tensor is decomposed into two parts as follows:

$$A_{ijkl}^{\prime }=I_{ijkl}+A_{d ijkl}$$

(44)

where $I_{ijkl}$ is the initial isotropic tensor, and the distortion of the yield surface is dominated by $A_{d ijkl}$:

$$A_{d ijkl}=A_0{\epsilon }_{ij}{\epsilon }_{kl}$$

(45)

where $A_0$ controls the shape of the yield surface. For negative or positive values of $A_0$, the yield surface heterogeneously expands or contracts, respectively. However, the distortion is consistent in the loading and reverse loading directions; thus, the yield surface remains a symmetrical ellipse throughout the deformation process. Early relevant studies constructed high-order anisotropy tensors to model the “egg-shaped” evolution of the yield surface. Voyiadjis and Foroozesh [78] incorporated an additional loading direction-related tensor into the framework of Baltov 65 (Equation (44)) and successfully reproduced the asymmetric distortion of the yield surface. Three high-order tensors, namely a second-order tensor, a fourth-order tensor, and a sixth-order tensor, are used in the model proposed by Rees [95]. The aforementioned models inevitably use numerous coefficients, which heavily limits their practicability.

Newly proposed models describe the distorted yield surface based on the Schmitt factor, which results in a small number of model parameters being required. The asymmetry problem of the Baltov 65 model is overcome in the Rokhgireh 13 model by modifying Equation (45) as follows by using the Schmitt factor:

$$A_{d ijkl}=A_0\cos \mathsf{\chi }{\epsilon }_{ij}{\epsilon }_{kl}$$

(46)

When $A_0$ is determined to be a negative value, the “egg-shaped” distortion of the yield surface forms naturally. The formulation presented in Equation (46) is further enhanced in the Rokhgrieh 17 model by adding in four decomposed back stresses for the superior prediction of multiaxial ratchetting. However, the Baltov 65, Rokhgrieh 13, and Rokhgrieh 17 models are only suitable for materials without initial anisotropy. The Feigenbaum 07 and Pietryga 12 models are similar to the aforementioned models but use a distortion-relevant component in the original isotropic or anisotropic tensor $\mathit{A}$. The distorted tensor ${\mathit{A}}^{\prime }$ in the Feigenbaum 07 model is expressed as follows:

$${\mathit{A}}^{\prime }=\mathit{A}+\left(\Vert \mathit{\alpha }\Vert \cos \mathsf{\chi }\right){\mathit{A}}_{\mathit{h}}$$

(47)

Here, the expression of $\cos \mathsf{\chi }$ is based on the definition of ${\mathit{N}}_{\mathit{s}-\mathit{\alpha }}$, in the parentheses of Equation (2). During the deformation process, the increment in the distortion-relevant tensor ${\mathit{A}}_{\mathit{h}}$ remains negative. Thus, $\left(\Vert \mathit{\alpha }\Vert \cos \mathsf{\chi }\right){\mathit{A}}_{\mathit{h}} \le 0$ as long as $\cos \mathsf{\chi } \ge 0$. Consequently, the front of the yield surface is elongated and forms a corner. If $\cos \mathsf{\chi }<0$, the yield surface is compressed, and a flat zone is developed at its rear (Figure 13a). However, as α is applied to describe the Bauschinger effect and “egg-shaped distortion” simultaneously, some flaws are inevitably introduced by this coupled method. For example, in reverse loading, once the back stress α evolves to 0, the distortion of the yield surface accordingly diminishes, which is unreasonable. The Pietryga 12 model uses a relatively simplified description of the aforementioned phenomena on the basis of Noman’s framework, which is expressed in Equations (34)–(36). The description in the Pietryga 12 model is expressed as follows:

$${\mathit{A}}^{\prime }=\mathit{A}+\left(-c_d\cos \mathsf{\chi }\right){\mathit{A}}_{\mathit{p}}$$

(48)

where $\mathit{A}$ denotes the tensor that considers orthogonal hardening in Noman’s model and ${\mathit{A}}_{\mathit{p}}$ denotes the parallel component. Equation (48) has similar influences on the yield surface as does Equation (47). Ortiz and Popov [96] conducted distortional hardening modeling without using anisotropy tensors. In their model, a Schmitt angle-related Fourier series is added to the flow stress curve. To simplify their hardening model, Feigenbaum and Dafalias [79] eliminated the fourth-order tensor from it. Their model describes the directional distortion merely by using a Schmitt factor-related polynomial. A certain class of hardening models reproduce the distorted yield surface according to the geometrical concept that was proposed by Phillips and Weng [97] and later extended by Eisenberg and Yen [98]. In the model developed by Yeh and Lin [99], the yield surface is divided into two parts, and the forward and rear curvatures of the yield surface are controlled using different ellipse functions. However, the aforementioned class of models has similar flaws. An additional translation of the yield surface center is introduced by the description of the intermediate distortion, which inevitably affects the description of the Bauschinger effect according to the back stress.

Using a distorted stress tensor is a simple and flexible skill to describe “egg-shaped” distortion. By using a distorted stress tensor, the descriptions of the Bauschinger effect and “egg-shaped” distortion can be completely decoupled. The François 01 model substitutes the original stress tensor $\mathit{s}$ (Equation (10)) with a distorted stress tensor ${\mathit{s}}_{\mathit{d}}$, which is expressed as follows:

$${\mathit{s}}_{\mathit{d}}=\mathit{s}+\mathit{d}=\mathit{s}+\frac{{\Vert {\mathit{s}}_{\mathit{o}}\Vert }^2}{2{\mathsf{\sigma }}_f\left(\overline{\epsilon }\right)}\frac{\mathit{\alpha }}{{\alpha }_{sat}}$$

(49)

where $\mathit{d}$ denotes the distortion component of the modified stress tensor ${\mathit{s}}_{\mathit{d}}$, and ${\alpha }_{sat}$ denotes the kinematic limit. Equation (49) is concise because the distortion degree is completely dominated by the evolution of $\mathit{\alpha }$ in the deformation process. Moreover, for a given deformation, the distortion component $\Vert \mathit{d}\Vert$ is supposed to be parabolic to the orthogonal component of stress $\Vert {\mathit{s}}_{\mathit{o}}\Vert$. Thus, as shown in Figure 13b, the distorted yield surface does not deviate from the yield point determined using a kinematic hardening model under reverse loading. Badreddine et al. [100] added a new term orthogonal to the kinematic component in the distorted stress for better reproducing the yield surface evolution examined by Khan et al. [29,30]. In the Zhang 20 model, Equation (49) is expanded to make the expression of distortion compatible with nonquadratic yield functions.

Finally, it is necessary to note that the aforementioned hardening models are specially constructed for steel or aluminum. As for other materials, for example magnesium alloy, the hardening behavior might be different. Thus, the distortion direction of the yield surface should not be always restricted in modeling. In addition to the fluctuating component ${\psi }_h$ in HAH models (Equation (14)), He 18 model includes two additional fluctuating components, namely ${\psi }_{t1}$ and ${\psi }_{t2}$, in the yield function for magnesium alloy. By using plane rotational tensors, yield surface distortion can be achieved along arbitrary directions.

5. Conclusions and Recommendations

This paper reviews the phenomenological modeling of deformation-induced anisotropic hardening behaviors. An overview of the anisotropic hardening behaviors is laid out according to their classification, the influence from the materials and forming processes, and the relevant experiments. The phenomenological modeling of the anisotropic hardening behaviors is reviewed from two aspects. First, the core modeling idea of anisotropic hardening behaviors, and typical mathematical description methods are extracted from different models. Second, the modeling implementations for each decomposed anisotropic hardening behavior (i.e., reverse loading hardening behaviors, orthogonal loading hardening behaviors, and intermediate distortion of the yield surface) are introduced. This paper tries to provide an instruction for the researchers to model for the anisotropic hardening behaviors according to their practical requirements. The conclusions of this review and the recommendations for the future development of the modeling for anisotropic hardening behaviors are listed as follows:

• A more comprehensive mathematical description of loading path changes: Both the decomposition and the modeling of hardening behaviors relies on the mathematical description of loading path changes. In the models reviewed in this paper, such descriptions are based on a scalar, which is typically represented by the Schmitt angle. Due to the conciseness of this scalar variable, its capability to measure the modes of loading path changes is restricted. The contraction operation of two tensors (Equation (3)) omits much information of the stress states before and after loading path changes, which might also closely relate to the hardening behaviors. For example, despite with the similar Schmitt angles, the hardening behaviors observed in the tension-tension tests and the rolling-tension tests are not completely the same [20]. Thus, a more comprehensive mathematical description of loading path changes to reserve more details of the loading path changes is needed.
• The trade-off between model completeness and simplicity according to the practical requirements: As mentioned in Section 3.2.1, the sufficient description of each decomposed nonlinear hardening behavior requires an independent state variable comprising at least two parameters. Thus, a detailed description of each anisotropic hardening behavior may result in the excessive parameters in a hardening model. For one thing, to ensure their completeness, the hardening models must meet the following requirements: (i) The descriptions of different hardening behaviors should be decoupled with each other, compatible to be insert into different models, and flexible to reproduce different details of the hardening behaviors (e.g., whether the phenomenon that the origin of the stress space lie outside the yield surface can be reproduced); (ii) global optimization should be avoided in parameter determination because the effects of different parameters may diminish over time; and (iii) more microstructure evolution mechanism from physics-based models should be incorporated (e.g., the description of orthogonal hardening in HAH 14 (Equation (42)) is inspired by the dislocation-based hardening function in Rauch 11). For another, a hardening model’s competence and accuracy can be sacrificed to some extent for easy use. Simplification can be conducted by describing multiple hardening behaviors with a single variable (e.g., the description of the Bauschinger effect and yield surface distortion in the François 01 and Feigenbaum 07 model); however, such simplification may introduce some unexpected problems. Thus, hardening modeling should be conducted in close accordance with the specified materials and forming processes.
• Incorporating the characteristics of the initial yield into the modeling for hardening behaviors: For one thing, since the hardening behaviors are modeled based on the yield functions, the compatibility between yield functions and hardening behaviors must be considered (e.g., HAH models can use yield functions of any order; while, François 01 model is not compatible with nonquadratic yield functions until it is expanded in Zhang 20). For another, the precision of a hardening model is improved when it is used together with a proper yield criterion. For example, as is discussed in HAH 20, the yield relevant stress-differential effect is in the same order of magnitude as the Bauschinger effect. Thus, to ensure the accuracy of hardening model, the influences of yield behaviors such as the stress-differential should be considered simultaneously. Furthermore, although anisotropic yield functions are specified to reproduce the initial anisotropy, using the initial value of the state variables to reproduce this feature is regarded as a more flexible approach. However, how to reflect the influences of pre-rolling on the current plasticity by using the initial values of state variables in the hardening models remains unclear [77].
• Modeling the anisotropic hardening behaviors for the materials of more kinds: Up to now, most of the anisotropic hardening models are constructed for conventional materials such as steel and aluminum. Among the models reviewed in this paper, only He 18 is specially applied for magnesium alloy, a HCP material. As can be seen in He’s work, the hardening behaviors of HCP materials and corresponding modeling methods are quite different. With the further requirement for the performance of metal materials, an increasing amount of advanced materials such as magnesium and titanium are used in sheet metal forming. Thus, anisotropic hardening models for such materials are in great need.