## 4 Numerical Simulation of Surface Deformation Induced by Backfill Mining
### Simulation Methods of Gangue Backfill Materials
\tag{14}\]  
where \(S_{\mathrm{ij}}^{N}\) and \(S_{\mathrm{ij}}^{O}\) separately represent the current and previous deviatoric stress tensors in a time step; \(e_{\mathrm{ij}}^{K_{\mathrm{k}}N}\) and \(e_{\mathrm{ij}}^{K_{\mathrm{k}}O}\) separately denote the current and previous deviatoric strain tensors in a time step; \(K_{\mathrm{k}}\) is the \(K_{\mathrm{k}}\)th part of the element model connected in tandem. By solving Equations (12)-(14) simultaneously, the following relationships of current and previous deviatoric strain tensors in a time step in the fractional [PERSON] model can be obtained:  
\[\left\{ \begin{array}{l}\Delta\epsilon_{\mathrm{ij}}^{e}=(S_{\mathrm{\bar{ \mathrm{\eta}}}}^{\mathrm{e}N}-S_{\mathrm{\bar{\mathrm{\eta}}}}^{\mathrm{e}O}) /2G_{\mathrm{e}}\\ \Delta\epsilon_{\mathrm{ij}}^{\mathrm{e}N}-(S_{\mathrm{\bar{\mathrm{\eta}}}}^{ \mathrm{e}N}+S_{\mathrm{\bar{\mathrm{\eta}}}}^{\mathrm{e}O})\Delta t^{\alpha} /4\eta_{\mathrm{v}}^{\alpha}\\ \Delta\epsilon_{\mathrm{ij}}^{\mathrm{e}V}=((S_{\mathrm{\bar{\mathrm{\eta}}}}^{ \mathrm{e}N}+S_{\mathrm{\bar{\mathrm{\eta}}}}^{\mathrm{e}VO})\Delta t^{\beta}+( D/C-1)e_{\mathrm{\bar{\mathrm{\eta}}}}^{\mathrm{e}VO})/C\end{array} \right. \tag{15}\]  
\(C\) and \(D\) in Equation (15) have the following relationship:  
\[\left\{ \begin{array}{l}C=4\eta_{\mathrm{e}V}^{\beta}+2G_{\mathrm{e}V}\Delta t ^{\beta}\\ D=4\eta_{\mathrm{e}V}^{\beta}-2G_{\mathrm{e}V}\Delta t^{\beta}\end{array} \right. \tag{16}\]  
By substituting Equation (15) into Equation (8) and combining the result with Equation (14), the relationship between the current and previous deviatoric stress tensors in a time step in the fractional Burgers' model can be derived:  
\[S_{\mathrm{\bar{\mathrm{\eta}}}}^{N}=\frac{1}{E}(\Delta\epsilon_{\mathrm{\bar {\mathrm{\eta}}}}+FS_{\mathrm{\bar{\mathrm{\eta}}}}^{O}-(D/C-1)e_{\mathrm{\bar{ \mathrm{\eta}}}}^{\mathrm{e}VO}) \tag{17}\]  
\(E\) and \(F\) in Equation (17) satisfy the following relationship:  
\[\left\{ \begin{array}{l}E=\frac{\Delta t^{\alpha}}{C}+\frac{1}{2G_{\mathrm{v} }}+\frac{\Delta t^{\alpha}}{4\eta_{\mathrm{v}}^{\alpha}}\\ F=\frac{1}{2G_{\mathrm{v}}}-\frac{\Delta t^{\alpha}}{C}-+\frac{\Delta t^{\alpha} }{4\eta_{\mathrm{v}}^{\alpha}}\end{array} \right. \tag{18}\]  
The fractional Burgers' model is realized using object-oriented C++ programming. After compilation, the dynamic link library (DLL) files are generated and embedded in the numerical simulation software, FLAC\({}^{\mathrm{3D}}\), to be invoked by users. Based on the built-in Burgers constitutive model of FLAC\({}^{\mathrm{3D}}\), the fractional Burgers' model is then developed, including three parts: modifying header files (.h files), modifying program source files (.cpp files), and generating DLL file (.dl). In FLAC\({}^{\mathrm{3D}}\), the newly developed fractional Burgers' model is invoked. The whole calculation is to calculate the new deviatoric stress, deviatoric strain, and spherical stress based on the initial stress using Equations (15) and (17); then, the unbalanced force, node rate, and displacement are calculated based on the new stress and strain and whether the unbalanced force is converged or not is judged. If the unbalanced force has converged, the calculation ends; otherwise, it returns to continue the next round of calculation. The secondary development flowchart of the fractional Burgers' model in FLAC\({}^{\mathrm{3D}}\) is illustrated in Figure 7.  
To test correctness of the secondarily developed fractional Burgers' model, example verification was performed through numerical simulation under conditions identical to those used in the creep compression tests. The shape and dimensions of the model were the same as the specimens in the creep compression tests. Horizontal constraints were applied to as boundary conditions of four-side walls, the vertical constraint was applied to the bottom boundary, and axial stress was applied to the top boundary. The axial stresses applied in different stages were 5, 10, 15, and 20 MPa, respectively.