## 3 Stress Integration Algorithms of the e-HAH Model for Elastic-Plastic FE Simulations  
### Motivation and General Statement  
The e-HAH model features a larger number of plastic state variables than the conventional isotropic hardening model. In the [PERSON] 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 (\(\mathrm{\overline{\varepsilon}}\) and \(\mathrm{g}_{\mathrm{i}}\) (\(i=1\sim 4\), L, S, R)) 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.