## 4 The Distributed Order Damage Creep Model of Rock Mass
### Establishment and Properties of the Distributed Order Damper Damage Model  
(1) Establishment of the distributed order damper damage model  
According to the fractional order theory, the fractional order damper damage model was constructed to show accelerated rock deformation, and the constitutive relation could be expressed by Equation (40). \[\sigma(t)=E\eta^{\gamma}e^{-\lambda t}0_{t}^{\gamma}[\varepsilon(t)]\ \ \ \ \ \ (0\leq\gamma\leq 1) \tag{40}\]  
where \(E\eta^{\gamma}e^{-\lambda t}\) denotes the time-dependent viscosity coefficient and \(\lambda\) denotes the damage factor. By the Laplace transform, Equation (41) can be expressed:  
\[\varepsilon(s+\lambda)=\frac{\sigma(s)}{E\eta^{\gamma}(s+\lambda)^{\gamma}}, \ 0\leq\gamma\leq 1 \tag{41}\]  
To show the anisotropy of the rock mass, multiple fractional order damage dampers with various orders were connected in series, as shown in Figure 7.  
According to the calculation methods of the component series, Equation (42) could be obtained. \[\varepsilon(s+\lambda)=\sum_{i=1}^{n}\varepsilon_{i}(s+\lambda)=\sum_{i=1}^{n }\frac{\sigma(s)}{E\eta^{\gamma_{i}}(s+\lambda)^{\gamma_{i}}},0\leq\gamma\leq 1 \tag{42}\]  
With the same method and the assumption that \(\Delta\gamma\) was small enough, for \(0\leq a\leq\gamma_{i}\leq b\leq 1\), \([a,b]\) could be converted into \(N=\frac{b-a}{\Delta\gamma}\), and Equation (42) could be transformed into Equation (43). \[\varepsilon(s+\lambda)=\sum_{i=1}^{n}\frac{\sigma(s)}{E\eta^{(i\Delta\gamma)}( s+\lambda)^{(i\Delta\gamma)}} \tag{43}\]  
According to the definition of the definite integral, Equation (44) could be obtained. \[\varepsilon(s+\lambda)=\int_{a}^{b}\frac{\omega(\gamma)}{(s+\lambda)^{\gamma} }d\gamma\sigma(s) \tag{44}\]  
Figure 7: The distributed order damage damper. (a) Standard form of the distributed order damage damper, (b) The simplified form. Specifically, when \(\omega(\gamma)\) is the impulse function, Equation (44) could be changed into Equation (45). \[\varepsilon(s+\lambda)=\int_{a}^{b}\frac{\delta(\gamma-\gamma_{0})}{(s+\lambda)^ {\gamma}}d\gamma= \frac{1}{(s+\lambda)^{\gamma_{0}}}\sigma(s) \tag{45}\]  
By the inverse Laplace transform, Equation (45) could be converted into Equation (46). \[\sigma(t)=e^{-\lambda t}{}_{0}\mathbb{D}_{t}^{\gamma_{0}}[\varepsilon(t)] \tag{46}\]  
This was the fractional order damage damper. By the inverse Laplace transform of Equation (44), Equation (47) could be obtained. \[\varepsilon(t)=e^{\lambda t}L^{-1}\biggl{\{}\int_{a}^{b}\frac{\omega(\gamma)} {(s+\lambda)^{\gamma}}d\gamma\sigma(s)\biggr{\}},\;0\leq\gamma\leq 1 \tag{47}\]  
With N harmonics, the Fourier series approximation of \(\omega(\gamma)\) in Equation (47) was expressed by Equation (48). \[\varepsilon(s+\lambda)=\frac{\sigma(s)}{E\eta^{\gamma}(s+\lambda)^{\gamma}}, \ 0\leq\gamma\leq 1 \tag{41}\]0  
where  
\[\left\{ \begin{array}{ll}a_{n}=\frac{2}{b-a}\int_{a}^{b}\omega(\gamma)\cos( \tau_{n}\gamma)d\gamma\\ b_{n}=\frac{2}{b-a}\int_{a}^{b}\omega(\gamma)\sin(\tau_{n}\gamma)d\gamma\\ \tau_{n}=\frac{2\pi\gamma}{b-a}\end{array} \right.\;\;n=0,1,2,\cdots,N \tag{49}\]  
By further calculation, Equation (50) could be obtained. \[ \begin{array}{ll}\varepsilon(t)&=e^{\lambda t}L^{-1}\biggl{\{}\int_{a}^{b} \frac{\omega(\gamma)}{(s+\lambda)^{\gamma}}d\gamma\sigma(s)\biggr{\}}\\ &=e^{\lambda t}L^{-1}\biggl{\{}\int_{a}^{b}\left\{\frac{a_{0}}{2}+\sum\limits_ {n=1}^{N}\left[a_{n}\cos(\tau_{n}\gamma)+b_{n}\sin(\tau_{n}\gamma)\right] \right\}\frac{1}{(s+\lambda)^{\gamma}}d\gamma\sigma(s)\biggr{\}}\end{array}  \tag{50}\]  
Equation (50) could be solved by the inverse Laplace transform of \(\int_{a}^{b}\frac{\omega(\gamma)}{(s+\lambda)^{\gamma}}d\gamma\).