## 4 Fractional Creep Model of Gangue-Roo Combination
### Fractional Creep Constitutive Model of Gangue-Roof Combination  
Analysis of the time-strain curve from stepwise loading of gangue-roof composites reveals that the creep process exhibits distinct instantaneous creep, decelerating creep, and steady-state creep stages, with no accelerated creep phase observed. As depicted in Figure 13, a fractional-order viscoelastic creep model is developed. The model consists of a Poynting-Thomson model in series with a Hookean element, where the integer-order dashpot in the Poynting-Thomson model is substituted by an Abelian fractional-order dashpot [34; 35; 36; 37]. The detailed model configuration is illustrated in Figure 20.  
The constitutive relationship between the components in the model is as follows:  
\[\left\{ \begin{array}{l}\sigma_{1}=E_{1}\varepsilon_{1},\sigma_{2}=E_{2} \varepsilon_{2}(t),\sigma_{3}=E_{3}\varepsilon_{3}(t),\sigma_{4}=\eta^{\alpha}D ^{\alpha}[\varepsilon_{4}(t)]\\ \sigma=\sigma_{1}=\sigma_{2}+\sigma_{3},\sigma_{3}=\sigma_{4},\varepsilon(t)= \varepsilon_{1}+\varepsilon_{2}(t),\varepsilon_{2}(t)=\varepsilon_{3}+ \varepsilon_{4}(t)\end{array} \right. \tag{1}\]  
The Laplace transform of Equation (1) is carried out to obtain the following:  
\[\left\{ \begin{array}{l}L[\varepsilon(s)]=\frac{L[\sigma_{1}(s)]}{E_{1}}+ \frac{L[\sigma_{2}(s)]}{E_{2}}\\ L[\varepsilon(s)]=\frac{L[\sigma_{1}(s)]}{E_{1}}+\frac{L[\sigma_{2}(s)]}{E_{3}} +\frac{L[\sigma_{3}(s)]}{\eta^{\alpha}s^{\beta}}\end{array} \right. \tag{2}\]  
Simplify Equation (2) to the following:  
\[\frac{E_{2}+E_{3}}{E_{3}}L[\varepsilon(s)]+\frac{E_{2}}{\eta^{\alpha}s^{\alpha }}L[\varepsilon(s)]=\frac{E_{1}+E_{2}+E_{3}}{E_{1}E_{3}}L[\sigma(s)]+\frac{E_ {1}+E_{2}}{E_{1}(\eta^{\alpha}s^{\alpha}}L[\sigma(s)] \tag{3}\]  
The constitutive relation of the model can be obtained by Laplace inverse transformation of Equation (3), written as follows:  
\[\frac{E_{1}+E_{2}+E_{3}}{E_{1}E_{3}}\eta^{\alpha}D^{\alpha}[\sigma(t)]+\frac{ E_{1}+E_{2}}{E_{1}}\sigma(t)=\frac{E_{2}+E_{3}}{E_{3}}\eta^{\alpha}D^{\alpha}[ \varepsilon(t)]+E_{2}\varepsilon(t) \tag{4}\]  
When the stress \(\sigma(t)\) is a constant \(\sigma_{c}\), Laplace transform is performed on Equation (4) to obtain the following:  
\[\frac{E_{1}+E_{2}+E_{3}}{E_{1}E_{3}}\frac{\eta^{\alpha}\sigma_{c}}{s^{1- \alpha}}+\frac{E_{1}+E_{2}}{E_{1}}\frac{\sigma_{c}}{s}=L[\varepsilon(t)]\left( \frac{E_{2}+E_{3}}{E_{3}}\eta^{\alpha}s^{\alpha}+E_{2}\right) \tag{5}\]  
Change Equation (5) into a two-parameter Mittag-Leffler function expression, written as follows:  
\[L[\varepsilon(t)]=\frac{\sigma_{c}\cdot\frac{E_{1}+E_{2}+E_{3}}{E_{1}(E_{2}+E _{3})}\cdot s^{\alpha-1}}{s^{\alpha}+\frac{E_{2}E_{3}}{(E_{2}+E_{3})\eta^{ \alpha}}}+\frac{\sigma_{c}\cdot\frac{(E_{1}+E_{2})E_{3}}{E_{1}(E_{2}+E_{3}) \eta^{\alpha}}\cdot s^{-1}}{s^{\alpha}+\frac{E_{2}E_{3}}{(E_{2}+E_{3})\eta^{ \alpha}}} \tag{6}\]  
Figure 20: Fractional creep constitutive model of gangue–roof combination. The inverse Laplace transform of the two-parameter Mittag-Leffler function is performed on Equation (6), and the creep relationship of the model can be obtained as follows:  
\[\epsilon(t)=\sigma_{c}\bigg{\{}\frac{E_{1}+E_{2}+E_{3}}{E_{1}(E_{2}+E_{3})}E_{ \alpha,1}\bigg{[}-\frac{E_{2}E_{3}}{(E_{2}+E_{3})\eta^{\alpha}}t^{\alpha}\bigg{]} +\frac{(E_{1}+E_{2})E_{3}}{E_{1}(E_{2}+E_{3})\eta^{\alpha}}t^{\alpha}E_{\alpha, A+1}\bigg{[}-\frac{E_{2}E_{3}}{(E_{2}+E_{3})\eta^{\alpha}}t^{\alpha}\bigg{]} \bigg{\}} \tag{7}\]  
The Mittag-Leffler function satisfies the following properties:  
\[\left\{ \begin{array}{l}E_{1,1}(z)=e^{z}\\ E_{\alpha,\beta}(z)=\frac{1}{\Gamma(\beta)}+zE_{\alpha,A+\beta}(z)\end{array} \right.