## 3 Stress Integration Algorithms of the e-HAH Model for Elastic-Plastic FE Simulations
### Stress Integration Algorithms for the e-HAH Model: Review of Existing Models  
#### 3.3.1 Cutting-Plane Method (CPM)  
[PERSON] 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{\varepsilon}\), is the only independent variable. For each iteration of the current time step \(\text{n}+1\), \(\delta\Delta\overline{\varepsilon}_{\text{n}+1}^{(\text{k}+1)}\) is calculated from Equations (16)-(18), using linearized Equation (20). \[\mathbf{R}_{\text{n}+1}^{(\text{k})}+\left(\frac{\partial\mathbf{R}}{\partial \Delta\overline{\varepsilon}}\right)_{\text{n}+1}^{(\text{k})}\delta\Big{(} \Delta\overline{\varepsilon}_{\text{n}+1}^{(\text{k}+1)}\Big{)}=0\,, \tag{21}\]\[\left(\frac{\partial\mathbf{R}}{\partial\Delta\bar{\varepsilon}}\right)_{ \mathrm{n+1}}^{\mathrm{(k)}}=\left(\frac{\partial\mathbf{\Phi}}{\partial\mathbf{ \sigma}}:\frac{\partial\mathbf{\sigma}}{\partial\Delta\bar{\varepsilon}}\right)_{ \mathrm{n+1}}^{\mathrm{(k)}}-\left(\frac{\partial\overline{\sigma}}{\partial \Delta\bar{\varepsilon}}\right)_{\mathrm{n+1}}^{\mathrm{(k)}}, \tag{22}\] \[\left(\frac{\partial\mathbf{\sigma}}{\partial\Delta\bar{\varepsilon}} \right)_{\mathrm{n+1}}^{\mathrm{(k)}}=-\mathbf{C}:\left(\frac{\partial\Phi}{ \partial\mathbf{\sigma}}\right)_{\mathrm{n+1}}^{\mathrm{(k)}},\] (23) \[\delta(\Delta\bar{\varepsilon})_{\mathrm{n+1}}^{\mathrm{(k+1)}}= \frac{\mathrm{R}_{\mathrm{n+1}}^{\mathrm{(k)}}}{\left(\frac{\partial\Phi}{ \partial\mathbf{\sigma}}:\mathbf{C}:\frac{\partial\Phi}{\partial\mathbf{\sigma}} \right)_{\mathrm{n+1}}^{\mathrm{(k)}}+\left(\frac{\partial\overline{\sigma}}{ \partial\Delta\bar{\varepsilon}}\right)_{\mathrm{n+1}}^{\mathrm{(k)}}}. \tag{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\bar{\varepsilon}_{\mathrm{n+1}}^{\mathrm{(k+1)}}=\Delta\bar{\varepsilon }_{\mathrm{n+1}}^{\mathrm{(k)}}+\delta\Big{(}\Delta\bar{\varepsilon}_{\mathrm{ n+1}}^{\mathrm{(k+1)}}\Big{)}, \tag{25}\]  
\[\mathrm{g}_{\mathrm{i,n+1}}^{\mathrm{(k+1)}}=\mathrm{g}_{\mathrm{i,n+1}}^{ \mathrm{(k)}}+\delta\Big{(}\mathrm{g}_{\mathrm{i,n+1}}\Big{)}^{\mathrm{(k)}}= \mathrm{g}_{\mathrm{i,n+1}}^{\mathrm{(k)}}+\left(\frac{\mathrm{d}\mathrm{g}_{ \mathrm{i,n+1}}}{\mathrm{d}\bar{\varepsilon}}\right)^{\mathrm{(k)}}\cdot \delta(\Delta\bar{\varepsilon}_{\mathrm{n+1}})^{\mathrm{(k+1)}},\,(\,\mathrm{ i}=1\,\sim 4,\,\mathrm{L},\,\mathrm{S},\,\mathrm{R}), \tag{26}\]  
\[\hat{\mathbf{h}}_{\mathrm{n+1}}^{\mathrm{s(k+1)}}=\hat{\mathbf{h}}_{\mathrm{n +1}}^{\mathrm{s(k)}}+\delta\Big{(}\hat{\mathbf{h}}_{\mathrm{n+1}}^{\mathrm{s (k)}}\Big{)}^{\mathrm{(k)}}=\hat{\mathbf{h}}_{\mathrm{n+1}}^{\mathrm{s(k)}}+ \left(\frac{\mathrm{d}\hat{\mathbf{h}}_{\mathrm{n+1}}^{\mathrm{s}}}{\mathrm{d }\bar{\varepsilon}}\right)^{\mathrm{(k)}}\cdot\delta(\Delta\bar{\varepsilon}_ {\mathrm{n+1}})^{\mathrm{(k+1)}}, \tag{27}\]  
\[\Delta\mathbf{\sigma}_{\mathrm{n+1}}^{\mathrm{(k+1)}}=\mathbf{C}:\Big{(}\Delta \bar{\varepsilon}_{\mathrm{n+1}}-\Delta\bar{\varepsilon}_{\mathrm{n+1}}^{ \mathrm{(k+1)}}\Big{)}, \tag{28}\]  
The iterative procedure continues until the following condition is satisfied:  
\[\Phi\Big{(}\mathbf{\sigma}_{\mathrm{n+1}}^{\mathrm{(k+1)}},\,\hat{\mathbf{g}}_{ \mathrm{i,n+1}}^{\mathrm{(k+1)}},\,\hat{\mathbf{h}}_{\mathrm{n+1}}^{\mathrm{s (k+1)}}\Big{)}-\overline{\sigma}\Big{(}\bar{\varepsilon}_{\mathrm{n+1}}^{ \mathrm{(k+1)}}\Big{)}<\overline{\sigma}(\bar{\varepsilon}_{\mathrm{n}})\cdot \mathrm{Tol},\,\left(\mathrm{Tol}=10^{-6}\right). \tag{29}\]  
Therefore, in the CPM approach, \(\overline{\varepsilon}\) is the only unknown variable, and all other variables are treated as a function of \(\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\Phi}{\partial\mathbf{\sigma}}=\Phi^{1-\mathrm{q}}\left[