## 3 Stress Integration Algorithms of the e-HAH Model for Elastic-Plastic FE Simulations
### Stress Integration Algorithms for the e-HAH Model: Proposed Models
#### 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}^{(k)}_{1,n+1}=\Phi\Big{(}\sigma^{(k)}_{n+1},\,\mathrm{g}^{(k)}_{i,n+1}, \,\hat{\mathbf{h}}^{\mathbf{s}(k)}_{i,n+1}\Big{)}-\overline{\sigma}\Big{(} \overline{\mathrm{z}}^{(k)}_{n+1}\Big{)}\,,\,(i=1-7,\,5\,=\,\mathrm{L},\,6= \mathrm{S},\,7=\mathrm{R}), \tag{44}\]  
\[\mathrm{R}^{(k)}_{2,n+1}=-\Delta\varepsilon^{\mathrm{p}(k)}_{n+1}+\Delta \gamma\,\partial_{\sigma}\Phi^{(k)}_{n+1}, \tag{45}\]  
\[\mathrm{R}^{(k)}_{2+i,n+1}=\delta^{(k)}_{i,n+1}-\delta_{i,n}-\left(\frac{ \mathrm{dg}_{i,n+1}}{\mathrm{d}\overline{\varepsilon}}\right)^{(k)}\cdot \Delta\overline{\varepsilon}^{(k)}_{n+1}\,, \tag{46}\]  
\[\mathrm{R}^{(k)}_{10,n+1}=\Delta\hat{\mathbf{h}}^{\mathbf{s}(k)}_{n+1}-\Delta \hat{\mathbf{h}}^{\mathbf{s}}_{n}-\left(\frac{\mathrm{d}\hat{\mathbf{h}}^{ \mathbf{s}}_{n+1}}{\mathrm{d}\overline{\varepsilon}}\right)^{(k)}\cdot \Delta\overline{\varepsilon}^{(k)}_{n+1}, \tag{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, \(\hat{\mathbf{h}}^{\mathbf{s}}\). Linearization of Equations (44)-(47) leads to the following equations:  
\[ \begin{array}{l}\mathrm{R}^{(k)}_{1,n+1}-\partial_{\overline{\varepsilon}} \overline{\sigma}\Big{(}\overline{\varepsilon}^{(k)}_{n+1}\Big{)}\cdot \delta(\Delta\overline{\varepsilon}_{n+1})^{(k+1)}+\partial_{\sigma}\Phi^{(k)} _{n+1}:\delta(\Delta\sigma_{n+1})^{(k+1)}\\ +\sum\limits_{\hat{\mathbf{i}}}\partial_{\mathrm{g}_{i}}^{2}\Phi^{(k)}_{n+1} \cdot\delta\Big{(}\Delta\mathrm{g}_{i,n+1}\Big{)}^{(k+1)}+\partial_{\hat{ \mathbf{i}}^{\mathbf{s}}}\Phi^{(k)}_{n+1}:\delta\Big{(}\Delta\hat{\mathbf{h} }^{\mathbf{s}}_{n+1}\Big{)}^{(k+1)}=0\,,\\ \end{array}  \tag{48}\]  
\[ \begin{array}{l}\mathrm{R}^{(k)}_{2,n+1}+\partial_{\sigma}\[PERSON]^{(k)}_{n+1} \cdot\delta(\Delta\overline{\varepsilon}_{n+1})^{(k+1)}+\[PERSON]^{-1}_{n+1}:\delta (\Delta\sigma_{n+1})^{(k+1)}\\ +\sum\limits_{\hat{\mathbf{g}}}\partial_{\sigma_{\hat{\mathbf{g}}_{i}}}^{2} \[PERSON]^{(k)}_{n+1}\cdot\delta\Big{(}\Delta\mathrm{g}_{i,n+1}\Big{)}^{(k+1)} \cdot\Delta\overline{\varepsilon}^{(k)}_{n+1}+\partial_{\sigma\hat{\mathbf{h} }^{\mathbf{s}}}\Phi^{(k)}_{n+1}:\delta\Big{(}\Delta\hat{\mathbf{h}}^{\mathbf{ s}}_{n+1}\Big{)}^{(k)}\cdot\Delta\overline{\varepsilon}^{(k+1)}_{n+1}=0\\ \Xi^{-1}_{n+1}=\left(\mathbf{C}^{-1}+\Delta\overline{\varepsilon}^{(k)}_{n+1} \,\partial_{\sigma\sigma}^{2}\Phi^{(k)}_{n+1}\right)^{-1}\\ \end{array}  \tag{49}\]  
\[