Fully Implicit Stress Update Algorithm for Distortion-Based Anisotropic Hardening with Cross-Loading Effect: Comparative Algorithmic Study and Application to Large-Size Forming Problem

Abstract

A fully implicit stress integration algorithm is developed for the distortional hardening model, namely the e-HAH model, capable of simulating cross-hardening/softening under orthogonal loading path changes. The implicit algorithm solves a complete set of residuals as nonlinear functions of stress, a microstructure deviator, and plastic state variables of the constitutive model, and provides a consistent tangent modulus. The number of residuals is set to be 20 or 14 for the continuum or shell elements, respectively. Comprehensive comparison programs are presented regarding the predictive accuracy and stability with different numerical algorithms, strain increments, material properties, and loading conditions. The flow stress and r-value evolutions under reverse/cross-loading conditions prove that the algorithm is robust and accurate, even with large strain increments. By contrast, the cutting-plane method and partially implicit Euler backward method, which are characterized by a reduced number of residuals, result in unstable responses under abrupt loading path changes. Finally, the algorithm is implemented into the finite element modeling of large-size, S-rail forming and the springback for two automotive steel sheets, which is often solved by a hybrid dynamic explicit-implicit scheme. The fully implicit algorithm performs well for the whole simulation with the solely static implicit scheme.

1. Introduction

Automotive, non-ferrous alloys and advanced high-strength steels have been extensively investigated by many researchers owing to their lightweight and excellent mechanical properties. Lightweight alloys, such as aluminum and magnesium alloys, exhibit superior specific strengths to conventional steels. By contrast, the increased strength, which can exceed 1 GPa, enables the advanced high-strength steel (AHSS) sheets to have a thinner gauge, thus, rendering them more competitive compared to low-density, lightweight alloys. The excellent strength and enhanced ductility of the AHSS can be achieved owing to their well-controlled microstructure and deformation mechanisms. For example, dual-phase steels or transformation induced steels introduce a hard martensite phase embedded in the softer ferrite phase or transforming martensite from meta-stable austenite. The strength of lightweight alloys can be achieved through the dispersion of optimally designed second phases, such as precipitations.

However, even in the presence of these improved plasticity properties, such as tensile strength and elongation at fracture, the intricate microstructure leads to a further complexity in the mechanical properties, especially when the deformation paths become convoluted. The complex mechanical responses include the Bauschinger effect and transient flow behavior under reversed loading path [1,2,3,4,5]. These anisotropic hardening behaviors of sheet metals have been reported to be critical factors for the springback of automotive sheet parts [6,7,8,9]. More recently, the anisotropic hardening behavior of advanced sheet metals has been further investigated under more complex loading paths than simple reverse loading. For instance, the deformation of ferritic-phased steels shows significant stress overshooting (or larger stress than the monotonic flow stress) when the loading path changes to being orthogonal to the previous loading path [10,11,12,13,14,15]. Interestingly, dual-phase steels with mixed ferrite and martensite phases exhibit a clear softening under the same loading path changes [15,16,17,18]. These plastic behaviors under orthogonal loading path changes are called cross-hardening and cross-softening, respectively. Experimental observations on these complex mechanical behaviors, which cannot be explained through the simple isotropic plasticity model, require the implementation of new anisotropic models in the field of metal forming and plasticity. Indeed, this complex behavior can influence the formability and springback of advanced sheet metals [9].

In the literature, a significant number of works have proposed models for anisotropic hardening behavior. Kinematic hardening is a representative concept, which explains the Bauschinger effect and transient flow hardening at load reversal by introducing yield surface translation during plastic deformation. The kinematic hardening models were pioneered by Prager [19] and Ziegler [20], and further extended by adding nonlinear terms [21] or by coupling with the isotropic hardening model. The series of isotropic and kinematic hardening models was well implemented into finite element (FE) simulations for sheet metals, especially when metals exhibit the Bauschinger effect, transient behavior, and permanent softening under reversed loading paths [22,23,24,25].

The kinematic hardening—based anisotropic hardening models were further extended by combing the distortional hardening concept to reproduce mechanical responses under rather complex loading conditions, beyond the simple loading-reverse loading path. Ortiz and Popov [26] introduced the distortion of the yield surface by controlling the size of the effective stress. Feigenbaum and Dafalias [27,28] employed the fourth-rank anisotropic tensor as a function of plastic deformation, but the fundamental basis they used remained that of isotropic-kinematic hardening. Teodosiu and Hu [29] introduced new effective values into the yield condition related to the structure and interactions of dislocations as a major plasticity mechanism. François [30] expressed the egg-shaped yield function by decomposing deviatoric stress into its collinear and orthogonal parts with respect to the back stress of kinematic hardening. Badreddine et al. [31] developed a non-associated elastoplastic anisotropic hardening model, which is coupled with a damage model based on François’s approach [30]. Qin et al. [32] suggested a model to represent the Bauschinger effect with the kinematic hardening component, while other strain-path change effects could be expressed through the distortion of the yield surface.

Besides the kinematic hardening or combined kinematic-distortional hardening, several anisotropic hardening models are solely based on the distortional hardening approach. Some of the distortional hardening models were developed to express the yield surface evolution by using the anisotropic coefficients as a function of the plastic work or equivalent plastic strain [33,34,35,36]. However, these models do not take into account the loading path change effect. Barlat and one the of co-authors of the present study [37,38,39] proposed a series of anisotropic hardening models without yield surface translation, which they named the homogeneous yield function based anisotropic hardening (HAH) models. The main concept of the original HAH model consists of using the distortion of the yield surface along a designated loading path, which is called a microstructure deviator, and the distortional amounts are controlled via adequately defined plastic state variables. Later, the model was extended for reproducing latent hardening, work hardening stagnation [38], as well as cross-hardening and softening under more general loading conditions [39]. This improved model was named the extended HAH model (or e-HAH, in this study). The performance of the distortional hardening models was validated through many applications, including the springback in U-draw bending and industrial S-rail forming [40,41,42,43,44,45,46]. The model is also applicable to subjects for which non-linear strain path effects are important, such as the fracture behaviors of metal after pre-deformation [47,48,49].

As the constitutive models of plasticity have advanced, especially the hardening laws, their implementation into FE simulations have become increasingly challenging. This is mainly attributed to the advanced constitutive models having more state variables, which are non-linearly cross-related. These challenges have led to more robust numerical formulations and implementations of the constitutive models, which eventually determine the accuracy and robustness of FE simulations. Numerous numerical algorithms were proposed to take into account the stress integration or stress update, using the elastic-plasticity constitutive models. The details of the general theoretical studies on these stress-integration algorithms for nonlinear plasticity are well documented in a book by Simo [50].

The basic principle of the stress update algorithm for the classical rate-independent elastoplastic model consists of locating the stress state on the yield surface described in the six-dimensional stress space, which is consistent with the material hardening. The hardening of the material is often represented by a uniaxial stress-plastic-strain curve as a reference stress state. Most of the stress update algorithms were developed in the elastic-predictor and plastic-corrector schemes. For example, the closest point projection method (CPPM) [51,52,53] and the cutting-plane method (CPM) [54] are employed as popular stress integration algorithms in FE models. The CPPM is often based on the Euler backward method (EBM), and thus it is an implicit method requiring the first and second derivatives of the yield surface to satisfy both the consistency and flow rule (or normality rule). In contrast, the CPM only satisfies the consistency condition without requiring the second derivative of the yield function; thus, this approach is also referred to as a semi-explicit method.

Regarding the stress integration algorithms on the distortional hardening models, similar approaches, based on either the CPPM or CPM schemes, have been reported. Lee et al. [55] implemented the first version of the HAH (denoted as original HAH hereafter) model into commercial ABAQUS software, using both the CPPM and CPM approaches. This implementation was also extended to the e-HAH model; however, in this case, only the CPM scheme was applied in combination with a sub-stepping numerical method [56]. Recently, Choi and Yoon [57] also reported the implementation of the HAH model, using the Euler backward scheme, but they calculated the derivatives of the yield function using numerical finite differences. More recently, Yoon et al. [58] coupled the CPPM scheme with a line-search algorithm for an updated version of the HAH model [59] to improve numerical efficiency. However, the existing numerical integration schemes do not fully exploit the implicit CPPM scheme. In other words, these approaches only account for the partial number of residuals in order to attain the solution process of linearized equations under the Newton-Raphson method. In this sense, these methods can be regarded as a partial or semi-implicit integration algorithm based on the Euler backward scheme.

In this study, a fully implicit stress integration algorithm based on the Euler backward scheme is proposed for the e-HAH model. To the best of our knowledge, the present work is the first to report this novel scheme. The present algorithmic study differs from previous studies [55,56,57,58] in that the complete set of residuals is established from the stress tensor, microstructure deviator, and other variables controlling the distortion of the yield function. Furthermore, the calculations of the first and second derivatives of the distorting yield function are both analytically expressed. Note that the CPPM or Euler backward scheme presented in previous studies adopted only the consistency condition and flow rule for the residuals. This may lead to non-physical evolutions of the distortional hardening state variables when there is an abrupt change in the loading path. For the comparative study on the effect of the employed stress integration algorithms, four different algorithms were implemented into the FE simulation. These algorithms consist of the CPM and EBM schemes with different numbers of residuals. Again, the EBM scheme is further classified based on either analytical derivatives or finite-difference—based numerical derivatives of the yield function. Comprehensive validations are presented in this study. Firstly, one-element analyses with two different material sets are presented for the investigated stress integration algorithms. The flow stress and anisotropy in deformation were evaluated under two distinctive loading conditions, i.e., reverse and cross-loading paths. Finally, the industrial-size S-rail forming and springback simulations were performed to evaluate the numerical accuracy and computational efficiency of the proposed algorithm. Note that even this large-scale industrial problem was solved using the purely static implicit FE solver, whereas most previous studies employed a dynamic explicit method for the forming process and a static implicit method for the springback, due to the divergence problem resulted from complex contact changes during forming.

The manuscript is organized as follows. In Section 2, a summary of the e-HAH model is presented. In Section 3, the numerical integration algorithms investigated in this study are introduced, alongside a summary of previous studies. Then, in Section 4, the accuracy and stability of the individual algorithms are comparatively evaluated via one-element analyses under two distinctive loading path changes. Finally, a large-scale industrial automotive forming and springback simulation is provided as a benchmark problem to assess the validity of the proposed, fully implicit numerical algorithm.

2. The e-HAH Model

The enhanced homogeneous yield function based anisotropic hardening (e-HAH) model was proposed to simulate the path-dependent anisotropic hardening responses of metallic materials. This model is based on the distortional hardening concept, without depending on kinematic hardening, and is capable of reproducing complex flow behaviors, such as cross-hardening or cross-softening [37,38,39].

The distortion of the initially isotropic yield function is controlled by a second-order tensor and eight scalar variables. The second-order tensor, named the microstructure deviator, ${\widehat{h}}^{\mathbf{s}}$ monitors the loading path changes and determines the direction of the distortion. The scalar variables control the magnitude of the distortions in terms of the loading path changes. In this section, the main equations of the constitutive model, along with the related evolution of the state variables, are briefly summarized. The complete derivations of the constitutive model can be found in [37,38,39].

The yield surface ($\mathsf{\Phi }$) of the e-HAH model is described as the following:

$$\mathsf{\Phi }={\left\{{\left(\sqrt{\mathsf{\psi }{(\mathbf{s})}^2+\mathsf{\psi }{({\mathbf{s}}_{\mathrm{P}})}^2}\right)}^{\mathrm{q}}+{\mathrm{f}}_1^{\mathrm{q}}{\left|{\widehat{h}}^{\mathbf{s}}:\mathbf{s}-\left|{\widehat{h}}^{\mathbf{s}}:s\right|\right|}^{\mathrm{q}}+{\mathrm{f}}_2^{\mathrm{q}}{\left|{\widehat{h}}^{\mathbf{s}}:\mathbf{s}+\left|{\widehat{h}}^{\mathbf{s}}:s\right|\right|}^{\mathrm{q}}\right\}}^{1/\mathrm{q}}=\overline{\mathsf{\sigma }}$$

(1)

where $\overline{\mathsf{\sigma }}$ is the equivalent stress, often identified from the uniaxial stress-plastic-strain curve as a reference state, and $\mathrm{q}$ is a constant. The two variables ${\mathrm{f}}_{\mathrm{i}}={\left({\mathrm{g}}_{\mathrm{i}}^{-\mathrm{q}}-1\right)}^{\mathrm{q}}$ with $\mathrm{i}=1\hspace{0.17em},\hspace{0.17em}2$ control the distortion of the yield surface along the direction of ${\widehat{h}}^{\mathbf{s}}$. At the start of the simulation, ${\widehat{h}}^{\mathbf{s}}$ is initialized as ${\widehat{h}}_{ini}^{\mathbf{s}}=\frac{{\mathbf{s}}_{ini}}{\sqrt{\mathrm{H}{\mathbf{s}}_{ini}:{\mathbf{s}}_{ini}}}$. Here, ${\mathbf{s}}_{ini}$ is the deviatoric stress corresponding to the initiation of the plastic deformation, and $\mathrm{H}$ is a constant, which is set to 8/3 in this work.

In Equation (1), $\sqrt{\mathsf{\psi }{(\mathbf{s})}^2+\mathsf{\psi }{({\mathbf{s}}_{\mathrm{P}})}^2}$ represents the distorted yield function, which is used to reproduce the cross-hardening or softening; in the original HAH model, it is replaced by any first-order homogeneous yield function $\mathsf{\varphi }$ [37]. $\mathsf{\psi }(\mathbf{s})$ and $\mathsf{\psi }({\mathbf{s}}_{\mathrm{P}})$ are defined in Equations (2) and (3), respectively, using the collinear (${\mathbf{s}}_{\mathrm{c}}$) and orthogonal (${\mathbf{s}}_{\mathrm{o}}$) components of the deviatoric stress tensor ($\mathbf{s}$) with respect to ${\widehat{h}}^{\mathbf{s}}$.

$$\mathsf{\psi }(\mathbf{s})=\mathsf{\varphi }({\mathbf{s}}^{″})=\mathsf{\varphi }\left({\mathbf{s}}_{\mathrm{c}}+\frac{1}{{\mathrm{g}}_{\mathrm{L}}}{\mathbf{s}}_{\mathrm{o}}\right),$$

(2)

$$\mathsf{\psi }({\mathbf{s}}_{\mathrm{p}})=\mathsf{\varphi }({{\mathbf{s}}_{\mathrm{p}}}^{″})=\mathsf{\varphi }\left(\frac{4\left(1-{\mathrm{g}}_{\mathbf{s}}\right)}{{\mathrm{g}}_{\mathrm{L}}}{\mathbf{s}}_{\mathrm{o}}\right),$$

(3)

where ${\mathbf{s}}^{″}={\mathbf{s}}_{\mathrm{c}}+\frac{1}{{\mathrm{g}}_{\mathrm{L}}}{\mathbf{s}}_{\mathrm{o}}$ and ${\mathbf{s}}_{\mathrm{p}}=4\left(1-{\mathrm{g}}_{\mathbf{s}}\right){\mathbf{s}}_{\mathrm{o}}$ are the two tensors, and ${\mathrm{g}}_{\mathrm{L}}$ and ${\mathrm{g}}_{\mathrm{S}}$ are the two corresponding state variables. The tensors ${\mathbf{s}}^{″}$ and ${\mathbf{s}}_{\mathrm{p}}$ are used for controlling the cross-hardening and softening, respectively. The two components are expressed as follows:

$${\mathbf{s}}_{\mathrm{c}}=\mathrm{H}\left({\widehat{h}}^{\mathbf{s}}:\mathbf{s}\right){\widehat{h}}^{\mathbf{s}},$$

(4)

$${\mathbf{s}}_{\mathrm{o}}=\mathbf{s}-{\mathbf{s}}_{\mathrm{c}}=\mathbf{s}-\mathrm{H}\left({\widehat{h}}^{\mathbf{s}}:\mathbf{s}\right){\widehat{h}}^{\mathbf{s}},$$

(5)

To measure the loading path change, the indicator $\cos \mathsf{\chi }$ is introduced according to the following:

$$\cos \mathsf{\chi }=\frac{{\widehat{h}}^{\mathbf{s}}:\widehat{\mathbf{s}}}{\sqrt{{\widehat{h}}^{\mathbf{s}}:{\widehat{h}}^{\mathbf{s}}}\sqrt{\widehat{\mathbf{s}}:\widehat{\mathbf{s}}}}=\mathrm{H}\cdot \left({\widehat{h}}^{\mathbf{s}}:\widehat{\mathbf{s}}\right),$$

(6)

where the normalized deviatoric stress tensor is defined as $\widehat{\mathbf{s}}=\frac{\mathbf{s}}{\sqrt{\mathrm{H}\mathbf{s}:\mathbf{s}}}$. The value of $\cos \mathsf{\chi }$ indicates the current loading state against the previous loading path. For example, values of 1, 0, and −1 denote monotonic loading, cross-loading, and reverse loading, respectively.

The evolution laws for the state variables are summarized in three groups. The first group (Equations (7)–(10)) is related to the distortion along ${\widehat{h}}^{\mathbf{s}}$:

$$\frac{\mathrm{d}{\mathrm{g}}_1}{\mathrm{d}\overline{\mathsf{\epsilon }}}={\mathrm{u}}_0\cdot {\mathrm{k}}_2\left({\mathrm{k}}_3\frac{\overline{\mathsf{\sigma }}(0)}{\overline{\mathsf{\sigma }}(\overline{\mathsf{\epsilon }})}-{\mathrm{g}}_1\right)+\left(1-{\mathrm{u}}_0\right)\cdot {\mathrm{k}}_1\frac{{\mathrm{g}}_4-{\mathrm{g}}_1}{{\mathrm{g}}_1},$$

(7)

$$\frac{\mathrm{d}{\mathrm{g}}_2}{\mathrm{d}\overline{\mathsf{\epsilon }}}={\mathrm{u}}_0\cdot {\mathrm{k}}_1\frac{{\mathrm{g}}_3-{\mathrm{g}}_2}{{\mathrm{g}}_2}+\left(1-{\mathrm{u}}_0\right)\cdot {\mathrm{k}}_2\left({\mathrm{k}}_3\frac{\overline{\mathsf{\sigma }}(0)}{\overline{\mathsf{\sigma }}(\overline{\mathsf{\epsilon }})}-{\mathrm{g}}_2\right),$$

(8)

$$\frac{\mathrm{d}{\mathrm{g}}_3}{\mathrm{d}\overline{\mathsf{\epsilon }}}=\left(1-{\mathrm{u}}_0\right)\cdot {\mathrm{k}}_5\left({\mathrm{k}}_4-{\mathrm{g}}_3\right),$$

(9)

$$\frac{\mathrm{d}{\mathrm{g}}_4}{\mathrm{d}\overline{\mathsf{\epsilon }}}={\mathrm{u}}_0\cdot {\mathrm{k}}_5\left({\mathrm{k}}_4-{\mathrm{g}}_4\right),$$

(10)

Here, ${\mathrm{u}}_0$ is a unit step function with respect to $\cos \mathsf{\chi }$. In other words, ${\mathrm{u}}_0=1$ when $\cos \mathsf{\chi }\ge 0$ and ${\mathrm{u}}_0=0$ when $\cos \mathsf{\chi }<0$. Depending on the sign of $\cos \mathsf{\chi }$, the different anisotropic hardening behaviors are expressed using ${\mathrm{g}}_1$ and ${\mathrm{g}}_2$. For example, when $\cos \mathsf{\chi }\ge 0$, ${\mathrm{g}}_1$ and ${\mathrm{g}}_2$ determine the magnitude of the Bauschinger effect and transient behavior, respectively, whereas, ${\mathrm{g}}_3$ and ${\mathrm{g}}_4$ determine the saturating values of ${\mathrm{g}}_1$ and ${\mathrm{g}}_2$, respectively. In the simulation, these four state variables are initialized to take the value of 1. In Equations (7)–(10), ${\mathrm{k}}_{1~5}$ are material constants related to the evolution rates and amounts of ${\mathrm{g}}_{1~4}$. Specifically, ${\mathrm{k}}_1$ controls the recovery rate of the transient behavior, ${\mathrm{k}}_2$ and ${\mathrm{k}}_3$ control the evolution rate and amount of the Bauschinger effect, respectively, and ${\mathrm{k}}_4$ and ${\mathrm{k}}_5$ determine the magnitude and rate of permanent softening, respectively. The model parameters, ${\mathrm{k}}_{1~5}$, can be identified using the load reversal experiments, i.e., tension-compression test and cyclic shear test.

The second group is related to the distortion orthogonal to the ${\widehat{h}}^{\mathbf{s}}$ direction. Two additional state variables, namely ${\mathrm{g}}_{\mathrm{L}}$ and ${\mathrm{g}}_{\mathrm{S}}$, are introduced:

$$\frac{\mathrm{d}{\mathrm{g}}_{\mathrm{L}}}{\mathrm{d}\overline{\mathsf{\epsilon }}}={\mathrm{k}}_{\mathrm{L}}\left\{\frac{\overline{\mathsf{\sigma }}(\overline{\mathsf{\epsilon }})-\overline{\mathsf{\sigma }}(0)}{\overline{\mathsf{\sigma }}(\overline{\mathsf{\epsilon }})}\left(\sqrt{\mathrm{L}\left(1-\mathrm{c}\mathrm{o}{\mathbf{s}}^2\mathsf{\chi }\right)+\mathrm{c}\mathrm{o}{\mathbf{s}}^2\mathsf{\chi }}-1\right)\hspace{0.17em}+\hspace{0.17em}1\hspace{0.17em}-{\mathrm{g}}_{\mathrm{L}}\right\},$$

(11)

$$\frac{\mathrm{d}{\mathrm{g}}_{\mathrm{S}}}{\mathrm{d}\overline{\mathsf{\epsilon }}}={\mathrm{k}}_{\mathrm{S}}\left\{1+\left(\mathrm{S}\hspace{0.17em}-\hspace{0.17em}1\right)\mathrm{c}\mathrm{o}{\mathbf{s}}^2\mathsf{\chi }\hspace{0.17em}-{\mathrm{g}}_{\mathrm{S}}\right\},$$

(12)

where ${\mathrm{k}}_{\mathrm{L}}$ and ${\mathrm{k}}_{\mathrm{S}}$ are the constants controlling the evolution rates of ${\mathrm{g}}_{\mathrm{L}}$ and ${\mathrm{g}}_{\mathrm{S}}$, respectively, whereas the parameters $\mathrm{L}$ and $\mathrm{S}$ are introduced to determine the magnitudes of cross-hardening and softening, respectively. The initial values for the two state variables are ${\mathrm{g}}_{\mathrm{L}}={\mathrm{g}}_{\mathrm{S}}=1$. From Equations (11) and (12), it can be seen that no latent hardening effect is present without loading path change or for $\cos \mathsf{\chi }=1$. In contrast, the yield function contracts along the cross-loading direction, even without a loading path change, and the evolution rate is maximized at the cross-loading state or for $\cos \mathsf{\chi }=0$.

The last group is related to the rotation of ${\widehat{h}}^{\mathbf{s}}$:

$$\begin{array}{c}\frac{\mathrm{d}{\widehat{h}}^{\mathbf{s}}}{\mathrm{d}\overline{\mathsf{\epsilon }}}={\mathrm{u}}_0\cdot \mathrm{k}\cdot \left\{{\left|\cos \mathsf{\chi }\right|}^{1/\mathrm{z}}+{\mathrm{g}}_{\mathrm{R}}\right\}\left\{\widehat{\mathbf{s}}-\mathrm{H}\left({\widehat{h}}^{\mathbf{s}}:\widehat{\mathbf{s}}\right){\widehat{h}}^{\mathbf{s}}\right\}\\ +\left(1-{\mathrm{u}}_0\right)\cdot \mathrm{k}\cdot \left\{{\left|\cos \mathsf{\chi }\right|}^{1/\mathrm{z}}+{\mathrm{g}}_{\mathrm{R}}\right\}\left\{-\widehat{\mathbf{s}}+\mathrm{H}\left({\widehat{h}}^{\mathbf{s}}:\widehat{\mathbf{s}}\right){\widehat{h}}^{\mathbf{s}}\right\}\end{array},$$

(13)

$$\frac{\mathrm{d}{\mathrm{g}}_{\mathrm{R}}}{\mathrm{d}\overline{\mathsf{\epsilon }}}={\mathrm{k}}_{\mathrm{R}}\left\{{\mathrm{k}}_{\mathrm{R}}^{\prime }\left(1-{\cos }^2\mathsf{\chi }\right)-{\mathrm{g}}_{\mathrm{R}}\right\}.$$

(14)

Depending on the sign of $\cos \mathsf{\chi }$, the rotation direction of ${\widehat{h}}^{\mathbf{s}}$ is different, for example, when $\cos \mathsf{\chi }\ge 0$, ${\widehat{h}}^{\mathbf{s}}$ rotates toward the direction of $\widehat{\mathbf{s}}$. The rate of ${\widehat{h}}^{\mathbf{s}}$ is controlled by the parameter $\mathrm{k}$. A state variable ${\mathrm{g}}_{\mathrm{R}}$ is introduced in the e-HAH model to prevent the singular evolution rate while maximizing the evolution rate at the cross-loading condition [39]. The suggested values for $\mathrm{z}$, ${\mathrm{k}}_{\mathrm{R}}$ and ${\mathrm{k}}_{\mathrm{R}}^{\prime }$ are 5, 15, and 0.2, respectively.

The model parameters, ${\mathrm{k}}_{\mathrm{L}}$, $\mathrm{L}$, ${\mathrm{k}}_{\mathrm{S}}$, $\mathrm{S}$ and $\mathrm{k}$ can be calibrated using the experiments with loading path changes other than load reversal experiments, such as multi-step tension test [12,13,15].

3. Stress Integration Algorithms of the e-HAH Model for Elastic-Plastic FE Simulations

3.1. Motivation and General Statement

The e-HAH model features a larger number of plastic state variables than the conventional isotropic hardening model. In the von Mises isotropic model, the only state variable used is the equivalent plastic strain. This increased number of state variables is inevitable for the efficient modeling of path-dependent anisotropic hardening; however, the overall numerical procedure becomes more complex. The stress and microstructure deviator tensors, as well as the eight plastic state variables ($\overline{\mathsf{\epsilon }}$ and ${\mathrm{g}}_{\mathrm{i}}\hspace{0.17em}\left(\mathrm{i}=1~4,\hspace{0.17em}\mathrm{L},\hspace{0.17em}\mathrm{S},\hspace{0.17em}\mathrm{R}\right)$) are associated with the numerical stress integration of the e-HAH model. Therefore, a total of 20 unknown variables (six for the stress tensor, six for the microstructure deviator, and eight state variables) are required to be solved by the stress integration algorithm for continuum elements. For the shell element under the plane stress assumption, the number of unknown variables is reduced to 14. Moreover, the state variables are cross related since the investigated distortional yield surface evolves as a function of the equivalent plastic strain. Therefore, if all the state variables are treated as independent, a total of 20 (or 14) nonlinear equations should be simultaneously solved for continuum (or shell) elements. However, if only a portion of the state variables is assumed to be independent, the exact evolution of the dependent variables should be included in the residuals of the independent variables. For the e-HAH model, this is not an easy task since the evolution laws cannot be explicitly expressed as functions of other independent variables. In this case, an additional iterative process is then required to determine the state variables in the stress update algorithms [55].

Another challenge in the computational modeling of e-HAH plasticity is represented by the continuous and anisotropic distortions of the yield surface during plastic deformation. The distortional characteristics may lead to abrupt changes in the first and/or second gradients of the yield function, which are also complex functions of the plastic deformation. The gradients are key factors in the flow rule of the elastic-plasticity theory, and the accurate numerical integration algorithm implemented into the Euler backward scheme is essential for robust FE simulations. Due to these difficulties, previous studies employed numerically calculated derivatives based on the finite difference method [57,58,60].

Motivated from the above issues in the FE modeling for the distortional e-HAH model, in-depth comparative studies on stress-integration algorithms were carried out in this work. For this purpose, numerical investigations on two common algorithms were considered, namely the CPM and EBM schemes. Only the plane stress condition was studied, but the overall approach can be directly extended to the general stress state. Moreover, the EBM algorithm was formulated with different numbers of residuals for updating the state variables: EBM with 4 residuals (EBM-4R) and EBM with 14 residuals (EBM-14R). The EBM-4R algorithm includes residuals for three stress components and an equivalent plastic strain, whereas the EBM-14R algorithm uses additional 10 residuals for three microstructure deviator components and seven state variables on the yield function distortions. In order to study the algorithm effect for calculating the derivatives of the evolving yield surface, two different methods for updating the yield surface gradients were then considered for the EBM-14R algorithm: analytically derived derivatives (EBM-14R/AD) and numerical derivatives based on finite difference (EBM-14R/ND). Note that only the analytical derivatives were used for the EBM-4R algorithm. In summary, four different stress integration algorithms were comparatively studied, i.e., CPM, EBM-4R, EBM-14R/AD, and EBM-14R/ND. In Figure 1, a schematic interpretation of the four different algorithms is presented. In Figure 1c, it can be seen that the directions normal to the yield surface are not identical between the EBM algorithms with analytical and numerical derivatives. This is a consequence of the truncation error that takes place from the finite difference used for the numerical derivatives. The detailed derivations for the analytical and numerical derivatives of the e-HAH model are presented in Appendix A and Appendix B.

3.2. Stress Update Algorithms for Elastic-Plasticity

In the implementation of small strain elastic-plasticity, the total strain increment, $\Delta \mathsf{\epsilon }$, is additively decomposed into its elastic component, $\Delta {\mathsf{\epsilon }}^{\mathrm{e}}$, and plastic component, $\Delta {\mathsf{\epsilon }}_{}^{\mathrm{p}}$.

$$\Delta \mathsf{\epsilon }=\Delta {\mathsf{\epsilon }}^{\mathrm{e}}+\Delta {\mathsf{\epsilon }}^{\mathrm{p}}.$$

(15)

If the associated flow rule is applied, $\Delta {\mathsf{\epsilon }}_{}^{\mathrm{p}}$ is calculated using the gradient of the yield function,$\mathsf{\Phi }$, and a scalar multiplier, $\Delta \mathsf{\gamma }$.

$$\Delta {\mathsf{\epsilon }}^{\mathrm{p}}=\Delta \mathsf{\gamma }\hspace{0.17em}\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}.$$

(16)

The stress increment, $\Delta \mathsf{\sigma }$, can be expressed as for the isotropic linear elastic metals:

$$\Delta \mathsf{\sigma }=\mathbf{C}\hspace{0.17em}:\left(\Delta \mathsf{\epsilon }-\Delta {\mathsf{\epsilon }}^{\mathrm{p}}\right)=\mathbf{C}\hspace{0.17em}:\left(\Delta \mathsf{\epsilon }-\Delta \mathsf{\gamma }\hspace{0.17em}\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}\right),$$

(17)

where $\mathbf{C}$ is the fourth-order elastic stiffness matrix. The relation between the plastic multiplier, $\Delta \mathsf{\gamma }$, and the equivalent plastic strain increment, $\Delta \overline{\mathsf{\epsilon }}$, is obtained from the 1st order homogeneous function and plastic work equivalence principle.

$$\Delta \overline{\mathsf{\epsilon }}=\frac{\mathsf{\sigma }:\Delta {\mathsf{\epsilon }}^{\mathrm{p}}}{\mathsf{\Phi }}=\frac{\mathsf{\sigma }:\Delta \mathsf{\gamma }\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}}{\mathsf{\Phi }}=\hspace{0.17em}\Delta \mathsf{\gamma }.$$

(18)

Therefore, it is noted that the stress integration in Equation (17) is used to update the equivalent plastic strain increment in Equation (18). Depending on the algorithmic treatment in the above general procedure, various algorithmic procedures have been proposed for numerically accurate and efficient integration methods.

In general, for most constitutive laws and their corresponding numerical implementations, the following predictor-corrector algorithm is employed:

$${\mathsf{\sigma }}_{\mathrm{n}+1}^{\mathrm{T}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}={\mathsf{\sigma }}_{\mathrm{n}}+\mathbf{C}\hspace{0.17em}:\Delta {\mathsf{\epsilon }}_{\mathrm{n}+1},$$

(19)

$$\mathrm{R}=\mathsf{\Phi }\left({\mathsf{\sigma }}_{\mathrm{n}+1}^{\mathrm{T}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}},\hspace{0.17em}{\mathrm{g}}_{\mathrm{i},\mathrm{n}},\hspace{0.17em}{\widehat{h}}_{\mathrm{n}}^{\mathbf{s}}\right)-\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}}\right)<0\hspace{0.17em},\hspace{0.17em}\left(\mathrm{i}=1~4,\hspace{0.17em}\mathrm{L},\hspace{0.17em}\mathrm{S},\hspace{0.17em}\mathrm{R}\right),$$

(20)

where $\Delta {\mathsf{\epsilon }}_{\mathrm{n}+1}$ is the total strain increment at the current time step n + 1, and ${\mathsf{\sigma }}_{\mathrm{n}+1}^{\mathrm{T}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}$ is the trial stress as a predictor. The numerical algorithm starts by checking the condition in Equation (20) as an elastic process, and an iterative process is followed if the condition in Equation (20) is not satisfied. Note that Equation (20) is simplified to represent the yield function $\mathsf{\Phi }$ as a function of the state variables in the e-HAH model.

In the next sub-sections, existing numerical algorithms for the stress integration of the e-HAH model are firstly introduced, followed by the algorithmic development of the proposed modeling.

3.3. Stress Integration Algorithms for the e-HAH Model: Review of Existing Models

3.3.1. Cutting-Plane Method (CPM)

Lee et al. presented the first numerical algorithm for the original HAH model [37,55] and e-HAH model [56] based on the CPM and e-HAH model, respectively. In their CPM approach, the equivalent plastic strain, $\overline{\mathsf{\epsilon }}$, is the only independent variable. For each iteration of the current time step n + 1, $\mathsf{\delta }\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k}+1)}$ is calculated from Equations (16)–(18), using linearized Equation (20).

$${\mathrm{R}}_{\mathrm{n}+1}^{(\mathrm{k})}+{\left(\frac{\partial \mathrm{R}}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}\mathsf{\delta }\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k}+1)}\right)=0\hspace{0.17em},$$

(21)

$${\left(\frac{\partial \mathrm{R}}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}={\left(\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}:\frac{\partial \mathsf{\sigma }}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}-{\left(\frac{\partial \overline{\mathsf{\sigma }}}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})},$$

(22)

$${\left(\frac{\partial \mathsf{\sigma }}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}=-\mathbf{C}:{\left(\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}\right)}_{\mathrm{n}+1}^{(\mathrm{k})},$$

(23)

$$\mathsf{\delta }{\left(\Delta \overline{\mathsf{\epsilon }}\right)}_{\mathrm{n}+1}^{(\mathrm{k}+1)}=\hspace{0.17em}\frac{{\mathrm{R}}_{\mathrm{n}+1}^{(\mathrm{k})}}{{\left(\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}:\mathbf{C}:\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}+{\left(\frac{\partial \overline{\mathsf{\sigma }}}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}}.$$

(24)

The state variables and stress tensor for the k+1 iteration can then be explicitly updated using the value in Equation (24) as follows:

$$\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k}+1)}=\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}+\mathsf{\delta }\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k}+1)}\right),$$

(25)

$${\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}^{(\mathrm{k}+1)}={\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}^{(\mathrm{k})}+\mathsf{\delta }{\left({\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}\right)}^{(\mathrm{k})}={\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}^{(\mathrm{k})}+{\left(\frac{\mathrm{d}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)}^{(\mathrm{k})}\cdot \mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)},\hspace{0.17em}\left(\mathrm{i}=1~4,\hspace{0.17em}\mathrm{L},\hspace{0.17em}\mathrm{S},\hspace{0.17em}\mathrm{R}\right),$$

(26)

$${\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}^{(\mathrm{k}+1)}={\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}^{(\mathrm{k})}+\mathsf{\delta }{\left({\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}\right)}^{(\mathrm{k})}={\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}^{(\mathrm{k})}+{\left(\frac{\mathrm{d}{\widehat{h}}_{\mathrm{n}+1}^{\mathbf{s}}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)}^{(\mathrm{k})}\cdot \mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)},$$

(27)

$$\Delta {\mathsf{\sigma }}_{\mathrm{n}+1}^{(\mathrm{k}+1)}=\mathbf{C}\hspace{0.17em}:\left(\Delta {\mathsf{\epsilon }}_{\mathrm{n}+1}-\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k}+1)}\right),$$

(28)

The iterative procedure continues until the following condition is satisfied:

$$\mathsf{\Phi }\left({\mathsf{\sigma }}_{\mathrm{n}+1}^{(\mathrm{k}+1)},\hspace{0.17em}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}^{(\mathrm{k}+1)},\hspace{0.17em}{\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}^{(\mathrm{k}+1)}\right)-\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k}+1)}\right)<\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}}^{}\right)\cdot \mathrm{Tol}\hspace{0.17em},\hspace{0.17em}\left(\mathrm{Tol}={10}^{-6}\right)\hspace{0.17em}.$$

(29)

Therefore, in the CPM approach, $\overline{\mathsf{\epsilon }}$ is the only unknown variable, and all other variables are treated as a function of $\mathsf{\sigma }$. From an algorithmic point of view, this leads to the yield surface gradient (in Equation (24)) being expressed as a complex form:

$$\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}={\mathsf{\Phi }}^{1-\mathrm{q}}\left[\begin{array}{l}{\left(\mathsf{\psi }(\mathbf{s})+\mathsf{\psi }({\mathbf{s}}_{\mathrm{p}})\right)}^{\frac{\mathrm{q}}{2}-1}\left\{\mathsf{\psi }(\mathbf{s})\frac{\partial \mathsf{\psi }(\mathbf{s})}{\partial \mathsf{\sigma }}+\mathsf{\psi }({\mathbf{s}}_{\mathrm{p}})\frac{\partial \mathsf{\psi }({\mathbf{s}}_{\mathrm{p}})}{\partial \mathsf{\sigma }}\right\}\\ +\frac{\partial {\mathrm{f}}_2}{\partial {\mathrm{g}}_2}\cdot \frac{\partial {\mathrm{g}}_2}{\partial \overline{\mathsf{\epsilon }}}\cdot \frac{\partial \overline{\mathsf{\epsilon }}}{\partial \mathsf{\sigma }}\cdot {\mathrm{f}}_2^{\mathrm{q}-1}{(2{\widehat{h}}^{\mathbf{s}}:\mathbf{s})}^{\mathrm{q}}\\ +2\cdot {\mathrm{f}}_2^{\mathrm{q}}\cdot {(2{\widehat{h}}^{\mathbf{s}}:\mathbf{s})}^{\mathrm{q}-1}\left\{s:\left(\frac{\partial {\widehat{h}}^{\mathbf{s}}}{\partial \overline{\mathsf{\epsilon }}}\otimes \frac{\partial \overline{\mathsf{\epsilon }}}{\partial \mathsf{\sigma }}\right)+{\widehat{h}}^{\mathbf{s}}:\frac{\partial \mathbf{s}}{\partial \mathsf{\sigma }}\right\}\end{array}\right],\hspace{0.17em}\left(\mathrm{when}\hspace{0.17em}{\widehat{h}}^{\mathbf{s}}:\mathbf{s}>0\right).$$

(30)

Lee et al. [56] employed a Newton-Raphson procedure as an alternative method to the analytical derivation of the gradient, and an additional sub-stepping method was also used for numerical stability. Moreover, they suggested a further simplification of Equation (30) by ignoring higher-order terms as follows:

$$\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}={\mathsf{\Phi }}^{1-\mathrm{q}}\left[\begin{array}{l}{\left(\mathsf{\psi }(\mathbf{s})+\mathsf{\psi }({\mathbf{s}}_{\mathrm{p}})\right)}^{\frac{\mathrm{q}}{2}-1}\left\{\mathsf{\psi }(\mathbf{s})\frac{\partial \mathsf{\psi }(\mathbf{s})}{\partial \mathsf{\sigma }}+\mathsf{\psi }({\mathbf{s}}_{\mathrm{p}})\frac{\partial \mathsf{\psi }({\mathbf{s}}_{\mathrm{p}})}{\partial \mathsf{\sigma }}\right\}\\ +2\cdot {\mathrm{f}}_2^{\mathrm{q}}\cdot {(2{\widehat{h}}^{\mathbf{s}}:\mathbf{s})}^{\mathrm{q}-1}{\widehat{h}}^{\mathbf{s}}:\frac{\partial \mathbf{s}}{\partial \mathsf{\sigma }}\end{array}\right],\hspace{0.17em}\left(\mathrm{when}\hspace{0.17em}{\widehat{h}}^{\mathbf{s}}:\mathbf{s}>0\right).$$

(31)

In the implicit FE formulation, an algorithmic tangent modulus is commonly required to obtain a quadratic convergence rate. However, in the conventional CPM—based stress-integration algorithm, this cannot be readily applied. For example, an additional numerical technique was proposed for the tangent moduli to preserve the quadratic convergence rate [54], but this could not be applied to the distortional hardening models. For the CPM of the HAH or e-HAH approaches presented in [55,56], the continuum tangent modulus, ${\mathbf{C}}^{\mathrm{e}\mathrm{p}}$, was used as an alternative to the consistent tangent modulus:

$${\mathbf{C}}^{\mathrm{e}\mathrm{p}}=\frac{\mathrm{d}{\mathsf{\sigma }}_{\mathrm{n}+1}}{\mathrm{d}{\mathsf{\epsilon }}_{\mathrm{n}+1}}=\mathbf{C}-\left(\frac{\mathbf{C}:\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\mathsf{\sigma }}_{\mathrm{n}+1}}\otimes \mathbf{C}:\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\mathsf{\sigma }}_{\mathrm{n}+1}}}{\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\mathsf{\sigma }}_{\mathrm{n}+1}}:\mathbf{C}:\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\mathsf{\sigma }}_{\mathrm{n}+1}}+\frac{\partial {\overline{\mathsf{\sigma }}}_{\mathrm{n}+1}}{\partial {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}}}\right).$$

(32)

3.3.2. Euler Backward Method (EBM)

Previous works have focused on the development of the EBM—based stress integration algorithm of the distortional anisotropic hardening laws. However, these algorithms simplified the set of equations in the EBM scheme by reducing the number of residuals associated with the state variables. Residuals on the stress tensor and equivalent plastic strain were then required to satisfy the consistency condition in elastic-plasticity under the associated flow rule. However, this method is affected by similar problems as the CPM—based algorithm [49] since a portion of the state variables should be treated as functions of independent variables when the number of residuals is less than that of the independent variables in the e-HAH model. This may cause difficulty in calculating the first and second derivatives of the yield function during its evolution as a function of plastic deformation.

The EBM algorithm developed by Lee et al. [56] introduced the following residuals:

$${\mathrm{R}}_{1,\mathrm{n}+1}^{(\mathrm{k})}=\mathsf{\Phi }\left({\mathsf{\sigma }}_{\mathrm{n}+1}^{(\mathrm{k})},\hspace{0.17em}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}^{(\mathrm{k})},\hspace{0.17em}{\widehat{h}}_{\mathrm{i},\mathrm{n}+1}^{\mathbf{s}(\mathrm{k})}\right)-\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}\right)\hspace{0.17em},\hspace{0.17em}\left(\mathrm{i}=1-7,\hspace{0.17em}5\hspace{0.17em}=\mathrm{L},\hspace{0.17em}6=\mathrm{S},\hspace{0.17em}7=\mathrm{R}\right),$$

(33)

$${\mathbf{R}}_{2,\mathrm{n}+1}^{(\mathrm{k})}=-\Delta {\mathsf{\epsilon }}_{\mathrm{n}+1}^{\mathrm{p}(\mathrm{k})}+\Delta \mathsf{\gamma }\hspace{0.17em}{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}^{(\mathrm{k})}$$

(34)

To apply the multi-variable Newton-Raphson method, the residuals are linearized according to the following:

$${\mathrm{R}}_{1,\mathrm{n}+1}^{(\mathrm{k})}-{\partial }_{\overline{\mathsf{\epsilon }}}\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}\right)\cdot \mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}+{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}^{(\mathrm{k})}:\mathsf{\delta }{\left(\Delta {\mathsf{\sigma }}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}=0\hspace{0.17em},$$

(35)

$$\begin{array}{c}{\mathbf{R}}_{2,\mathrm{n}+1}^{(\mathrm{k})}+{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}^{(\mathrm{k})}\cdot \mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}+{\mathsf{\Xi }}^{‒1}{}_{\mathrm{n}+1}^{(\mathrm{k})}:\mathsf{\delta }{\left(\Delta {\mathsf{\sigma }}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}=0.\\ {\mathsf{\Xi }}^{‒1}{}_{\mathrm{n}+1}^{(\mathrm{k})}={\left({\mathbf{C}}^{‒1}+\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}\hspace{0.17em}{\partial }_{\mathsf{\sigma }\mathsf{\sigma }}^2{\mathsf{\Phi }}_{\mathrm{n}+1}^{(\mathrm{k})}\hspace{0.17em}\right)}^{‒1}\end{array}$$

(36)

The linearized equations can be calculated with respect to $\mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k})}$ and $\mathsf{\delta }{\left(\Delta {\mathsf{\sigma }}_{\mathrm{n}+1}\right)}^{(\mathrm{k})}$ by solving the following system:

$$\left\{\begin{array}{c}\mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\\ \mathsf{\delta }{\left(\Delta {\mathsf{\sigma }}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\end{array}\right\}=-{\left[\begin{array}{cc}{\left(\frac{\partial {\mathrm{R}}_1}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathrm{R}}_1}{\partial \Delta \mathsf{\sigma }}\right)}_{\mathrm{n}+1}^{\mathrm{k}}\\ {\left(\frac{\partial {\mathbf{R}}_2}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathbf{R}}_2}{\partial \Delta \mathsf{\sigma }}\right)}_{\mathrm{n}+1}^{\mathrm{k}}\end{array}\right]}^{-1}\cdot \left\{\begin{array}{c}{\mathrm{R}}_{1,\mathrm{n}+1}^{(\mathrm{k})}\\ {\mathbf{R}}_{2,\mathrm{n}+1}^{(\mathrm{k})}\end{array}\right\}.$$

(37)

Upon reaching convergence for the solutions of Equation (37), the other variables are updated explicitly using Equations (25)–(27). The iterative process stops when the following conditions are satisfied:

$$\Vert {\mathrm{R}}_{1+\mathrm{i},\mathrm{n}+1}^{(\mathrm{k})}\Vert <\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}}\right)\cdot \hspace{0.17em}\mathrm{T}\mathrm{o}\mathrm{l},\hspace{0.17em}\Vert {\mathbf{R}}_{2,\mathrm{n}+1}^{(\mathrm{k})}\Vert <\mathrm{T}\mathrm{o}\mathrm{l}\hspace{0.17em}.$$

(38)

As previously mentioned, the major difficulty in this algorithmic procedure lies in the calculation of the first and second derivatives, due to the anisotropically evolving yield function. Additionally, this evolution is closely associated with a multiple number of dependent variables in the e-HAH model. To overcome this issue, Lee et al. [55] employed the multi-step Newton-Raphson method by subdividing the variables $\hspace{0.17em}{\mathrm{g}}_{1\mathrm{o}\mathrm{r}2}$. Recently, Choi and Yoon [57] also applied the multi-stage EBM algorithm by using the sub-division of the strain increment, $\Delta {\mathsf{\epsilon }}_{\mathrm{n}+1}$, which was originally proposed by Yoon et al. [61]. Furthermore, Yoon et al. [58] used a line-search method for the step size control in Equation (36) as an alternative approach to increase the algorithm stability for the distortional hardening model.

Contrary to the CPM algorithm, the EBM—based stress integrations have a quadratic convergence, which is consistent with the global Newton-Raphson method. The consistent tangent modulus can be expressed as follows:

$${\mathbf{C}}^{\mathrm{e}\mathrm{p}}={\mathsf{\Xi }}_{\mathrm{n}+1}-\left(\frac{{\mathsf{\Xi }}_{\mathrm{n}+1}:{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}\otimes {\mathsf{\Xi }}_{\mathrm{n}+1}:{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}}{{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}:{\mathsf{\Xi }}_{\mathrm{n}+1}:{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}+{\partial }_{\overline{\mathsf{\epsilon }}}\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}\right).$$

(39)

In this Equation, the tangent modulus requires the calculations of the first and second derivatives of the yield function under the associated flow rule. As previously discussed, the existing EBM—based algorithms introduced simpler forms without presenting complete derivations of their analytical set.

3.4. Stress Integration Algorithms for the e-HAH Model: Proposed Models

3.4.1. Cutting-Plane Method (CPM)

In comparison with the CPM algorithm used by Lee et al. [55], whereby the effect of the dependent variables is included in the yield surface gradient, the present CPM scheme introduces the linearization of the yield condition in Equation (20) as follows:

$${\left(\frac{\partial \mathrm{R}}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}={\left(\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}:\frac{\partial \mathsf{\sigma }}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}+{\left(\frac{\partial \mathsf{\Phi }}{\partial {\widehat{h}}^{\mathbf{s}}}:\frac{\partial {\widehat{h}}^{\mathbf{s}}}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}+{\left({\displaystyle \sum _{\mathrm{i}}\frac{\partial \mathsf{\Phi }}{\partial {\mathrm{g}}_{\mathrm{i}}}\cdot \frac{\partial {\mathrm{g}}_{\mathrm{i}}}{\partial \Delta \overline{\mathsf{\epsilon }}}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}-{\left(\frac{\partial \overline{\mathsf{\sigma }}}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}.$$

(40)

After substituting Equations (40) and (23) into Equation (21), $\mathsf{\delta }{\left(\Delta \overline{\mathsf{\epsilon }}\right)}_{\mathrm{n}+1}^{(\mathrm{k}+1)}$ can be calculated:

$$\mathsf{\delta }{\left(\Delta \overline{\mathsf{\epsilon }}\right)}_{\mathrm{n}+1}^{(\mathrm{k}+1)}=\hspace{0.17em}\frac{{\mathrm{R}}_{\mathrm{n}+1}^{(\mathrm{k})}}{{\left(\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}:\mathbf{C}:\frac{\partial \mathsf{\Phi }}{\partial \mathsf{\sigma }}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}-{\left(\frac{\partial \mathsf{\Phi }}{\partial {\widehat{h}}^{\mathbf{s}}}:\frac{\partial {\widehat{h}}^{\mathbf{s}}}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}-{\left({\displaystyle \sum _{\mathrm{i}}\frac{\partial \mathsf{\Phi }}{\partial {\mathrm{g}}_{\mathrm{i}}}\cdot \frac{\partial {\mathrm{g}}_{\mathrm{i}}}{\partial \Delta \overline{\mathsf{\epsilon }}}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}+{\left(\frac{\partial \overline{\mathsf{\sigma }}}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{(\mathrm{k})}}.$$

(41)

The state variables and stress tensor for the $(\mathrm{k}+1)$th iteration at the current time step are iteratively updated using Equations (25)–(28) until Equation (29) is satisfied. The tangent modulus ${\mathbf{C}}^{\mathrm{e}\mathrm{p}}$ is also calculated using the consistency condition:

$$\mathrm{d}{\mathrm{R}}_{\mathrm{n}+1}=\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\mathsf{\sigma }}_{\mathrm{n}+1}}:\mathrm{d}{\mathsf{\sigma }}_{\mathrm{n}+1}+\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}}:\frac{\partial {\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}}{\partial {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}}\mathrm{d}{\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}+\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}\cdot \frac{\partial {\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}{\partial {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}}\mathrm{d}{\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}-\frac{\partial {\overline{\mathsf{\sigma }}}_{\mathrm{n}+1}}{\partial {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}}\mathrm{d}{\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}=0.$$

(42)

By substituting Equation (17) into Equation (42), the elastoplastic tangent modulus is expressed as follows:

$${\mathbf{C}}^{\mathrm{e}\mathrm{p}}=\frac{\mathrm{d}{\mathsf{\sigma }}_{\mathrm{n}+1}}{\mathrm{d}{\mathsf{\epsilon }}_{\mathrm{n}+1}}=\mathbf{C}-\left(\frac{\mathbf{C}:\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\mathsf{\sigma }}_{\mathrm{n}+1}}\otimes \mathbf{C}:\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\mathsf{\sigma }}_{\mathrm{n}+1}}}{\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\mathsf{\sigma }}_{\mathrm{n}+1}}:\mathbf{C}:\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\mathsf{\sigma }}_{\mathrm{n}+1}}-\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}}:\frac{\partial {\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}}{\partial {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}}-\frac{\partial {\mathsf{\Phi }}_{\mathrm{n}+1}}{\partial {\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}\cdot \frac{\partial {\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}{\partial {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}}+\frac{\partial {\overline{\mathsf{\sigma }}}_{\mathrm{n}+1}}{\partial {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}}}\right).$$

(43)

3.4.2. Fully Implicit Euler Backward Method (EBM)

The main characteristic of the e-HAH model is represented by the fact that the state variables are cross-related with each other. The complexity resulting from the constitutive model makes it challenging to construct the exact linearization of the governing equation with a limited number of independent variables. This is why the existing algorithms were developed with their own simplifying methods. In this study, our efforts are focused on deriving a fully implicit EBM algorithm by simultaneously solving the whole set of residuals defined for the stress tensor, microstructure deviator tensor, and eight state variables (or equivalently 14 unknowns) associated with the e-HAH scheme under plane stress conditions. To the best of our knowledge, this is the first time that an algorithmic implementation of the distortional anisotropic hardening e-HAH model is proposed.

In the following, a set of nonlinear equations defining the residuals of the e-HAH model is provided.

$${\mathrm{R}}_{1,\mathrm{n}+1}^{(\mathrm{k})}=\mathsf{\Phi }\left({\mathsf{\sigma }}_{\mathrm{n}+1}^{(\mathrm{k})},\hspace{0.17em}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}^{(\mathrm{k})},\hspace{0.17em}{\widehat{h}}_{\mathrm{i},\mathrm{n}+1}^{\mathbf{s}(\mathrm{k})}\right)-\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}\right)\hspace{0.17em},\hspace{0.17em}\left(\mathrm{i}=1-7,\hspace{0.17em}5\hspace{0.17em}=\mathrm{L},\hspace{0.17em}6=\mathrm{S},\hspace{0.17em}7=\mathrm{R}\right),$$

(44)

$${\mathbf{R}}_{2,\mathrm{n}+1}^{(\mathrm{k})}=-\Delta {\mathsf{\epsilon }}_{\mathrm{n}+1}^{\mathrm{p}(\mathrm{k})}+\Delta \mathsf{\gamma }\hspace{0.17em}{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}^{(\mathrm{k})},$$

(45)

$${\mathrm{R}}_{2+\mathrm{i},\mathrm{n}+1}^{(\mathrm{k})}={\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}^{(\mathrm{k})}-{\mathrm{g}}_{\mathrm{i},\mathrm{n}}-{\left(\frac{\mathrm{d}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)}^{(\mathrm{k})}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}\hspace{0.17em},$$

(46)

$${\mathbf{R}}_{10,\mathrm{n}+1}^{(\mathrm{k})}=\Delta {\widehat{h}}_{\mathrm{n}+1}^{\mathbf{s}(\mathrm{k})}-\Delta {\widehat{h}}_{\mathrm{n}}^{\mathbf{s}}-{\left(\frac{\mathrm{d}{\widehat{h}}_{\mathrm{n}+1}^{\mathbf{s}}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)}^{(\mathrm{k})}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})},$$

(47)

Equation (44) provides the consistency condition, Equation (45) represents the associated flow rule, Equation (46) expresses the state variables controlling the distortion of the yield surface, and Equation (47) is written to define the rotation of the microstructure deviator, ${\widehat{h}}^{\mathbf{s}}$.

Linearization of Equations (44)–(47) leads to the following equations:

$$\begin{array}{l}{\mathrm{R}}_{1,\mathrm{n}+1}^{(\mathrm{k})}-{\partial }_{\overline{\mathsf{\epsilon }}}\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}\right)\cdot \mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}+{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}^{(\mathrm{k})}:\mathsf{\delta }{\left(\Delta {\mathsf{\sigma }}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\\ +{\displaystyle \sum _{\mathrm{i}}{\partial }_{{\mathrm{g}}_{\mathrm{i}}}{\mathrm{f}}_{\mathrm{n}+1}^{(\mathrm{k})}\cdot {\left(\Delta {\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}}+{\partial }_{{\widehat{h}}^{\mathbf{s}}}{\mathsf{\Phi }}_{\mathrm{n}+1}^{(\mathrm{k})}:\mathsf{\delta }{\left(\Delta {\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}=0\hspace{0.17em},\end{array}$$

(48)

$$\begin{array}{c}{\mathbf{R}}_{2,\mathrm{n}+1}^{(\mathrm{k})}+{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}^{(\mathrm{k})}\cdot \mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}+{\mathsf{\Xi }}^{‒1}{}_{\mathrm{n}+1}^{(\mathrm{k})}:\mathsf{\delta }{\left(\Delta {\mathsf{\sigma }}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\\ +{\displaystyle \sum _{{\mathrm{g}}_{\mathrm{i}}}{\partial }_{\mathsf{\sigma }{\mathrm{g}}_{\mathrm{i}}}^2}{\mathsf{\Phi }}_{\mathrm{n}+1}^{(\mathrm{k})}\hspace{0.17em}\cdot \mathsf{\delta }{\left(\Delta {\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}+\hspace{0.17em}{\partial }_{\mathsf{\sigma }{\widehat{h}}^{\mathbf{s}}}^2{\mathsf{\Phi }}_{\mathrm{n}+1}^{(\mathrm{k})}\hspace{0.17em}:\mathsf{\delta }{\left(\Delta {\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}\right)}^{(\mathrm{k})}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k}+1)}=0\\ {\mathsf{\Xi }}^{‒1}{}_{\mathrm{n}+1}^{(\mathrm{k})}={\left({\mathbf{C}}^{‒1}+\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}\hspace{0.17em}{\partial }_{\mathsf{\sigma }\mathsf{\sigma }}^2{\mathsf{\Phi }}_{\mathrm{n}+1}^{(\mathrm{k})}\hspace{0.17em}\right)}^{‒1}\end{array}$$

(49)

$$\begin{array}{c}{\mathrm{R}}_{2+\mathrm{i},\mathrm{n}+1}^{(\mathrm{k})}-\left(\frac{\mathrm{d}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)\cdot \mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}-{\partial }_{\mathsf{\sigma }}\left(\frac{\mathrm{d}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)\cdot \mathsf{\delta }{\left(\Delta {\mathsf{\sigma }}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}\\ -{\partial }_{{\mathrm{g}}_{\mathrm{i}}}\left(\frac{\mathrm{d}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)\cdot \mathsf{\delta }{\left(\Delta {\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}-{\partial }_{{\widehat{h}}^{\mathbf{s}}}\left(\frac{\mathrm{d}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)\cdot \mathsf{\delta }{\left(\Delta {\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}=0\end{array}$$

(50)

$$\begin{array}{c}{\mathbf{R}}_{10,\mathrm{n}+1}^{(\mathrm{k})}+\Delta {\widehat{h}}_{\mathrm{n}+1}^{\mathbf{s}(\mathrm{k})}-{\left(\frac{\mathrm{d}{\widehat{h}}_{\mathrm{n}+1}^{\mathbf{s}}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)}^{(\mathrm{k})}\cdot \mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}-{\partial }_{\mathsf{\sigma }}{\left(\frac{\mathrm{d}{\widehat{h}}_{\mathrm{n}+1}^{\mathbf{s}}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)}^{(\mathrm{k})}:\mathsf{\delta }{\left(\Delta {\mathsf{\sigma }}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}\\ -{\partial }_{{\mathrm{g}}_{\mathrm{i}}}{\left(\frac{\mathrm{d}{\widehat{h}}_{\mathrm{n}+1}^{\mathbf{s}}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)}^{(\mathrm{k})}\cdot \mathsf{\delta }{\left(\Delta {\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}-{\partial }_{{\widehat{h}}^{\mathbf{s}}}{\left(\frac{\mathrm{d}{\widehat{h}}_{\mathrm{n}+1}^{\mathbf{s}}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\right)}^{(\mathrm{k})}:\mathsf{\delta }{\left(\Delta {\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}^{(\mathrm{k})}\hspace{0.17em}=0\end{array}.$$

(51)

The solutions of the above linearized equations are then obtained for $\mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}$, $\mathsf{\delta }{\left(\Delta {\mathsf{\sigma }}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}$, $\mathsf{\delta }{\left(\Delta {\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}$, and $\mathsf{\delta }{\left(\Delta {\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}$ by solving the following matrix equations:

$$\left\{\begin{array}{c}\mathsf{\delta }{\left(\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\\ \mathsf{\delta }{\left(\Delta {\mathsf{\sigma }}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\\ \mathsf{\delta }{\left(\Delta {\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\\ \mathsf{\delta }{\left(\Delta {\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}\right)}^{(\mathrm{k}+1)}\end{array}\right\}=-{\left[\begin{array}{cccc}{\left(\frac{\partial {\mathrm{R}}_1}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathrm{R}}_1}{\partial \Delta \mathsf{\sigma }}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathrm{R}}_1}{\partial \Delta {\mathrm{g}}_{\mathrm{i}}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathrm{R}}_1}{\partial \Delta {\widehat{h}}^{\mathbf{s}}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}\\ {\left(\frac{\partial {\mathbf{R}}_2}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathbf{R}}_2}{\partial \Delta \mathsf{\sigma }}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathbf{R}}_2}{\partial \Delta {\mathrm{g}}_{\mathrm{i}}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathbf{R}}_2}{\partial \Delta {\widehat{h}}^{\mathbf{s}}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}\\ {\left(\frac{\partial {\mathrm{R}}_{2+\mathrm{i}}}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathrm{R}}_{2+\mathrm{i}}}{\partial \Delta \mathsf{\sigma }}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathrm{R}}_{2+\mathrm{i}}}{\partial \Delta {\mathrm{g}}_{\mathrm{i}}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathrm{R}}_{2+\mathrm{i}}}{\partial \Delta {\widehat{h}}^{\mathbf{s}}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}\\ {\left(\frac{\partial {\mathbf{R}}_{10}}{\partial \Delta \overline{\mathsf{\epsilon }}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathbf{R}}_{10}}{\partial \Delta \mathsf{\sigma }}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathbf{R}}_{10}}{\partial \Delta {\mathrm{g}}_{\mathrm{i}}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}& {\left(\frac{\partial {\mathbf{R}}_{10}}{\partial \Delta {\widehat{h}}^{\mathbf{s}}}\right)}_{\mathrm{n}+1}^{\mathrm{k}}\end{array}\right]}^{-1}\cdot \left\{\begin{array}{c}{\mathrm{R}}_{1,\mathrm{n}+1}^{(\mathrm{k})}\\ {\mathbf{R}}_{2,\mathrm{n}+1}^{(\mathrm{k})}\\ {\mathrm{R}}_{2+\mathrm{i},\mathrm{n}+1}^{(\mathrm{k})}\\ {\mathbf{R}}_{10,\mathrm{n}+1}^{(\mathrm{k})}\end{array}\right\}$$

(52)

Equation (52) shows that the Jacobian matrix has a size of 14 × 14 when it is expressed using the Voigt notation under the plane stress condition. By applying Newton’s method, all state variables are updated for the $\left(\mathrm{k}+1\right)$th iteration, and this process continues until the residuals satisfy the following criteria:

$$\begin{array}{l}\Vert {\mathrm{R}}_{1+\mathrm{i},\mathrm{n}+1}^{(\mathrm{k})}\Vert <\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}}\right)\cdot \hspace{0.17em}\mathrm{T}\mathrm{o}\mathrm{l},\hspace{0.17em}\Vert {\mathbf{R}}_{2,\mathrm{n}+1}^{(\mathrm{k})}\Vert <\Vert \Delta {\mathsf{\epsilon }}_{\mathrm{n}+1}\Vert \hspace{0.17em}\cdot \mathrm{T}\mathrm{o}\mathrm{l}\hspace{0.17em},\hspace{0.17em}\\ \Vert {\mathrm{R}}_{2+\mathrm{i},\mathrm{n}+1}^{(\mathrm{k})}\Vert <\mathrm{T}\mathrm{o}\mathrm{l},\mathrm{a}\mathrm{n}\mathrm{d}\hspace{0.17em}\Vert {\mathbf{R}}_{10,\mathrm{n}+1}^{(\mathrm{k})}\Vert <{\mathrm{H}}^{-\frac{1}{2}}\cdot \mathrm{T}\mathrm{o}\mathrm{l}\end{array}$$

(53)

In this study, the tolerance, Tol, is set to be ${10}^{-6}$ for all the residuals with adequate normalizations.

The linearized flow rule in Equation (17) and isotropic linear elasticity in Equation (16) lead to the calculation of the consistent tangent modulus ${\mathbf{C}}^{\mathrm{e}\mathrm{p}}$ for the EBM of the e-HAH model as follows:

$$\mathrm{d}{\mathsf{\sigma }}_{\mathrm{n}+1}=\mathbf{C}:\left(\mathrm{d}{\mathsf{\epsilon }}_{\mathrm{n}+1}-\mathrm{d}{\mathsf{\epsilon }}_{\mathrm{n}+1}^{\mathrm{p}}\right),$$

(54)

$$\begin{array}{c}\mathrm{d}{\mathsf{\epsilon }}_{\mathrm{n}+1}^{\mathrm{p}}=\mathrm{d}\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\cdot {\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}+\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\cdot {\partial }_{\mathsf{\sigma }\mathsf{\sigma }}^2{\mathsf{\Phi }}_{\mathrm{n}+1}:\mathrm{d}{\mathsf{\sigma }}_{\mathrm{n}+1}\\ +\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\cdot {\partial }_{\mathsf{\sigma }{\mathrm{g}}_{\mathrm{i}}}^2{\mathsf{\Phi }}_{\mathrm{n}+1}\cdot \mathrm{d}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}+\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\cdot {\partial }_{\mathsf{\sigma }{\widehat{h}}^{\mathbf{s}}}^2{\mathsf{\Phi }}_{\mathrm{n}+1}:\mathrm{d}{\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}\end{array}.$$

(55)

Substituting Equation (55) into Equation (54), one obtains the following:

$$\begin{array}{c}\mathrm{d}{\mathsf{\sigma }}_{\mathrm{n}+1}={\mathsf{\Xi }}_{\mathrm{n}+1}:\left(\mathrm{d}{\mathsf{\epsilon }}_{\mathrm{n}+1}-\mathrm{d}\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\cdot {\mathbf{A}}_{\mathrm{n}+1}\right)\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\\ {\mathbf{A}}_{\mathrm{n}+1}=\left({\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}+{\partial }_{\mathsf{\sigma }{\mathrm{g}}_{\mathrm{i}}}^2{\mathsf{\Phi }}_{\mathrm{n}+1}\cdot \frac{\mathrm{d}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}+{\partial }_{\mathsf{\sigma }{\widehat{h}}^{\mathbf{s}}}^2{\mathsf{\Phi }}_{\mathrm{n}+1}:\frac{\mathrm{d}{\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}\cdot \Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)\end{array}.$$

(56)

The consistency condition with Equation (54) yields the following:

$$\mathrm{d}\Delta {\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}=\frac{{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}:{\mathsf{\Xi }}_{\mathrm{n}+1}:\mathrm{d}{\mathsf{\epsilon }}_{\mathrm{n}+1}}{{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}:{\mathsf{\Xi }}_{\mathrm{n}+1}:{\mathbf{A}}_{\mathrm{n}+1}-{\partial }_{{\mathrm{g}}_{\mathrm{i}}}{\mathsf{\Phi }}_{\mathrm{n}+1}\cdot \frac{\mathrm{d}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}-{\partial }_{{\widehat{h}}^{\mathbf{s}}}{\mathsf{\Phi }}_{\mathrm{n}+1}:\frac{\mathrm{d}{\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}+{\partial }_{\overline{\mathsf{\epsilon }}}\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}.$$

(57)

Finally, ${\mathbf{C}}^{\mathrm{e}\mathrm{p}}$ can be calculated from Equations (56) and (57) according to the following:

$${\mathbf{C}}^{\mathrm{e}\mathrm{p}}={\mathsf{\Xi }}_{\mathrm{n}+1}-\left(\frac{{\mathsf{\Xi }}_{\mathrm{n}+1}:{\mathbf{A}}_{\mathrm{n}+1}\otimes {\mathsf{\Xi }}_{\mathrm{n}+1}:{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}}{{\partial }_{\mathsf{\sigma }}{\mathsf{\Phi }}_{\mathrm{n}+1}:{\mathsf{\Xi }}_{\mathrm{n}+1}:{\mathbf{A}}_{\mathrm{n}+1}-{\partial }_{{\mathrm{g}}_{\mathrm{i}}}{\mathsf{\Phi }}_{\mathrm{n}+1}\cdot \frac{\mathrm{d}{\mathrm{g}}_{\mathrm{i},\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}-{\partial }_{{\widehat{h}}^{\mathbf{s}}}{\mathsf{\Phi }}_{\mathrm{n}+1}:\frac{\mathrm{d}{\widehat{h}}^{\mathbf{s}}{}_{\mathrm{n}+1}}{\mathrm{d}\overline{\mathsf{\epsilon }}}+{\partial }_{\overline{\mathsf{\epsilon }}}\overline{\mathsf{\sigma }}\left({\overline{\mathsf{\epsilon }}}_{\mathrm{n}+1}\right)}\right).$$

(58)

It is important to note that the previously formulated EBM—based algorithms with four residuals [55,56,58] can be retrieved if the residuals in the current formulation are properly reduced. Moreover, the numerical algorithm can be simply extended to a general stress state. The flow charts of the proposed, fully implicit EBM and CPM algorithms are shown in Figure 2.

4. Evaluation of Stress Integration Algorithms for the e-HAH Model

In this section, the numerical accuracy and stability of the proposed stress integration algorithms are analyzed for different case problems. The algorithms described in the previous sections were implemented into the static implicit FE ABAQUS/Standard software using the UMAT user material subroutine. For the comparative study, the four algorithms presented in Section 3.1 were implemented into the FE model. Note that the EBM-14R approach represents the fully implicit Euler backward algorithm implemented into the e-HAH model for the first time.

Two levels of the evaluation procedure are presented in the following sections. The first level consists of a very detailed fundamental analysis on the accuracy of the different algorithms for the investigated e-HAH model. For this, the anisotropic characteristics of the flow stress and r-value are predicted and evaluated under two different loading paths. Only a single element is used in this level to rule out other numerical effects. To clarify the effect of the anisotropic hardening responses of the e-HAH model under different loading paths, different model materials were selectively compared. The investigated model materials exhibit high or low evolution rates under reverse loading and hardening or softening under cross-loading. To represent the initial anisotropy of the material, the non-quadratic anisotropic yield function, Yld2000-2d [62,63], was employed. For the summary of the initial (undistorted) yield function, Yld2000-2d, the reader is referred to in Appendix C.

The second level comprises a real-scale simulation based on the S-rail forming and springback process [45], which was proposed as a benchmark problem and often utilized for the analysis of the constitutive model and numerical algorithm in the sheet metal forming community. The reason for choosing a benchmark is that most of the existing FE simulations are based on hybrid explicit and implicit algorithms. In other words, the forming process is solved via the dynamic explicit FE model as a quasi-static problem, whereas the springback is calculated using the static implicit algorithm. The first is applied to avoid divergence issues typically encountered in contact problems between complex tools and sheet metal, whereas the implicit algorithm is the optimal approach for the springback unloading process to reduce the computational time. In this study, all the stress integration algorithms were implemented into the static implicit FE solver, and the numerical performance of each algorithm was comparatively studied to assess their applicability to industrial-size problems. Two materials, namely stainless steel (STS) and dual-phase steel (DP), were investigated since they have both been used for real automotive parts and exhibit distinctive anisotropic hardening behaviors under loading path changes.

4.1. One Element Analysis

In the one-element analysis, the loading condition consists of a compression-tension (C-T) test. The accuracies of the predicted flow stress and r-value were compared for the investigated numerical algorithms. Two loading paths were simulated for the C-T test. The first path consists of a 5% compression along the rolling direction (RD) followed by a 10% tension along the same direction. This case is denoted as “C5T10R,” where “R” represents the fact that the loading path is a “Reversed” loading condition. During the compression along the RD, the rate of the yield function distortion represented by the transient behavior and Bauschinger effect is maximal in the opposite side of the compressive loading, or tension, along the RD.

The second loading path consists of a 5% compression along the RD followed by a 10% tension at an angle of 54° to the RD. Note that the 54° angle corresponds to the cross-loading state, where the condition ${\widehat{h}}^{\mathbf{s}}:\widehat{s}=0$ is satisfied (the exact value of the angle for the cross-loading condition is 54.74°). This case is denoted as “C5T10CR,” where “CR” stands for “cross” loading. This loading path was selected because the cross-loading effect with rotation of the microstructure deviator becomes maximized at the loading direction with the condition.

A four-node shell element with reduced integration in the ABAQUS/Standard, also denoted S4R, was used. The boundary conditions for the above cases are schematically shown in Figure 3.

The r-value along the θ angle with respect to the RD is defined in Equation (59) and calculated via the nodal displacements of an element, according to the following:

$${\mathrm{r}}_{\mathsf{\theta }}=\frac{{\mathsf{\epsilon }}_{\mathrm{w}\mathrm{i}\mathrm{d}\mathrm{t}\mathrm{h}}^{\mathrm{p}}}{{\mathsf{\epsilon }}_{\mathrm{t}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{k}\mathrm{n}\mathrm{e}\mathbf{s}\mathbf{s}}^{\mathrm{p}}}=-\frac{{\mathsf{\epsilon }}_2^{\mathrm{p}}}{{\mathsf{\epsilon }}_1^{\mathrm{p}}+{\mathsf{\epsilon }}_2^{\mathrm{p}}},\hspace{0.17em}\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\hspace{0.17em}\left\{{\mathsf{\epsilon }}_1^{\mathrm{p}}\hspace{0.17em}=\ln \left(1+\Delta {\mathrm{u}}_1^{\mathrm{p}}\right)\hspace{0.17em}\mathrm{a}\mathrm{n}\mathrm{d}\hspace{0.17em}{\mathsf{\epsilon }}_2^{\mathrm{p}}=\hspace{0.17em}\ln \left(1+\Delta {\mathrm{u}}_2^{\mathrm{p}}\right)\right\},$$

(59)

where $\Delta {\mathrm{u}}_1^{\mathrm{p}}$ and $\Delta {\mathrm{u}}_2^{\mathrm{p}}$ are the element nodal displacements along the RD and transverse direction, respectively, during tension, and the superscript “$\mathrm{p}$” represents the plastic part of the strain or displacement.

The mechanical properties of the two model materials and the constitutive parameters for the initial yield function and e-HAH model are listed in

Table 1. Anisotropic mechanical properties. The subscripts 0, 45, 90 represent rolling direction (RD), diagonal direction (DD), and transverse direction (TD), respectively.

Material	σ0/σ0	σ45/σ0	σ90/σ0	σb/σ0
Mat 1,2	1	1	1	1
	r0	r45	r90	rb
	0.5	1.0	1.5	1.0

and

Table 2. Constitutive parameters of isotropic linear elasticity and e-HAH model.

Elastic property and isotropic hardening law parameters
Material	Young’s modulus (GPa)			Poisson’s ratio			K (MPa)			e0	n
MAT 1,2	70.0			0.33			500.0			0.01	0.2
Yld2000-2d coefficients
Material	m	α1		α2		α3	α4	α5	α6	α7	α8
MAT 1,2	8.0	0.7826		1.1778		1.1075	0.9753	1.0154	0.9028	0.9989	0.9990
e-HAH model parameters
Material	q	k	k1	k2	k3	k4	k5	L	kL	S	kS
MAT 1	2.0	250.0	300.0	250.0	0.25	0.80	50	1.0	0.0	0.80	75.0
MAT 2	2.0	25.0	60	50.0	0.90	0.80	50	1.65	384.0	1.0	0.0

, respectively. The model materials “MAT 1” and “MAT 2” have distinctive e-HAH parameters, but other properties, such as elasticity, initial yield function, and isotropic hardening, are identical. For isotropic hardening, the Swift hardening power law, $\overline{\mathsf{\sigma }}\left(\overline{\mathsf{\epsilon }}\right)=\mathrm{K}{\left({\mathrm{e}}_0+\overline{\mathsf{\epsilon }}\right)}^{\mathrm{n}}$, is applied for both materials.

Figure 4 and Figure 5 show the reference flow stress curves and the corresponding evolution of the e-HAH surfaces for MAT 1 and MAT 2, respectively. Since the exact solutions are not available, a reference curve with a sufficiently small time-step is assumed to be exact. The reference flow curves were obtained using the CPM scheme implemented into the dynamic explicit FE ABAQUS/Explicit software with the VUMAT user material subroutine. The average strain increment was set to $2.4\times {10}^{-8}$. Figure 4a and Figure 5a illustrate the flow stress curves under the C5T10R loading path, whereas Figure 4c and Figure 5c display the same curves under the C5T10CR loading path. The evolutions of the yield surfaces are provided in Figure 4b and Figure 5b for C5T10R, and in Figure 4d and Figure 5d for C510CR. For comparison, the evolutions of the isotropic yield surfaces are also included in the figures. In each figure, three points are marked, namely A, B, and C, where A and C represent the initial and final stress states during the second loading, respectively, whereas B is selected between the two loading points to compare the transient behavior in the second loading path.

Figure 4a,b shows that a significant amount of the Bauschinger effect and permanent softening (${\mathrm{g}}_4\cdot \overline{\mathsf{\sigma }}=0.8224\cdot \overline{\mathsf{\sigma }}$) are represented from the e-HAH parameters of MAT 1 under C5T10R. For C5T10CR, a similar transient behavior and permanent softening were calculated, but the Bauschinger effect (contraction) was found to be much less than for the C5T10R loading path. For MAT 2, the Bauschinger effect is less than in MAT 1 for the C5T10R loading path, though the permanent softening is also pronounced. However, the C5T10CR loading condition showed a very noticeable stress overshoot upon changing the loading path and subsequent softening. Note that MAT 1 and MAT 2 were virtually designed to represent the characteristics of the e-HAH model.

In the following, the accuracy of each investigated stress integration algorithm is assessed with different time increments. The effect of the time-step size on the accuracy of the stress update is evaluated via one-element simulations under the two loading paths. Three different strain increments during tensile loading were considered, i.e., $5.0\times {10}^{-3}$, $5.0\times {10}^{-4}$, and $5.0\times {10}^{-6}$.

From the simulations, the evolution of the flow stress and r-value at the second loading step were evaluated, and the averaged relative errors are reported, using the following equations:

$$\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}{\mathrm{r}}_{\mathsf{\sigma }}=\frac{1}{\mathrm{N}}{\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{N}}\frac{\left|{\mathsf{\sigma }}_{\mathrm{m}}^{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}-{\mathsf{\sigma }}_{\mathrm{m}}\right|}{{\mathsf{\sigma }}_{\mathrm{m}}^{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}}\times 100(\%)}\hspace{0.17em}\mathrm{a}\mathrm{n}\mathrm{d}\hspace{0.17em}\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}{\mathrm{r}}_{\mathrm{r}}=\frac{1}{\mathrm{N}}{\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{N}}\frac{\left|{\mathrm{r}}_{\mathsf{\theta },\mathrm{m}}^{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}-{\mathrm{r}}_{\mathsf{\theta },\mathrm{m}}\right|}{{\mathrm{r}}_{\mathsf{\theta },\mathrm{m}}^{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}}\times 100(\%)},$$

(60)

where $\mathrm{N}$ is the total number of data points used for the error estimations.

4.1.1. C5T10R Loading Condition

Figure 6 and Figure 7 show the flow stress curves and r-value evolutions of MAT 1 and MAT 2, respectively, under the C5T10R loading path. For this analysis, three different strain increments were applied to each stress integration algorithm. The detailed values, including relative errors and CPU times, are also listed in

Table 3. Average relative errors in the flow curves and r-values of MAT 1 under C5T10R.

Algorithm	$\Delta {\mathsf{\epsilon }}_{\mathbf{x}\mathbf{x}}$	$\mathbf{E}\mathbf{r}\mathbf{r}\mathbf{o}{\mathbf{r}}_{\mathsf{\sigma }}(\%)$	$\mathbf{E}\mathbf{r}\mathbf{r}\mathbf{o}{\mathbf{r}}_{\mathbf{r}}(\%)$	Total CPU Time (Sec)
CPM	$5\times {10}^{-3}$	2.38	316	1.1
	$5\times {10}^{-4}$	0.42	5.83	5.7
	$5\times {10}^{-6}$	0.14	0.048	205
EBM-4R	$5\times {10}^{-3}$	4.34	5.23	0.5
	$5\times {10}^{-4}$	0.95	1.80	5.7
	$5\times {10}^{-6}$	0.14	0.035	232
EBM-14R/AD	$5\times {10}^{-3}$	1.62	0.63	0.6
	$5\times {10}^{-4}$	0.32	0.10	5.8
	$5\times {10}^{-6}$	0.14	0.010	231
EBM-14R/ND	$5\times {10}^{-3}$	1.62	0.63	0.7
	$5\times {10}^{-4}$	0.32	0.10	5.9
	$5\times {10}^{-6}$	0.14	0.010	249

and

Table 4. Average relative errors in the flow curves and r-values of MAT 2 under C5T10R.

Algorithm	$\Delta {\mathsf{\epsilon }}_{\mathbf{x}\mathbf{x}}$	$\mathbf{E}\mathbf{r}\mathbf{r}\mathbf{o}{\mathbf{r}}_{\mathsf{\sigma }}(\%)$	$\mathbf{E}\mathbf{r}\mathbf{r}\mathbf{o}{\mathbf{r}}_{\mathbf{r}}(\%)$	Total CPU Time (Sec)
CPM	$5\times {10}^{-3}$	0.42	124	0.8
	$5\times {10}^{-4}$	0.18	4.88	5.9
	$5\times {10}^{-6}$	0.21	0.048	206
EBM-4R	$5\times {10}^{-3}$	0.27	3.43	0.7
	$5\times {10}^{-4}$	0.17	0.19	6.0
	$5\times {10}^{-6}$	0.21	0.022	202
EBM-14R/AD	$5\times {10}^{-3}$	0.49	0.30	0.6
	$5\times {10}^{-4}$	0.25	0.05	5.7
	$5\times {10}^{-6}$	0.22	0.021	201
EBM-14R/ND	$5\times {10}^{-3}$	0.49	0.29	0.7
	$5\times {10}^{-4}$	0.25	0.04	6.0
	$5\times {10}^{-6}$	0.22	0.021	227

for MAT 1 and MAT 2, respectively. At a first glance, all the numerical algorithms investigated seem to yield rather similar results. However, some distinctive features can be observed in the predicted flow curves and r-values. Regarding the method for the yield surface gradient calculation, both approaches based on analytical derivatives (EBM-14R/AD) and finite difference (EBM-14R/ND) predicted almost the same level of accuracy when being implemented into the fully implicit algorithm with 14 residuals. Small differences exist between the two cases in the computational cost measured by the CPU time. Indeed, the EBM-14R/ND scheme using the finite difference method takes ~8% longer CPU time than the EBM-14R/AD scheme using the analytical derivatives.

The effect of the investigated algorithms on the accuracy and stability of the simulations with the e-HAH model grows considerably when the evolution of the yield surface becomes drastic within a given time increment. The flow stress curves predicted by the CPM and EBM-4R schemes show significant oscillations at the early strain range of MAT 1, when the strain increment is large, with $\Delta {\mathsf{\epsilon }}_{\mathrm{x}\mathrm{x}}=5\times {10}^{-3}$, as shown in Figure 6a. In contrast, the two algorithms with the same strain increment are stable for MAT 2, which exhibits a lower evolution rate than MAT 1 (Figure 7a). The state variable ${\mathrm{g}}_1$ controls the transient rate, and Figure 8a shows its oscillations for MAT 1 when the CPM and EBM-4R schemes are used. As expected, the ${\mathrm{g}}_1$ value of MAT 2 exhibits a stable evolution, even for the CPM and EBM-4R algorithms. In terms of the r-value evolution, the values predicted by the CPM approach show a much lower accuracy for a strain increment of $\Delta {\mathsf{\epsilon }}_{\mathrm{x}\mathrm{x}}=5\times {10}^{-3}$. This inaccuracy is attributed to the lack of knowledge of the flow rule as a residual during the algorithmic treatment in the CPM scheme.

4.1.2. C5T10CR Loading Condition

Figure 9 and Figure 10 show the flow stress and r-value evolutions under the C5T10CR cross-loading path for MAT 1 and MAT 2, respectively. The accuracy and stability of the investigated stress-integration algorithms for MAT 1 and MAT 2 are summarized in

Table 5. Average relative errors in the flow curves and r-values of MAT 1 under C5T10CR.

Algorithm	$\Delta {\mathsf{\epsilon }}_{\mathbf{x}\mathbf{x}}$	$\mathbf{E}\mathbf{r}\mathbf{r}\mathbf{o}{\mathbf{r}}_{\mathsf{\sigma }}(\%)$	$\mathbf{E}\mathbf{r}\mathbf{r}\mathbf{o}{\mathbf{r}}_{\mathbf{r}}(\%)$	Total CPU Time (Sec)
CPM	$5\times {10}^{-3}$	22.2	155	1.6
	$5\times {10}^{-4}$	0.56	19.4	5.9
	$5\times {10}^{-6}$	0.65	0.57	196
EBM-4R	$5\times {10}^{-3}$	22.4	43.8	0.5
	$5\times {10}^{-4}$	0.61	2.86	6.2
	$5\times {10}^{-6}$	0.65	1.60	210
EBM-14R/AD	$5\times {10}^{-3}$	0.38	2.78	0.8
	$5\times {10}^{-4}$	0.88	1.89	6.0
	$5\times {10}^{-6}$	0.66	1.57	222
EBM-14R/ND	$5\times {10}^{-3}$	0.44	4.77	1.1
	$5\times {10}^{-4}$	0.88	1.84	7.6
	$5\times {10}^{-6}$	0.66	1.58	260

and

Table 6. Average relative errors in the flow curves and r-values of MAT 2 under C5T10CR.

Algorithm	$\Delta {\mathsf{\epsilon }}_{\mathbf{x}\mathbf{x}}$	$\mathbf{E}\mathbf{r}\mathbf{r}\mathbf{o}{\mathbf{r}}_{\mathsf{\sigma }}(\%)$	$\mathbf{E}\mathbf{r}\mathbf{r}\mathbf{o}{\mathbf{r}}_{\mathbf{r}}(\%)$	Total CPU Time (Sec)
CPM	$5\times {10}^{-3}$	22.2	155	1.6
	$5\times {10}^{-4}$	0.56	19.4	5.9
	$5\times {10}^{-6}$	0.65	0.57	196
EBM-4R	$5\times {10}^{-3}$	22.4	43.8	0.5
	$5\times {10}^{-4}$	0.61	2.86	6.2
	$5\times {10}^{-6}$	0.65	1.60	210
EBM-14R/AD	$5\times {10}^{-3}$	0.38	2.78	0.8
	$5\times {10}^{-4}$	0.88	1.89	6.0
	$5\times {10}^{-6}$	0.66	1.57	222
EBM-14R/ND	$5\times {10}^{-3}$	0.44	4.77	1.1
	$5\times {10}^{-4}$	0.88	1.84	7.6
	$5\times {10}^{-6}$	0.66	1.58	260

, respectively. Similar to the C5T10R reverse loading path, the overall accuracy for the flow stress and r-values increases upon decreasing the strain increment as expected. However, due to the anisotropic hardening and its orthogonal distortion in the e-HAH model, distinct features can be also noticed, which mostly occur under the large strain increment $\Delta {\mathsf{\epsilon }}_{\mathrm{x}\mathrm{x}}=5\times {10}^{-3}$. The evolution of the r-value predicted by the CPM scheme shows an abnormal behavior as can be observed in Figure 9a. This is because the semi-explicit algorithm does not account for the residual minimization for a flow rule. Moreover, as shown in Figure 9a,b, the flow stress and r-value predicted by the CPM and EBM-4R schemes show a different trend with respect to the algorithms using 14 residual—based EBMs. The evolution of the two state variables, ${\mathrm{g}}_1$ and ${\mathrm{g}}_2$, under a strain increment $\Delta {\mathsf{\epsilon }}_{\mathrm{x}\mathrm{x}}=5\times {10}^{-3}$ are presented in Figure 11a,b. Under this loading path, the transient behavior of CPM and EBM-4R is controlled by ${\mathrm{g}}_2$, as the loading path indicator $\cos \mathsf{\chi }$ is negative for both algorithms (Figure 11c). However, for the other algorithms that include the exact solution, $\cos \mathsf{\chi }$ is positive, which leads ${\mathrm{g}}_1$ to have a major effect on the transient behavior. This circumstance represents the dominant cause behind the abnormally predicted flow stress curve and r-value evolutions for the two algorithms. In other words, the undesirable evolution of the state variables results in significantly deviated stress updates in the e-HAH model, as shown in Figure 11d.

Similar to MAT 1, oscillating flow curves are also predicted by the CPM and EBM-4R approaches for MAT 2 when the time increment is not small enough. In this loading path, the fluctuating flow behavior is attributed to the state variable ${\mathrm{g}}_{\mathrm{L}}$ associated with latent hardening (Figure 12). Note that the MAT 2 flow curve is influenced by ${\mathrm{g}}_{\mathrm{L}}$ more significantly than for MAT 1, due to the slower evolution of ${\widehat{h}}^{\mathbf{s}}$.

Regarding the computational time, all investigated numerical algorithms present rather similar values, especially when compared with the much larger differences encountered in the accuracy values, particularly under large strain increments. However, it is important to note that under the C5T10CR loading path, the finite difference gradient—based algorithm (EBM-14R/ND) is characterized by a ~10–15% longer CPU time, compared to the analytical gradient algorithm (EBM-14R/AD).

4.2. Forming Problem: S-Rail Forming and Springback

In this section, the proposed stress integration algorithms are applied to the simulation of large-scale industrial part forming with purely implicit FE software. Though the implicit FE solver was previously employed for large-scale models, this was often limited to simple material constitutive laws, such as isotropic and kinematic hardening. Therefore, this is the first time that the distortional hardening model with cross-hardening or softening (that is, the e-HAH model) is applied to the simulation of industrial forming processes with a static implicit FE solver. By contrast, numerous previous studies have employed the combined dynamic explicit and implicit solvers for complex forming and elastic-driven springback simulations, respectively. In particular, this happens for the case in which the elastic-plasticity constitutive laws become more complex, such as in the present anisotropic hardening model (e-HAH).

In this study, the S-rail part forming and springback simulations were performed as a benchmark for the industrial forming process [45]. The simulations were carried out using the static implicit solver, i.e., the ABAQUS/Standard. From the simulation results, the springback, equilibrium iterations, time increment, and computational time during the forming step were comparatively analyzed using the different CPM, EBM-14R/AD, and EBM-14R/ND algorithms. Real automotive sheet metals, composed of dual-phase steel (DP780) and stainless steel (STS), were employed in the simulations. Note that DP780 exhibits a significant cross-contraction (or softening) behavior under loading path changes, due to the large Baushinger effect induced by its martensitic islands embedded in the ferrite matrix. In contrast, the STS sheet exhibits stress over-shooting or cross-hardening behavior under loading path changes. The material properties and their related model parameters are listed in

Table 7. Material properties and model parameters of DP780 [41] and STS [64] for e-HAH.

Elastic property and isotropic hardening law parameters
Material	Young’s modulus (GPa)			Poisson’s ratio			K (MPa)			e0	n
DP780	194.0			0.33			1295.0			0.0008	0.142
STS	223.9			0.33			882.6			0.0056	0.224
Yld2000-2d coefficients
Material	M	α1		α2		α3	α4	α5	α6	α7	α8
DP780	6	0.946		1.022		1.015	1.000	1.011	0.968	1.009	1.006
STS	2	1.0		1.0		1.0	1.0	1.0	1.0	1.0	1.0
e-HAH model parameters
Material	q	k	k1	k2	k3	k4	k5	L	kL	S	kS
DP780	2	120.0	150.0	78.5	0.25	1.0	0.0	1.0	0.0	0.8	75.0
STS	2	18.7	580.0	10.0	0.90	0.88	240	1.65	384.0	1.0	0.0

. The FE model set-up, including tool and sheet dimensions, is presented in Figure 13.

The thickness of the blank sheet was set to 1.2 mm, and a blank holding force of 200 kN was applied during the forming simulation. The friction coefficient between the blank and tools was set to 0.05. A total punch displacement of 37 mm was applied. As for the blank elements, 4264 four-node shell elements with reduced integration (S4R) were used. The simulation consisted of the holding, forming, and springback steps. Springback simulations were conducted by removing the tools after applying constraints on three specific nodes on the X’-Z’ plane to prevent rigid body motion. In other words, the center node was kept fixed, one node on the Z’ axis was constrained along the X’ and Y’ directions, and the last node on the X’ axis was constrained along the Z’ direction. The simulation time for each step was set to one (however, this does not have a physical meaning due to the static implicit FE algorithm. This is a relative measure for the strain increment control for different stress-integration algorithms). The reference, minimum, and maximum time increments for the forming step were set to $5.0\times {10}^{-3}$, $1.0\times {10}^{-5}$, and $1.0\times {10}^{-2}$, respectively.

The sections A, B, and C of the springback are depicted in Figure 13c. Section A and C are on the Z-plane, $\pm 115\hspace{0.17em}\mathrm{m}\mathrm{m}$ away from the center. Section B lies on the Z’-plane, which is rotated by 30° around the Y-axis.

The springback results for section A, B, and C are presented for DP780 and STS in Figure 14. Almost no differences can be observed in the springback results obtained using the different algorithms. This may be due to the automatic time stepping built into the ABAQUS routine and to the fact that the time increment is small enough to have accurate solutions.

The averaged equilibrium iteration number, averaged time increment ($\Delta {\mathrm{t}}_{\mathrm{a}\mathrm{v}\mathrm{g}}$), and relative wallclock time are listed in

Table 8. The averaged equilibrium iteration number, the averaged time increment and the relative wallclock time during the forming step.

Algorithm	Average Equilibrium Iteration Number		$\Delta {\mathbf{t}}_{\mathbf{a}\mathbf{v}\mathbf{g}}$		Relative Wallclock Time
	DP780	STS	DP780	STS	DP780	STS
CPM	13.42	13.71	$1.21\times$ 10−3	$5.62\times$ 10−4	1.0	1.0
EBM-14R/AD	9.47	12.72	$2.23\times$ 10−3	$1.10\times$ 10−3	0.59	0.57
EBM-14R/ND	9.90	12.53	$2.26\times$ 10−3	$1.12\times$ 10−3	0.84	0.54

. The relative wallclock time is normalized with the value calculated using the CPM approach. Pronounced differences in the time increment and calculation time obtained from the different algorithms can be observed. The averaged time increment for CPM is almost half of that of the EBM algorithms. It should be considered that the tangent modulus calculated with CPM is not consistent with the Newton-Raphson method in ABAQUS, which requires a rather smaller time step, compared to the EBM algorithms.

For both materials, the EBM-14R/AD and EBM-14R/ND approaches show a negligible difference in the averaged time increment and average equilibrium iteration number. In other words, the approximate computational cost for global equilibrium calculations in the FE software is similar for the two algorithms. It is thus presumed that the two algorithms have a very similar accuracy under a given time increment. For the STS case, the EBM-14R/AD and EBM-14R/ND methods show similar wallclock times, whereas these times are different for the DP780 case. Indeed, the EBM-14R/ND scheme is approximately 40% slower than the EBM-14R/AD scheme. The difference in convergence speed may be due to the different exponents of the non-quadratic yield function used for the DP780 and STS materials, which also determine the nonlinearity of the yield function evolution.

5. Summary

In this study, a fully implicit stress integration algorithm was developed for the e-HAH anisotropic hardening model. The proposed algorithm is capable of reproducing both cross-hardening and softening under complex loading path changes by introducing the distortional hardening concept. Particularly, the proposed algorithm solves for the complete set of residuals defined from the e-HAH model in the context of the EBM. The major differences between the developed model and previously proposed stress integration algorithms can be summarized as follows.

The developed EBM algorithm is formulated based on a total of 14 residuals for the stress tensor, microstructure deviator tensor, and the whole state variables associated with the e-HAH model. On the contrary, the previous algorithms for the HAH or e-HAH model introduced a partially implicit scheme by considering only a limited number of residuals for greater simplicity. This led to unstable and inaccurate evolutions of the plastic state variables related to the distortions of the yield function under loading path changes.

For calculating the first and second derivatives of the yield surface, which are inevitably required in the formulation of the common predictor-corrector numerical schemes, the present model provides both analytical (EBM-AD) and finite difference—based numerical methods (EBM-ND). The analytical expressions of the e-HAH yield surface are given in the appendix.

The accuracy and robustness of the developed algorithm and its implementation were validated by a one-element analysis and application to a large-size industrial problem. For the one-element analysis, the compression-tension test and compression-cross-tension test were conducted for two representative materials. For the large-size problem, the S-rail benchmark forming and springback simulations were conducted. These validations can be summarized as follows.

The existing semi-explicit CPM and implicit EBM schemes with partially introduced residuals (EBM-4R) resulted in an abnormal evolution of the stress-strain curve or r-value for both the reversed and cross-loading conditions when the strain increment became larger. In contrast, the proposed fully implicit algorithm with the complete set of residuals (EBM-14R) showed stable and accurate results regardless of the investigated strain increments. Moreover, both methods based on analytical and numerical derivatives provided virtually the same accuracy for the EBM-14R algorithm.

For the cross-loading path, the lack of knowledge of the residuals for the evolution of the e-HAH yield surface led to the abnormally wrong rotation of the microstructure deviator ${\widehat{h}}^{\mathbf{s}}$, which caused undesirable evolutions of the other state variables. This justifies the better accuracy and robustness of the present fully implicit algorithm based on full considerations of the residuals of the e-HAH model.

For the S-rail forming and springback, the investigated algorithms were all successful without divergence, even with the static implicit solver and complicated e-HAH model. However, the EBM—based algorithms required less computational time than the CPM—based one. This is due to the smaller time step determined by the automatic time increment, which resulted from the non-consistent tangent modulus of the CPM algorithm.

Noticeable differences in the averaged time increments and equilibrium iteration number were observed between the EBM-AD and EBM-ND approaches in the S-rail forming simulations. However, the yield surface exponent, which determines the sharpness of the non-quadratic yield function at the bi-axial stress state, played a more dominant role in determining the computational cost in the numerical, derivative—based EBM algorithm.