A Homogeneous Anisotropic Hardening Model in Plane Stress State for Sheet Metal under Nonlinear Loading Paths

In sheet metal forming, the material is usually subjected to a complex nonlinear loading process, and the anisotropic hardening behavior of the material must be considered in order to accurately predict the deformation of the sheet. In recent years, the homogeneous anisotropic hardening (HAH) model has been applied in the simulation of sheet metal forming. However, the existing HAH model is established in the second-order stress deviator space, which makes the calculation complicated and costly, especially for a plane stress problem such as sheet metal forming. In an attempt to reduce the computational cost, an HAH model in plane stress state is proposed, and called the HAH-2d model in this paper. In the HAH-2d model, both the stress vector and microstructure vector contain only three in-plane components, so the calculation is significantly simplified. The characteristics of the model under typical nonlinear loading paths are analyzed. Additionally, the feasibility of the model is verified by the stress–strain responses of DP780 and EDDQ steel sheets under different two-step uniaxial tension tests. The results show that the HAH-2d model can reasonably reflect the Bauschinger effect and the permanent softening effect in reverse loading, and the latent hardening effect in cross loading, while the predictive accuracy for cross-loading softening remains to be improved. In the future, the HAH-2d model can be further modified to describe more anisotropic hardening behaviors and applied to numerical simulations.


Introduction
With increasing requirements for lightweight materials in the aerospace, aviation, and automotive industries, integrated thin-walled components with complex shapes are finding broader applications [1,2]. During the forming of an integrated component, the material often undergoes complex nonlinear loading paths, even including multi-step pre-deformation. The deformation behavior of sheet metal often exhibits strong path dependence [3][4][5] and pronounced anisotropic hardening behaviors in nonlinear loading paths; these may be, for example, the Bauschinger effect, the permanent softening effect and the latent hardening effect [6,7]. The Bauschinger effect is the phenomenon wherein the yield stress under reverse loading decreases after a certain amount of pre-deformation [8,9]; the permanent softening effect is another phenomenon during reverse loading wherein the reloading flow stress remains lower than the flow stress of monotonic loading [10]. Meanwhile, there is a wide range of strain path changes between monotonic loading and reverse loading. For some materials, the reloading flow stress at an angle with respect to the first loading may overshoot the monotonic stress-strain curve; this phenomenon is referred to the latent hardening effect [11,12].
Practically, the anisotropic hardening behavior significantly affects the spring-back and forming limit [13][14][15], and the reduction of spring-back and the avoidance of necking or fracture are of paramount importance in the forming of an integrated component. In order to guide process formulation and reduce the cost of process development, it is necessary to predict the deformation process and optimize the forming process by simulation. At present, the isotropic hardening assumption is commonly used in general simulations of metal forming, i.e., the subsequent yield surface expands uniformly based on the initial yield surface. However, an isotropic hardening model cannot reasonably capture the anisotropic hardening behavior under nonlinear loading paths [16,17]. Therefore, establishing advanced anisotropic hardening models to accurately predict the plastic deformation behavior of thin-walled metals under nonlinear loading conditions is an important research direction [18].
A variety of anisotropic hardening models have been proposed [19,20], for example, the kinematic hardening model established by introducing a back stress tensor into the isotropic hardening model. The kinematic hardening model was first proposed by Prager to capture the Bauschinger effect, wherein the yield surface translates with the deformation but its shape remains the same [21]. Subsequently, some kinematic hardening models with a nonlinear relationship between the increments of back stress and plastic strain were proposed to more accurately capture the flow stress under reverse loading or cyclical loading conditions [22,23]. In fact, the typical hardening behavior of metals is often a mixture of isotropic hardening and kinematic hardening, i.e., the yield surface changes both in size and position. Hence, the mixed hardening mode was proposed [24,25].
As an alternative to the kinematic hardening model, Barlat et al. [26] proposed the homogeneous anisotropic hardening (HAH) model, which is a distortional hardening model without a back stress tensor. In the HAH model, a fluctuating component is added into the traditional yield function to change the shape of yield surface, and a microstructure deviator related to the loading history is introduced to reflect the state of the microstructure. In order to further describe the latent hardening, more variables were introduced into the model [27,28]. More importantly, the HAH model can be regarded as a theoretical framework for the anisotropic hardening model, in which any homogeneous yield function and hardening law can be used, and the yield surface can be distorted depending on the specific loading path [17]. Some researchers have successfully implemented the HAH model in numerical simulations, and the prediction accuracy has improved clearly [29][30][31].
More recently, in order to capture complex and varied anisotropic hardening behavior and improve prediction accuracy, more parameters were introduced, for example, the HAH-20 model (using, at most, 23 coefficients [17]) and the HEXAH model (using, at most, 15 coefficients [32]). On the other hand, in order to reduce the computational cost in numerical simulation, some efficient algorithms were implemented, including a multi-step return-mapping algorithm for the HAH model [33], a fully implicit numerical algorithm that can solve a complete set of residuals [34], and a fast and robust stress-update algorithm based on the general cutting-plane method [35].
It should be pointed out that the HAH model and these improvements were proposed based on the second-order stress deviator with six independent components, corresponding to the fully 3D (three dimensional) stress state. However, in the plastic forming of thinwalled components, the normal stress is usually much smaller than the in-plane principal stresses, which means that the stress condition can be simplified to the plane stress state with only three components. Therefore, shell elements and material constitutive models for plane stress states are widely used in numerical simulations of sheet metal forming. When the HAH model is selected to analyze the deformation of sheet metals under nonlinear loading paths, it is required to transform the in-plane stress components to a stress deviator, then calculate the strain increments and new stress deviator after a deformation increment using an iterative method, and lastly, transform the new stress deviation tensor into inplane stress components for output. Obviously, the plane stress problem is transformed into a more complex 3D stress problem, which makes the calculation complicated and costly. Therefore, it is of great significance to establish an anisotropic hardening model in plane stress state to analyze the deformation of sheet metals subjected to complex loading.
In this paper, an anisotropic hardening model for plane stress problems such as thin-walled metal forming will be proposed based on the framework of the HAH model. For convenience, the proposed model is named the HAH-2d model in this work. In the meantime, the performance of the HAH-2d model in predicting the evolution of yield loci and the stress-strain curve will be analyzed under typical nonlinear loading paths. Moreover, the new model will be verified by experimental results from two-step uniaxial tension tests.

Fundamentals of the HAH Model
The original HAH model was proposed by Barlat et al. in the second-order stress deviator space [26]. The yield function, as well as the plastic potential, is as follows: where, s andĥ are the stress deviator and microstructure deviator, respectively, and ":" denotes the double dot product. The microstructure deviatorĥ, a normalized tensorial state variable, is proposed to capture the deformation history and reflect the microstructure evolution of the material. It defines an axis, not a direction; namely,ĥ and −ĥ represent the same microstructural state.
Its initial valueĥ 0 is defined as the normalized stress deviatorŝ 0 corresponding to the initial yield, as in the following equation, where the factor 8/3 is used for convenience [26,27].
If the material is reloaded in a different stress state,ĥ will rotate gradually towards the new stress deviator s when cos χ ≥ 0, or towards −s when cos χ < 0. Additionally, cos χ determined by Equation (3) is the cosine of the angle between s andĥ, representing the variation of the loading path [27,28]. Monotonic, reverse and cross-loading sequences are, respectively, represented by cos χ = 1, −1 and 0.
where, s = √ s : s = √ s ij s ij denotes the norm of the second-order tensor s.

Stress Vector and Microstructure Vector
During the forming process of sheet metal, the material is mainly in plane stress state; the three in-plane stress components can be simply expressed as a stress vector σ = (σ 11 , σ 22 , σ 12 ) T , while the other stress components such as σ 33 , σ 13 and σ 23 are equal to 0. In order to describe the microstructure evolution of the sheet metal under plane stress conditions, a normalized microstructure vectorĤ = (Ĥ 11 ,Ĥ 22 ,Ĥ 12 ) T is defined in terms of stress vector in this paper. Its initial valueĤ 0 also corresponds to the normalized stress vectorσ 0 leading to the first increment of plastic deformation.
Considering that the plastic deformation of metals is driven by the stress deviator s, it is still necessary to discuss the variation of loading path based on stress deviator when analyzing plastic deformation processes under plane stress conditions. The second-order deviatoric tensors, s and h, corresponding to the stress vector σ and the microstructure vectorĤ, can be obtained through the following two equations, respectively.
Then, the double dot product s : h can be expressed as: where the matrix D = 1 Therefore, the value of cos χ can be calculated through vectors σ andĤ, as follows: Moreover,Ĥ defines an axis in the vector space of plane stress. When the sheet metal undergoes loading-path changes,Ĥ will remain the same if cos χ = ±1, or rotate gradually towards the new stress vector σ if 0 ≤ cos χ < 1, or rotate towards −σ if −1 ≤ cos χ < 0.
In this work, a possible evolution forĤ is given as follows: in whichσ is the normalized stress vector of σ, k is a constant that controls the rate of rotation, and d − ε is the increment of equivalent strain − ε , which can be calculated according to the plastic work increment dW p : where dε is the plastic strain increment vector, dε 11 and dε 22 are strain increments along directions 1 and 2, and dγ 12 is the shear strain increment.
The evolution ofĤ is shown schematically in Figure 1.Ĥ rotates in the plane determined byĤ i andσ, and finally,Ĥ takes the direction ofσ or −σ.  (11) in which σ is the normalized stress vector of σ , k is a constant that controls the rate of rotation, and dε is the increment of equivalent strain ε , which can be calculated according to the plastic work increment p dW : where dε is the plastic strain increment vector, 11 dε and 22 dε are strain increments along directions 1 and 2, and 12 dγ is the shear strain increment.
The evolution of Ĥ is shown schematically in Figure 1. Ĥ rotates in the plane determined by ˆi H and σ , and finally, Ĥ takes the direction of σ or −σ .

Distortional Yield Function
The distortional yield function of the HAH model in the plane stress space (HAH-2d model) is as follows:

( )
Φ σ is also a homogeneous function of degree 1,

Distortional Yield Function
The distortional yield function of the HAH model in the plane stress space (HAH-2d model) is as follows: The stable component φ can be any homogeneous yield function of degree 1 for plane stress. The fluctuating component φ h is adopted to cover the Bauschinger and permanent softening effects, and the multiplier component 1/F L is used to capture the latent hardening effect. Correspondingly, Φ(σ) is also a homogeneous function of degree 1, i.e., Φ(aσ) = aΦ(σ) holds for any real number a. In addition, σ p is the projection of the stress deviator s in the direction of the microstructure deviator h.
It is important to note that the associated flow rule is applied in this paper, the yield function also serves as the plastic potential. In other words, the plastic strain increment vector dε is defined as where dλ ≥ 0 is the plastic multiplier.

Bauschinger Effect and Permanent Softening Effect
A schematic diagram of yield loci in the σ 11 − σ 22 plane with f 1 > 0, f 2 > 0 and F L = 1 is shown in Figure 2. σ f and σ r are yield stress vectors on the current distortional yield locus along the positive and negative directions ofĤ, respectively, and σ iso is the yield stress vector on the isotropic yield locus. State variables g 1 and g 2 are defined as Equation (17) to represent the relationship between σ f , σ r and σ iso , and the parameters f 1 and f 2 in Equation (14) can be expressed as Equation (18). i.e., σ holds for any real number a. In addition, p σ is the projection of the stress deviator s in the direction of the microstructure deviator h.
It is important to note that the associated flow rule is applied in this paper, the yield function also serves as the plastic potential. In other words, the plastic strain increment vector dε is defined as where 0 dλ ≥ is the plastic multiplier.

Bauschinger Effect and Permanent Softening Effect
A schematic diagram of yield loci in the 11 Figure 2. f σ and r σ are yield stress vectors on the current distortional yield locus along the positive and negative directions of Ĥ , respectively, and iso σ is the yield stress vector on the isotropic yield locus. State variables 1 g and 2 g are defined as Equation (17) to represent the relationship between f σ , r σ and iso σ , and the parameters 1 f and 2 f in Equation (14) can be expressed as Equation (18). For a sheet metal exhibiting the Bauschinger effect, the normalized yield surface in the opposite direction of loading tends to contract during forward loading, and the value of 1 f should increase from 0 towards a finite value. As a result, the yield stress of reverse loading is reduced compared with the final forward-loading stress, which is the Bauschinger effect. During reverse loading, the flat yield surface at the side of reverse loading tends to recover the yield surface determined by the isotropic hardening function ( ) . At the same time, the yield surface opposite to the reverse loading contracts, For a sheet metal exhibiting the Bauschinger effect, the normalized yield surface in the opposite direction of loading tends to contract during forward loading, and the value of f 1 should increase from 0 towards a finite value. As a result, the yield stress of reverse loading is reduced compared with the final forward-loading stress, which is the Bauschinger effect. During reverse loading, the flat yield surface at the side of reverse loading tends to recover the yield surface determined by the isotropic hardening function φ(σ) = − σ. At the same time, the yield surface opposite to the reverse loading contracts, and f 2 starts to increase in the same way as f 1 . If the yield surface at the side of reverse loading cannot recover to the level of isotropic hardening, the permanent softening effect will occur during reverse loading correspondingly.
A possible evolution for g 1 and g 2 is given in Equation (19) which has been successfully applied in existing HAH models [26][27][28]. State variables g 3 and g 4 represent the maximum saturation values of g 2 and g 1 , respectively. Material coefficients k 1 and k 2 control the evolution rate of g 1 and g 2 , and k 3 controls their lower bound; k 4 controls the minimum value of g 3 and g 4 , and k 5 controls their evaluation rate. Meanwhile, the value of k 4 is usually slightly less than 1, according to the permanent softening effect. When the initial values of g 1 and g 2 are both 1 and k 1 = k 2 = 0, the ability to describe the Bauschinger effect will be suppressed, and when g 3 = g 4 = 1 and k 4 = 1 (or k 5 = 0), the ability to describe the permanent softening effect will be suppressed.

Latent Hardening Effect
For some materials, the flow stress during reloading may overshoot the stress during monotonic loading under the condition of the same equivalent strain. The phenomenon of overshooting was explained as the latent hardening effect. Additionally, the degree of overshooting is usually maximum for cross loading, in which the angle between the two corresponding normalized stress deviatorsŝ 1 andŝ 2 is 90 • (ŝ 1 :ŝ 2 = 0). As the degree of plastic deformation increases in the second step,ĥ rotates gradually towardsŝ 2 or −ŝ 2 .
Suppose s v and s ' v are the stress deviators vertical toĥ on the distortional yield surface and the isotropic one, respectively. Additionally, the corresponding plane stress vectors are σ v and σ ' v . Then, a state variable g L can be defined as Equation (20), which represents that the distortional yield surface at cos χ = 0 expands to g L times the yield surface under monotonic loading, and its minimum value is 1.
Because the degree of the latent hardening effect is different under different strain path changes, it is important to determine the yield surface expansion for arbitrary χ. The stress deviator s can be decomposed into s parallel toĥ and s ⊥ perpendicular toĥ, as shown in Figure 3 [28]. Considering that the latent hardening effect is mainly caused by the vertical component s ⊥ , a new deviator s ' is composed as Then, the expansion ratio of the yield surface for arbitrary χ can be represented by ε under any strain path change except cos χ = ±1. Furthermore, the evolution equation of g L is as follows [27]: where L and k L are constants controlling the upper bound and evolution rate of g L , respectively. When the initial value of g L is 1 and L = 1, F L will be identical to 1, and the latent hardening effect is suppressed.  Then, the expansion ratio of the yield surface for arbitrary χ can be represented by Thus, the surface of ( ) ( ) under any strain path change except cos 1 χ = ± . Furthermore, the evolution equation of L g is as follows [27]: where L and L k are constants controlling the upper bound and evolution rate of L g , respectively. When the initial value of L g is 1 and L = 1, L F will be identical to 1, and the latent hardening effect is suppressed.

Coefficient Identification
It is recommended to identify the coefficients in the HAH-2d model in the following order: yield function ( ) φ σ , isotropic hardening law ( ) σ ε , and distortional hardening [27]. The coefficients of ( ) φ σ and ( ) σ ε can be determined independently in a conventional way; for example, the coefficients of ( ) φ σ can be calculated based on the yield stresses and anisotropy coefficients along different directions of the sheet, and the coefficients of ( ) σ ε can be obtained by fitting a monotonic stress-strain curve. Additionally, the nine coefficients related to the distortion, namely, q, k, k1, k2, k3, k4, k5, kL and L, can be identified using an optimization method in which the input experimental data should come from reverse-loading and cross-loading tests.

Coefficient Identification
It is recommended to identify the coefficients in the HAH-2d model in the following order: yield function φ(σ), isotropic hardening law − σ − ε , and distortional hardening [27].
The coefficients of φ(σ) and − σ − ε can be determined independently in a conventional way; for example, the coefficients of φ(σ) can be calculated based on the yield stresses and anisotropy coefficients along different directions of the sheet, and the coefficients of − σ − ε can be obtained by fitting a monotonic stress-strain curve. Additionally, the nine coefficients related to the distortion, namely, q, k, k 1 , k 2 , k 3 , k 4 , k 5 , k L and L, can be identified using an optimization method in which the input experimental data should come from reverse-loading and cross-loading tests.
When the stress-strain data of forward-reverse loading are available, coefficients k 1 , k 2 , k 3 , k 4 and k 5 can be determined independently. Otherwise, the five coefficients have to be evaluated together with k, k L , L and q. In particular, the Bauschinger, permanent softening and latent hardening effects can be expressed simultaneously in the stress-strain curve from a two-step tensile test with −1 < cos χ < 0, for example, re-tension in the direction orthogonal to the fist tensile direction (cos χ = −0.5). Given that the orthogonal two-step tensile test is easy to carry out, all the coefficients associated with the distortion could be identified based on the stress-strain curve.

Hardening Behavior Predicted under Typical Nonlinear Loading Paths
In this section, the stress-strain curves and yield locus evolutions of a generic sheet material under several typical nonlinear loading paths were predicted by the HAH-2d model. The typical nonlinear loading paths in the range of −1 ≤ cos χ < 1 include reverse loading (cos χ = −1), cross loading (cos χ = 0) and two-step tensions conducted at 45 • (cos χ = 0.25) and 90 • (cos χ = −0.5) from the first loading direction.
The initial yield condition of the generic sheet material is assumed to satisfy the Mises isotropic yield criterion, see Equation (24), and the hardening satisfies the Swift hardening law [36], see Equation (25). The material has a Bauschinger effect, permanent softening effect and latent hardening effect simultaneously. The coefficients of the swift hardening law and coefficients associated with the distortion are listed in Table 1.
where K, ε 0 and n are material coefficients.

Reverse Loading
The evolution of the normalized yield locus of the generic material, predicted by the HAH-2d model under a reverse loading, namely, uniaxial compression (UC) to true strain of 0.10, followed by uniaxial tension (UT) in the reverse direction, is shown in Figure 4.
The microstructure vectorĤ maintains (−1, 0, 0) during the whole deformation. The normalized yield loci in σ 11 − σ 22 plane are divided into left and right sides by the line σ 22 /σ 11 = 2.0. This line is chosen because the two stress deviators corresponding to stress vectors (1, 2, 0) T and (1, 0, 0) T are orthogonal. During the compression step, the yield locus on the left side remains unchanged, while the yield locus on the right side contracts inward, which will result in a reduction in re-loading yield stress when tension occurs in the opposite direction. As the tensile strain increases in the second step, the yield locus on the left side contracts inward instead. The yield locus on the right side recovers towards the isotropic hardening yield, but not completely, because of the permanent softening effect.  The predicted stress-strain curve during reverse loading is given in Figure  reloading yield stress after compression pre-deformation is significantly lower th The predicted stress-strain curve during reverse loading is given in Figure 5. The reloading yield stress after compression pre-deformation is significantly lower than the monotonic flow stress under the same deformation level, which is the Bauschinger effect. When k 4 = 0.9, the stress will not reach the monotonic uniaxial tensile curve, which is the permanent softening effect. Additionally, the permanent softening phenomenon disappears when k 4 = 1.0. Therefore, the HAH-2d model can describe the Bauschinger effect and permanent softening effect during reverse loading. The predicted stress-strain curve during reverse loading is given in Figure 5. The reloading yield stress after compression pre-deformation is significantly lower than the monotonic flow stress under the same deformation level, which is the Bauschinger effect. When k4 = 0.9, the stress will not reach the monotonic uniaxial tensile curve, which is the permanent softening effect. Additionally, the permanent softening phenomenon disappears when k4 = 1.0. Therefore, the HAH-2d model can describe the Bauschinger effect and permanent softening effect during reverse loading. The sensitivity of the HAH-2d model to parameters k1, k2, and k3 should be considered when predicting the Bauschinger effect during reverse loading. The state of the yield surface associated with the Bauschinger effect is reflected by state variables g1 and g2, and their evolutions follow the same rule. Therefore, the effects of k1, k2, and k3 on the predicted responses, such as the evaluations of g1 and the stress-strain curves, were analyzed as The sensitivity of the HAH-2d model to parameters k 1 , k 2 , and k 3 should be considered when predicting the Bauschinger effect during reverse loading. The state of the yield surface associated with the Bauschinger effect is reflected by state variables g 1 and g 2 , and their evolutions follow the same rule. Therefore, the effects of k 1 , k 2 , and k 3 on the predicted responses, such as the evaluations of g 1 and the stress-strain curves, were analyzed as shown in Figure 6. It can be found from Figure 6a,b that k 1 affects the evaluation rate of g 1 during reloading. The greater k 1 is, the greater the increase rate of g 1 , and the faster the flow stress reaches the saturation state. Figure 6c shows that k 2 affects the evaluation rate of g 1 during preloading. However, the corresponding stress-strain curves of reloading under the condition of 30 ≤ k 2 ≤ 50, k 1 = 100 and k 3 = 0.5 are almost coincident, as shown in Figure 6d. The reason is that there is little difference in g 1 after a pre-strain of 0.1, and the difference is quickly covered by the evolution during reverse loading. For the same reason, the difference among the curves in Figure 6f under different k 3 is not obvious too. Even so, it should be noted that the value of g 1 during preloading and the reloading yield stress are significantly affected by k 3 ; the smaller the value of k 3 , the smaller the value of g 1 and the reloading yield stress. Figure 7 shows the effects of k 4 and k 5 on the predicted responses of the permanent softening effect during reverse loading, where the state variable g 4 represents the maximum saturation value of the reloading stress. It can be seen that the permanent softening effect is mainly affected by k 4 , and the smaller the value of k 4 , the lower the value of g 4 and the saturation value of reloading stress. In the meantime, when using a value of k 4 less than 1.0, the larger the value of k 5 , the faster the value of g 4 approaches the saturation value, and the lower the reloading stress-strain curve will be after a certain pre-strain. Of course, if the pre-strain is large enough, the influence of k 5 on the reloading stress-strain curve will be significantly reduced or may even disappear.  mum saturation value of the reloading stress. It can be seen that the permanent softening effect is mainly affected by k4, and the smaller the value of k4, the lower the value of g4 and the saturation value of reloading stress. In the meantime, when using a value of k4 less than 1.0, the larger the value of k5, the faster the value of g4 approaches the saturation value, and the lower the reloading stress-strain curve will be after a certain pre-strain. Of course, if the pre-strain is large enough, the influence of k5 on the reloading stress-strain curve will be significantly reduced or may even disappear.

Cross Loading
It is known that the two stress deviators corresponding to ，， are orthogonal to each other, and the corresponding deformation types are plane strain tension and uniaxial tension, respectively. Therefore, the loading sequence of a plane strain tension (PT) followed by a uniaxial tension at an angle of 90° from the first tension direction is a typical cross loading path, and cos 0 χ = at the start of the second loading. The evolution of the normalized yield locus of the generic material, predicted by the HAH-2d model during uniaxial tension after plane strain tension to equivalent strain of 0.10, is shown in Figure 8

Cross Loading
It is known that the two stress deviators corresponding toσ PT = (1/ √ 5, 2/ √ 5, 0) T andσ UT = (1, 0, 0) T are orthogonal to each other, and the corresponding deformation types are plane strain tension and uniaxial tension, respectively. Therefore, the loading sequence of a plane strain tension (PT) followed by a uniaxial tension at an angle of 90 • from the first tension direction is a typical cross loading path, and cos χ = 0 at the start of the second loading. The evolution of the normalized yield locus of the generic material, predicted by the HAH-2d model during uniaxial tension after plane strain tension to equivalent strain of 0.10, is shown in Figure 8. In the plane strain tension process, there isĤ =σ PT , the yield locus above the line σ 22 = 0 remains unchanged, and the yield locus below the line σ 22 = 0 contracts inward because of the Bauschinger effect. In the second loading step,Ĥ rotates gradually towardsσ UT with the increase in tensile strain, and the yield locus nearσ UT expands outward first and then returns to the isotropic hardening yield locus. locus above the line 22 0 σ = remains unchanged, and the yield loc 22 0 σ = contracts inward because of the Bauschinger effect. In the se Ĥ rotates gradually towards UT σ with the increase in tensile strain, near UT σ expands outward first and then returns to the isotropic hard The uniaxial tensile stress-strain curves in the second loading step after different degrees of plane strain pre-tension are shown in Figure 9. Obvious stress overshooting is observed, and the higher the pre-deformation the greater the overshooting. When the tensile strain reaches about 0.06 in the second step, the tensile stress returns to the level of single uniaxial tension. It indicates that the HAH-2d model can capture the latent hardening in cross loading. The uniaxial tensile stress-strain curves in the second loading step after diffe grees of plane strain pre-tension are shown in Figure 9. Obvious stress oversho observed, and the higher the pre-deformation the greater the overshooting. When sile strain reaches about 0.06 in the second step, the tensile stress returns to the single uniaxial tension. It indicates that the HAH-2d model can capture the latent h ing in cross loading. The latent hardening effect predicted by the HAH-2d model during cross loa affected mainly by parameters k, L, and kL. Therefore, the sensitivity of the model Figure 9. Stress-strain curves predicted by the HAH-2d model in plane strain tension-uniaxial tension with different pre-strain levels.
The latent hardening effect predicted by the HAH-2d model during cross loading is affected mainly by parameters k, L, and k L . Therefore, the sensitivity of the model to these three parameters was analyzed. Figure 10 shows the predicted responses of latent hardening effect under different values of k during a typical cross loading, PT0.1−UT. With the increase in the reloading strain, the value of χ gradually decreases from 90 • to 0 • , and the value of g L rapidly increases to a certain value and then gradually decreases to 1.0. As a result, the phenomenon of stress overshooting is observed. Meanwhile, the larger the value of k, the faster the value of χ and g 4 decrease, and the faster the overshoot stress recovers.  Figure 11 shows the effects of L and k L on the predicted responses of latent hardening effect during cross loading. It can be seen from Figure 11a,b that the larger the value L, the larger g L can be obtained and the more significant the stress overshooting. On the other hand, the latent hardening effect is also influenced by the value of k L , as shown in Figure 11c,d. Additionally, the larger the value of k L , the faster the stress value reaches its peak, and the larger the maximum stress value is. However, the sensitivity of the model to k L is much lower than that to L. In addition, neither L nor k L affects the strain span of the phenomenon of stress overshooting.

Two-Step Uniaxial Tension
Stress-strain curves in two-step uniaxial tensions with tensile axes at 45 • and 90 • from each other were predicted by the HAH-2d model and are shown in Figures 12 and 13, respectively. The curves with 45 • express only latent hardening, while the curves with 90 • express the Bauschinger effect, permanent softening and latent hardening at the same time. Therefore, it is better to optimize the coefficients of the HAH-2d model by using the stress-strain data of a two-step uniaxial tension with 90 • , or the orthogonal two-step tension mentioned in Section 3.3.

Determination of Material Coefficients
In this paper, the HAH-2d model will be applied to describe the hardening behavior of DP780 and EDDQ steel sheets with 1.2 mm thickness that were tested by Ha et al. using two-step uniaxial tension tests [6]. As shown in Figure 14, the sub-specimens were cut along different directions at the center of the pre-deformed big-specimen [6]. Additionally, the experimental results from tests consisted first of tension in the rolling direction (RD) and second of tension at 90° (TD, the transverse direction), 60° and 45° from RD were used in this work. Coefficients of the Yld2000-2d anisotropic yield function (see Appendix A) and the Swift hardening law (see Equation (25)) for the two materials are listed in Table  2, where m = 6 was adopted as recommended for BCC metals [6].

Determination of Material Coefficients
In this paper, the HAH-2d model will be applied to describe the hardening behavior of DP780 and EDDQ steel sheets with 1.2 mm thickness that were tested by Ha et al. using two-step uniaxial tension tests [6]. As shown in Figure 14, the sub-specimens were cut along different directions at the center of the pre-deformed big-specimen [6]. Additionally, the experimental results from tests consisted first of tension in the rolling direction (RD) and second of tension at 90 • (TD, the transverse direction), 60 • and 45 • from RD were used in this work. Coefficients of the Yld2000-2d anisotropic yield function (see Appendix A) and the Swift hardening law (see Equation (25)) for the two materials are listed in Table 2, where m = 6 was adopted as recommended for BCC metals [6].  To quantitatively evaluate the error between the predicted and experimental data, the root mean square error (RMSE) and the average absolute relative error (AARE) were used in this paper. The definition of RMSE and AARE can be expressed by Equations (26) and (27).  To quantitatively evaluate the error between the predicted and experimental data, the root mean square error (RMSE) and the average absolute relative error (AARE) were used in this paper. The definition of RMSE and AARE can be expressed by Equations (26) and (27).
where σ P is the predicted stress, σ E is the experimental stress, and N is the number of data. The coefficients associated with the distortion of the HAH-2d model for each material were identified with an optimization method using the stress-strain data from an orthogonal two-step tensile test with 10% pre-strain in RD, as shown in Figure 15. In the optimization method, the objective is to minimize the RMSE. Considering that the DP780 steel expresses no permanent softening and latent hardening, k 4 = 1.0, k 5 = 0, k L = 0 and L = 1.0 were set directly. Similarly, as the EDDQ steel expresses no permanent softening, k 4 = 1.0 and k 5 = 0 were set. In the meantime, the value of q was set to 2.0 for convenience. The obtained coefficients of the HAH-2d models for DP780 and EDDQ are listed in Table 3. Additionally, the predicted stress-strain curves are plotted in Figure 15a,b, both of which agree well with the experimental data. Additionally, the values of RMSE for DP780 and EDDQ are 9.81 and 1.85, respectively. Incidentally, the reason for using the experimental data from reference [27] in Figure 15 is that the experimental data used in reference [27] are also from reference [6], but more points are given for this testing condition.

Comparison of Prediction and Experimental Results
In order to validate the HAH-2d model, the hardening curves in two-step uniaxial tensions of DP780 and EDDQ sheets, which consisted of 4% or 10% pre-strain in RD and re-tension at 90 • (TD), 60 • or 45 • from RD, were predicted and compared with the experimental results. Incidentally, the two-step tension with tensile axes at 60 • from each other can be called "pseudo cross loading", because it is close to the ideal cross-loading condition, about 55 • between the first and second tensile directions.
The predicted stress-strain curves of the HAH-2d model for DP780 under different two-step uniaxial tensions are compared with the experimental curves in Figure 16, and the corresponding prediction errors are listed in Table 4. It can be seen from Figure 16a that the HAH-2d model captures the transient hardening behavior of DP780 fairly well for re-tension in TD after a 4% pre-strain in RD, and the AARE value is only 1.28%. However, the prediction at 60 • shown in Figure 16b slightly overestimates the re-loading yield stress and leads to a lower hardening rate. As a result, the values of RMSE and AARE increase to 37.85 and 3.44% under the condition of 10% pre-strain. On the other hand, as shown in Figure 16c, the phenomenon of cross-loading softening was observed in re-loading stressstrain curves at 45 • , but the predicted re-loading yield stresses at 45 • reach the monotonic curve immediately. Additionally, the high RMSE values of 24.49 and 47.13 are mainly due to the significant difference between the predicted stress and the experimental stress at the beginning of the second tension.   Figure 17a-c show the predicted reloading flow stress curves for EDDQ in different two-step uniaxial tensions with pre-tension in RD and re-tension at 90 • (TD), 60 • and 45 • from RD. The predicted curves capture all the experimental results well, even though there is a little overestimation of the flow stress in the re-tension at 60 • from RD. The corresponding prediction errors are listed in Table 5. The maximum values of RMSE and AARE, 9.86 and 2.00%, both occur in the condition of re-tension at 60 • after a 10% pre-tension. This indicates that the HAH-2d model can capture anisotropic hardening responses such as flow stress overshooting and strain-hardening stagnation.

Conclusions
In this work, an HAH model in plane stress state was proposed, and the evolution of the yield loci and the stress-strain curve of a generic material under some typical nonlinear loading paths was analyzed using this model. Furthermore, the model was validated by experimental results from two-step uniaxial tension tests. The conclusions are as follows: (1) The HAH-2d model was proposed in plane stress state. Both the loading stress vector and microstructure vector contain only three in-plane components. Compared with the HAH model in fully 3D space, the transformation between the plane stress vector and the second-order stress deviator and the complicated calculation in the stress deviator space can be avoided in the HAH-2d model. The computational cost in finite element applications is therefore reduced.
(2) The Bauschinger effect and permanent softening effect in reverse loading and the latent hardening effect in cross loading were investigated in this paper. The evolution of the yield surface and the stress-strain curve of the sheet metal with these effects can be predicted by the HAH-2d model.
The principal values of X are