# Creep Model of Weakly Cemented Soft Rock Considering Damage and Secondary Development in FLAC3D
## 3 Creep Damage Constitutive Model Considering Damage for Weakly Cemented Soft Rock
### Three-Dimensional Extension of Improved Nishihara Model  
As surrounding rock consistently experiences three-dimensional stress conditions, the model above derived under two-dimensional conditions may fail to accurately simulate the long-term strength characteristics. A three-dimensional extension of the model is required to enhance the model's practicality and facilitate the subsequent derivation of the finite difference schemes of the constitutive model [47; 48; 49; 50]. For the three-dimensional stress state model, the expressions for stress and strain are as follows:  
\[\sigma_{ij}=\delta_{ij}\sigma_{m}+S_{ij} \tag{7}\]  
\[\epsilon_{ij}=\delta_{ij}\epsilon_{m}+\epsilon_{ij} \tag{8}\]  
where \(S_{ij}\) and \(\sigma_{m}\) represent the deviatoric and spherical stresses, respectively; \(\epsilon_{m}\) and \(e_{ij}\) denote the deviatoric and spherical strains, respectively; and \(\delta_{ij}\) is the Christoffel symbol. According to the generalized Hooke's law, the Hooke three-dimensional constitutive relationship is as follows:  
\[\epsilon_{ij}^{\rm e}=\frac{1}{2G}s_{ij}+\frac{1}{3K_{0}}\delta_{ij}\sigma_{ \rm m} \tag{9}\]  
\(G\) and \(K_{0}\) are the bulk modulus and shear modulus of the Hook body, respectively. The three-dimensional stress criterion stipulates that spherical stress only induces instantaneous deformation, while rock creep is derived from deviatoric stress. By introducing a damage factor, Equation (8) is modified to a Hookean three-dimensional constitutive relationship considering damage:  
\[s_{ij}=2G(1-D(t))\epsilon_{ij}^{\rm e} \tag{10}\]  
According to research by [PERSON] et al. [51], when volume change is assumed elastic and rheological properties mainly manifest in shear deformations, the Kelvin-Voigt three-dimensional constitutive relationship considering damage is formulated as follows:  
\[\epsilon_{ij}^{\rm ev}=\frac{s_{ij}}{2K[1-D(t)]}(1-{\rm e}^{-\frac{K}{\eta^{ \rm m}}t}) \tag{11}\]  
Both the Bingham model and the winding element are activated upon reaching plastic failure; thus, the Bingham three-dimensional constitutive relationship considering damage is expressed as follows:  
\[\epsilon_{ij}^{\rm ep}=\frac{<F>}{2\eta^{\rm ep}[1-D(t)]}\frac{\partial g}{ \partial\sigma_{ij}}t \tag{12}\]  
The three-dimensional constitutive relationship for the winding element considering damage is delineated as follows:\[\varepsilon_{ij}^{\delta}=\frac{\lambda<F>}{4\rho h^{2}[1-D(t)]}\frac{\partial g}{ \partial\sigma_{ij}}t^{3} \tag{13}\]  
where <F> is the switching function, and \(\left\langle F\right\rangle=\left\{ \begin{array}{l}0,F\leqslant 0\\ F,F>0\end{array} \right..F\) is the yield function, and \(g\) is the plastic potential function. Based on the above equations, the total strain of the model is computed as the aggregate of the strains from its components, and the three-dimensional form of the creep damage model considering damage for weakly cemented soft rock is as follows:  
\[\left\{ \begin{array}{l}\varepsilon_{ij}=\frac{1}{2C[1-D(t)]}s_{ij}+\frac{1} {3K_{0}[1-D(t)]}\delta_{ij}\sigma_{m}+\frac{s_{ij}}{2K[1-D(t)]}(1-\mathrm{e}^{ -\frac{K}{\theta^{\mathrm{w}}}t}),\sigma_{\mathrm{D}}<\sigma_{\mathrm{s}}\\ \varepsilon_{ij}=\frac{1}{2C[1-D(t)]}s_{ij}+\frac{1}{3K_{0}[1-D(t)]}\delta_{ ij}\sigma_{m}+\frac{s_{ij}}{2K[1-D(t)]}(1-\mathrm{e}^{-\frac{K}{\theta^{ \mathrm{w}}}t})+\frac{<F>}{2\eta^{\mathrm{w}}[1-D(t)]}\frac{\partial g}{ \partial\sigma_{ij}}t+\frac{\lambda<F>}{4\rho h^{2}[1-D(t)]}\frac{\partial g} {\partial\sigma_{ij}}t^{3},\sigma_{\mathrm{D}}\geq\sigma_{\mathrm{s}}\end{array} \right. \tag{14}\]  
In laboratory experiments, confining pressure is generally applied under equal confining conditions (\(\sigma_{2}=\sigma_{3}\)), and the three-dimensional creep equation of the model is expressed by the following:  
\[\left\{