## 4 Fractional Nonlinear Damage Creep Model
### Establishment of the Fractional Nonlinear Creep Model  
The classical Burgers creep model, composed of a Kelvin model and a Maxwell model connected in series, can only describe the creep characteristics of some soft rocks such as clay, asphalt, and so on. However, it is unable to describe the accelerating creep stage. Based on previous research [23; 24; 33], we establish a modified Burgers model by replacing the Newton dashpot with an Abel dashpot. When the rock begins to enter the accelerating creep stage, a nonlinear dashpot is introduced and connected in parallel with a plastic switch element. Then, these two components are connected in series, as illustrated in Figure 14. This newly established fractional nonlinear model, abbreviated as FNC, is capable of appropriately depicting the entire creep process of rocks: (1) a Hooke spring to reflect transient creep; (2) a modified Kelvin model comprising a Hooke spring and an Abel dashpot to reflect decelerating creep; (3) an Abel dashpot to reflect steady-state creep; (4) the nonlinear dashpot describes the viscoplastic deformation when the switch is on. Based on the above analysis, the total strain can be expressed as follows:  
\[\varepsilon=\varepsilon_{e}+\varepsilon_{v}+\varepsilon_{kv}+\varepsilon_{n1} \tag{17}\]  
Figure 14: The fractional nonlinear creep model. Figure 13: Creep curves of Abel dashpot with different fractional orders. \(\varepsilon_{e}\) is the instantaneous elastic strain:  
\[\varepsilon_{e}=\frac{\sigma}{E_{1}} \tag{18}\]  
\(\varepsilon_{kv}\) is the strain of the decelerating creep stage, and the improved Kelvin model follows the following pattern during the creep process:  
\[\left\{ \begin{array}{l}\varepsilon_{kv}=\varepsilon_{H}=\varepsilon_{A}\\ \sigma=E_{2}\varepsilon_{kv}+\eta_{2}D^{n}(\varepsilon_{kv})\end{array} \right. \tag{19}\]  
where \(\varepsilon_{H}\) is the strain of the Hooke spring, \(\varepsilon_{A}\) is the strain of the Abel dashpot, \(E_{2}\) is the elastic modulus, and \(\eta_{2}\) is the viscosity coefficient of the Abel dashpot. Equation (19) can be expressed as:  
\[\frac{E_{2}}{\eta_{2}}\varepsilon_{kv}+D^{n}(\varepsilon_{kv})=\frac{\sigma }{\eta_{2}} \tag{20}\]  
Consider the initial condition of \(\varepsilon_{kv}\) (0) = 0 (\(t\) = 0) by means of the Laplace transform and inverse transformation:  
\[\varepsilon_{kv}=\frac{\sigma}{\eta_{2}}\sum_{k=0}^{\infty}\frac{(-(E_{2}/ \eta_{2}))^{k}t^{n(1+k)}}{n(1+k)\Gamma[(1+k)n]} \tag{21}\]  
[PERSON] et al. [21] used the Mittag-Leffler function to simplify Equation (21):  
\[\varepsilon_{kv}=\frac{\sigma}{E_{2}}-\frac{\sigma}{E_{2}}E_{n}\left(-\frac{E _{2}}{\eta_{2}}t^{n}\right) \tag{22}\]  
\(E_{n}\) represents the Mittag-Leffler function with the following expression:  
\[E_{n}(t)=\sum_{k=0}^{\infty}\frac{t^{k}}{\Gamma[1+nk]} \tag{23}\]  
\(\varepsilon_{v}\) is the steady creep strain:  
\[\varepsilon_{v}=\frac{\sigma}{\eta_{1}}\frac{t^{n}}{\Gamma(1+n)} \tag{24}\]  
where \(\eta_{2}\) is the viscosity coefficient of Abel dashpot. \(\varepsilon_{vp}\) is the deformation in the accelerated creep phase, triggered by the yield stress, \(\sigma_{s}\), which is the first level of stress before the accelerated creep stage. When we apply the Laplace transform to Equation (16), we achieve the following expression:  
\[\widetilde{\varepsilon}_{n1}(s)=\frac{\sigma}{\eta_{n1}s^{3}} \tag{25}\]  
where \(\widetilde{\varepsilon}_{n1}\) is the Laplace transform of \(\varepsilon\), and \(s\) is the complex variable of the Laplace transform. By applying the inverse Laplace transform to Equation (25), we obtain the following expression:  
\[\begin{cases}\varepsilon_{n1}=\frac{\sigma}{2\eta_{n1}}\  
\[\varepsilon(t)=\left\{