## 3 Numerical Algorithm of Elastic-Plastic Damage Coupled Model  
The proposed elastic-plastic damage coupled model of rock was implemented into the software LS-DYNA as a user-defined material model [23], shown in Figure 1.  
In Figure 1, \(\overline{\sigma}_{ij}^{n}\) is effective stress tensor; \(\overline{\sigma}_{ij}^{n\,+1}\) is trial effective stress tensor; \(\sigma_{ij}^{n\,+1}\), \(\sigma_{ij}^{n\,+1}\) are stress tensor with damage; \(C_{cjkl}\) is elastic tensor; \(\Delta\varepsilon_{kl}\) is strain increment; \(D_{\text{td}}\) is tension damage degree rock micro-body (\(0<D_{\text{td}}<1\)).  
The coupling process between the Cowper-Symonds hardening model and the proposed damage model (shown in Figure 1) may be summarized as follows.  
_Step 1_ Calculate the effective stress tensor \(\overline{\sigma}_{ij}^{n\,+1}\) at \(t=t_{n}\) by Equation (22).  
\[\overline{\sigma}_{ij}^{n\,}=\frac{\sigma_{ij}^{n\,}}{1-D} \tag{22}\]  
_Step 2_ Calculate the trial effective stress tensor \(\overline{\sigma}_{ij}^{n\,+1}\) at \(t=t_{n+1}\) by Equation (23).  
\[\overline{\sigma}_{ij}^{n\,+1}=\overline{\sigma}_{ij}^{n\,}+C_{cjkl}\Delta \varepsilon_{kl} \tag{23}\]  
_Step 3_ Judge whether or not to yield by the given yield condition \(\phi\) with Equation (24).  
\[\phi=\overline{\sigma}_{i}^{*\,2}-\overline{\sigma}_{y}^{2} \tag{24}\]  
where \(\overline{\sigma}_{y}\) is current yield stress; \(\overline{\sigma}_{i}^{*}\) is trial effective stress.  
If \(\phi\leq 0\), the rock is still in the elastic deformation zone. At this time, the trial strain is equal to the real strain, \(\overline{S}_{ij}^{*\,n+1}=\overline{S}_{ij}^{n+1}\). Here, \(\overline{S}_{ij}^{n+1}\) is deviatoric stress tensor; \(\overline{S}_{ij}^{n\,+1}\) is trial deviatoric stress tensor, and we may go to Step 2 and calculate \(\overline{\sigma}_{ij}^{n\,+1}\).  
If \(\phi>0\), the material is in plastic state, and plastic strain should be modified by Equation (25).  
\[\overline{\varepsilon}_{\text{p}}^{\text{eff}^{n+1}}=\overline{\varepsilon}_ {\text{p}}^{\text{eff}^{n}}+\Delta\varepsilon_{\text{p}}^{\text{eff}}= \overline{\varepsilon}_{\text{p}}^{\text{eff}^{n}}+\frac{\overline{\sigma}_ {i}-\overline{\sigma}_{y}^{n}}{3G+E_{\text{p}}} \tag{25}\]  
where \(\overline{\varepsilon}_{\text{p}}^{\text{eff}^{n}}\), \(\overline{\varepsilon}_{\text{p}}^{\text{eff}^{n+1}}\) are effective plastic deformation; \(\Delta\varepsilon_{\text{p}}^{\text{eff}}\) is effective plastic strain increment; \(\overline{\sigma}_{y}^{n}\) is current yield stress; \(G\) is the shear modulus.  
Figure 1: Coupling between Cowper–Symonds hardening model and proposed damage model.  
Meanwhile, proportionally reduce the deviator stress using Equation (26) and make it return to the yield surface.  
\[\overline{S}_{ij}^{n+1}=\overline{S}_{ij}^{*\,n+1}-\frac{3G\Delta \overline{\sigma}_{\text{P}}^{\text{eff}}}{\overline{\sigma}_{i}}\overline{S}_{ ij}^{*\,n+1} \tag{26}\]  
_Step 4_ Calculate the effective stress \(\overline{\sigma}_{ij}^{n+1}\) at \(t=t_{n+1}\) by Equation (27).  
\[\overline{\sigma}_{ij}^{n+1}=\overline{S}_{ij}^{n+1}+(p^{n+1}+q) \delta_{ij} \tag{27}\]  
where \(p^{n+1}=-k\ln(\frac{V_{n+1}}{V_{0}})\), \(k\) is the bulk modulus of elasticity, \(V_{n+1}\) is the volume of element at \(t=t_{n+1}\), \(V_{0}\) is the volume of element at the original time; \(q\) is artificial volume viscous damping; \(\delta_{ij}\) is Kronecker delta symbol.  
_Step 5_ Calculate the current damage weightings (\(\alpha_{\text{t}}\) and \(\alpha_{\text{c}}\)) by Equation (21), and the current damage degree \(D\) by Equation (20). Comparing the current damage degree with the previous damage degrees, the maximum damage degree may be obtained till the current step.  
_Step 6_ Calculate the nominal stress \(\sigma_{ij}^{n+1}\) at \(t=t_{n+1}\) by Equation (28).  
\[\sigma_{ij}^{n+1}=(1-D)\overline{\sigma}_{ij}^{n+1} \tag{28}\]