## 2 Numerical Models  
In the Euclidean space, the position vector \(\mathbf{x}\) and its time derivative define the velocity as follows:  
\[\mathbf{v}=\dot{\mathbf{x}}. \tag{1}\]  
The velocity gradient \(\mathbf{L}\), rate of deformation \(\mathbf{D}\), and rate of spin \(\mathbf{W}\) are then defined as  
\[\mathbf{L}=grad\mathbf{v}. \tag{2}\]\[\mathbf{L}=\mathbf{D}+\mathbf{W},\;\mathbf{D}=\frac{1}{2}\Big{(}\mathbf{L}+\mathbf{L }^{\mathsf{T}}\Big{)},\;\text{and}\;\mathbf{W}=\frac{1}{2}\Big{(}\mathbf{L}- \mathbf{L}^{\mathsf{T}}\Big{)}. \tag{3}\]  
They can be separated into elastic and plastic parts as  
\[\mathbf{L}=\mathbf{L}^{\epsilon}+\mathbf{L}^{p},\;\mathbf{D}=\mathbf{D}^{ \epsilon}+\mathbf{D}^{p},\;\text{and}\;\mathbf{W}=\mathbf{W}^{\epsilon}+ \mathbf{W}^{p}. \tag{4}\]  
The superscripts '\(e\)' and '\(p\)' denote the elastic and plastic parts, respectively. The deformation gradient tensor \(\mathbf{F}\) is also defined as [33; 34]  
\[\mathbf{F}=\mathbf{R}\mathbf{U}=\mathbf{F}^{\epsilon}\mathbf{F}^{p}, \tag{5}\]  
\[\dot{\mathbf{F}}=\mathbf{L}\mathbf{F}, \tag{6}\]  
\(\mathbf{F}^{\epsilon}\) and \(\mathbf{F}^{p}\) represent the elastic and plastic parts of \(\mathbf{F}\), respectively. \(\mathbf{R}\) and \(\mathbf{U}\) are rotational and stretch tensors, respectively. The three deformation gradient tensors define each dilatancy as  
\[J=det(\mathbf{F}). \tag{7}\]  
\[J^{\epsilon}=det(\mathbf{F}^{\epsilon}). \tag{8}\]  
\[J^{p}=det(\mathbf{F}^{p})=1. \tag{9}\]  
Note that the elastic dilatancy \(J^{\epsilon}\) is the same as the total dilatancy \(J\), because of the isochoric condition of plasticity, as shown in Equation (9). In addition, \(\mathbf{F}\) and \(\mathbf{F}^{\epsilon}\) define strain tensors as follows:  
\[\mathbf{E}=\frac{1}{2}\Big{(}\mathbf{F}^{\mathsf{T}}\mathbf{F}-\mathbf{I} \Big{)}, \tag{10}\]  
\[\mathbf{L}=grad\mathbf{v}. \tag{2}\]0  
\[\mathbf{L}=grad\mathbf{v}. \tag{2}\]1  
\(\mathbf{E}\) is the Green-Lagrangian strain tensor. \(\mathbf{E}^{\epsilon}\) and \(\mathbf{E}^{p}\) are the elastic and plastic parts in \(\mathbf{E}\). \(\mathbf{I}\) is the identity matrix. The Cauchy stress tensor \(\mathbf{T}\) is now defined in this study as  
\[\mathbf{L}=grad\mathbf{v}. \tag{2}\]2  
\[\mathbf{L}=grad\mathbf{v}. \tag{2}\]3  
\[\mathbf{L}=grad\mathbf{v}. \tag{2}\]4  
Next, to define the plastic flow, the dissipation inequality should be considered as [31; 32; 33]  
\[\mathbf{L}=grad\mathbf{v}. \tag{2}\]5  
The flow rule defines the rate of plastic strain [35],  
\[\mathbf{L}=grad\mathbf{v}. \tag{2}\]6  
where \(\overline{\sigma}\) is a plastic potential function, and \(\Gamma\) is the plastic multiplier. \(\frac{\partial\overline{\sigma}}{\partial\mathbf{T}}\) and \(\Gamma\) control the direction and magnitude of plastic flow, respectively. Based on \(J_{2}\) associated flow rule, this work defines the equivalent stress as  
\[\mathbf{L}=grad\mathbf{v}. \tag{2}\]7\[\mathbf{L}=grad\mathbf{v}. \tag{2}\]8  
\(S_{ij}\) represents the components of the deviatoric stress, given by  
\[\mathbf{L}=grad\mathbf{v}. \tag{2}\]9  
\(\mathbf{e}_{i}\) denotes the basis vectors of the Euclidean space. The yield condition is then specified as  
\[\mathbf{L}=\mathbf{D}+\mathbf{W},\;\mathbf{D}=\frac{1}{2}\Big{(}\mathbf{L}+\mathbf{L }^{\mathsf{T}}\Big{)},\;\text{and}\;\mathbf{W}=\frac{1}{2}\Big{(}\mathbf{L}- \mathbf{L}^{\mathsf{T}}\Big{)}. \tag{3}\]0  
\[\mathbf{L}=\mathbf{D}+\mathbf{W},\;\mathbf{D}=\frac{1}{2}\Big{(}\mathbf{L}+\mathbf{L }^{\mathsf{T}}\Big{)},\;\text{and}\;\mathbf{W}=\frac{1}{2}\Big{(}\mathbf{L}- \mathbf{L}^{\mathsf{T}}\Big{)}. \tag{3}\]1  
\(\sigma_{h}\) is the hardening stress. \(A\) and \(B\) are material constants, and \(n\) is the exponent value for hardening. Solving Equations (17), (21), and (22) simultaneously provides the value of \(\Gamma\). Since the numerical method to solve \(\Gamma\) has been explained in many works [36], this study does not repeat it.