The Seepage-Creep Numerical Simulation Model of Coal Measures Sandstone Based on Particle Discrete Element

The stability of the surrounding rock in the goaf is one of the main factors to ensure the safety of the mining production. In this study, a coupled numerical simulation model of particle flow code, computational fluid dynamics, and parallel-bounded stress corrosion (PFC-CFD-PSC model) was proposed to conduct the seepage-creep coupling simulation through PFC 5.0 software. Based on the mechanical results, the crack propagation law of the coal measure sandstone under the mining-induced seepagecreep coupling load was analyzed, and the distribution of pressure gradient in the fluid grid was studied. The results show that seepage pressure has a significant effect on the damage of the rock mass; the greater the seepage pressure on the coal measure sandstone model, the more penetration cracks are generated and the greater the degree of macroscopic damage. When the residual strength is 70% of the peak strength, the damage degree of the coal measure sandstone increases with the increase of the seepage pressure; except for the four corners, the pressure gradient of the fluid at the top of the rock mass model is relatively large. The maximum value of the fluid grid pressure gradient increases with the increasing seepage pressure; it increases slowly under the seepage pressure of less than 1.5MPa and increases sharply under the seepage pressure of greater than 1.5MPa.


Introduction
Deep coal mining has been widely used in the mining areas of China [1][2][3]. However, the deep geological environment is extremely complex and characterized by high ground temperature, high ground stress, and high karst water pressure and strong mining disturbance [2,4,5]. As a result, the stability of the overburden in the goaf is easily affected during the mining. To control the overburden movement and ensure the safety of mining production, it is important to study the mechanical behavior of the rock mass under the complex condition at the mining goaf.
Due to the existence of the underground water, the fluid will influence the mechanical behavior of overburden. The failure and the permeability change of the rock mass are closely related to the evolution of microscopic damage and the generation of macroscopic cracks [6][7][8][9][10]. To effectively understand the relationship between rock failure and water seepage, the fluid-solid coupling problem has been investigated from stress state analysis to failure process analysis. The existing studies have shown that the seepage pressure is inversely proportional to the crack initiation strength, while the confining pressure is proportional to the crack initiation strength; the crack length and the curvature radius of the crack tip are also inversely proportional to the crack initiation strength [11][12][13]. Besides, a generalized particle dynamics (GPD) method was proposed to conduct theoretical and numerical research on the crack propagation of the fractured rock mass under the coupling action of difference stress [14,15].
As a common numerical calculation method for data visualization, the computational fluid dynamics (CFD) can be used to simulate and analyze physical phenomena such as liquid movement by the finite volume method [16].
Although the traditional CFD method can numerically simulate the local mesoscopic of rock mass media, the boundary of the solid particles needs to be described by a body-fitted grid, and then, the Navier-Stokes equation needs to be solved. For porous media composed of many solid particles, a large number of body-fitted grids is required in the CFD method, resulting in the time-consuming calculation of the CFD model and the high cost in engineering application. Some researchers simulated the complex flow field in porous media by using the lattice Boltzmann method [17][18][19][20]. Through the combination of the lattice Boltzmann method and the discrete element method (LB-DEM), Feng et al. realized the two-dimensional coupling calculation of solid particles and fluids [21]. In this method, a structured background grid was used to simulate the particle boundary, an interaction model between the fluid and the particle boundary was established, and the characteristics of the local flow field near the particles were described accurately. Although the calculation amount of the LB-DEM is an order of magnitude smaller than that of the traditional CFD method, the calculation efficiency of the LB-DEM is still unacceptable when the number of solid particles reaches more than 10,000. Besides, the particle flow code and computational fluid dynamics (PFC-CFD) coupling method is also an effective method to achieve fluid and solid coupling [22][23][24]. Specifically, the rock skeleton is built by the particle and the uniform fluid grids are used to simulate the fluid flowing through the rock skeleton; then, the interaction between the mechanical field and seepage field is simulated.
Since the deformation of deep surrounding rock is a long-term process [25][26][27], creep and seepage phenomena need to be considered synchronously when studying the mechanical behavior of deep rock mass [28]. Based on the theory of stress erosion, Potyondy proposed a parallelbonded stress corrosion (PSC) model by the parallel bond model of PFC software and investigated the creep and damage evolution of granite specimens by the proposed PSC model [29]. After that, research on the creep behavior of the rock based on the PSC model has been carried out and the reliability of this model was proved through the comparison of experimental results of different models [30,31].
To study seepage-creep coupling behavior, a particle flow code, computational fluid dynamics, and parallelbounded stress corrosion (PFC-CFD-PSC) numerical simulation method was proposed based on the particle discrete element and computational fluid dynamics in this study. Then, the proposed model was used to simulate and analyze the coal measure sandstone of the surrounding rock in the deep goaf. The crack propagation law of the coal measure sandstone and the seepage characteristics under the combined action of seepage field and creep field were discussed. This study provides theoretical guidance for the deformation control of the surrounding rock and the development of underground water storage engineering.

Methodology
In this study, a 3D PFC-CFD-PSC rock model was established by the PFC3D5.0 software, and the fracture development of coal-measure sandstone under the coupled seepage-creep load was simulated.

Seepage Simulation in the PFC-CFD-PSC Model
2.1.1. Forces on Particles. The particles in fluids are subject to several forces, such as gravity, resistance, buoyancy, pressure gradient force, false mass force, and Basset force. Previous studies have shown that the latter three (pressure gradient force, false mass force, and Basset force) have little effect on particle movement than the former three (gravity, resistance, buoyancy, and pressure gradient force) [32].Therefore, the forces of gravity, resistance, and buoyancy on particles were considered in this simulation experiment.
The force on particles in fluids can be expressed as follows: where F f ,p is the overall force of the fluid on the particle, F g is the gravity of the particle, and F d and F f are the dragging force and buoyancy of the particle, respectively. Figure 1 shows the force on the particle in fluids.
Under the influence of fluid, the dragging force F d on the particle can be calculated as follows [33,34]: where C d is the dragging coefficient, and C d = 12μ f /rρ f u p ; ρ f is the density of the fluid; r p is the radius of the particle; μ f is the fluid viscosity coefficient; u f and u p are the velocity vectors of the fluid and particles; and f ðεÞ is a function related to the porosity. The buoyancy of the particle in fluids can be calculated as follows: where g is the acceleration of gravity.

Movement Velocity of Particles in Fluids.
Spherical particles can be used to analyze the movement of solids in fluids. The existing studies have found that the movement velocity of particles increases with the increase of radius, gravity, and the density difference between the fluid and the solid and decreases with the increase of the fluid viscosity coefficient [35]. The movement velocity of spherical particles in the fluid can be expressed by where u p is the movement velocity of the particle.

Calculation of Seepage Coefficient.
To ensure the accuracy and speed of calculation, a cube method was used to 2 Geofluids monitor the porosity of the coal measure sandstone during loading in the simulation model. For the cube method, when the particle is in different fluid units, its volume in a certain fluid unit can be calculated by the volume of the converted cube of the particle. Figure 2 shows the specific calculation method, in which the cube represents fluid elements and the circles represent particles. As shown in Figure 2, the volume of particle A covering different fluid elements can be calculated by the cube method, while the volume of particle B in one fluid element can be directly calculated. The mainly step of the cube method is converting the particle with radius R which is not at the same fluid element into a cube with the same position and the side length 2R, as shown in Figure 2. Then, the volume of the particle in the calculated fluid element can be expressed as follows: where V pf is the volume of the particle in the calculated fluid element, V cf is the volume of the overlap part of the transformed cube and the calculated fluid element, V p is the volume of the calculated particle, and V c is the volume of the transformed cube. The permeability coefficient of a certain fluid element m can be obtained as follows: where ε is the matrix porosity.
To ensure the rationality of the permeability coefficient and the continuity in the calculation of the numerical model, if the porosity is greater than 0.7, the porosity was also set as 0.7 in Equation (6). sandstone is affected by the chemical corrosion of fluids, long-term in situ stress, and other corrosive factors. Aiming at the creep behavior of surrounding rock in the deep mining area, the concept of stress corrosion was introduced, and the PSC model was established based on the parallel bonding of particle discrete elements. The principle of the PSC model was based on the phenomenon that the stress corrosion at the bond of the particle rock mass model preferentially occurred in the chemically active pore fluid and the higher stress area near the crack tip [29]. According to the theory of fracture mechanics, the PSC model was introduced into the original PFC-CFD model of coal-measure sandstone.

Creep
Based on the linear elastic fracture mechanics equation, the derivative of the bond radius of the numerical model with time can be expressed as follows [29]: where D is the bonded radius, σ is the stress of the calculated contact bond, σ c is the tensile strength of the bond, σ a is the threshold stress, and β 1 and β 2 are constants of the material and change with the temperature and chemical environment.

Sensitivity Analysis of Parameters in Creep Equation.
The sensitivity analysis of parameters β 1 , β 2 , and σ a was carried out to obtain their values in this study, and the results were shown in Figures 3 and 4. When one of the three parameters was analyzed, the other two parameters were fixed in this study. Referring to previous studies [36] and calibration results, the values of β 1 , β 2 , and σ a were finally determined as 5 × 10 −14 , 21, and 20 MPa, respectively. Figures 3 and 4 show the relationship between β 1 , β 2 , and σ a and the failure time t f of the model. As shown in Figure 3, β 1 and the logarithm of damage time t f can be fitted by a polynomial function lg ðt f Þ = 4:50 − 0:37β 1 + 0:014β 1 2 , and the determination coefficient R 2 is 0.98. With the increase of β 1 , the failure time of the model decreases sharply. Through the analysis in Figure 3 and Equation (7), when β 1 is smaller than 5:86 × 10 −14 , the rate of bond radius corrosion and damage development of the numerical model is slow, and the logarithm of failure time is more than 2.74.
As shown in Figure 3, the logarithmic fitting of β 2 and damage time t f can be fitted by the linear relation lg ðt f Þ = 10:4 − 0:36β 2 , and the determination coefficient is R 2 = 0:99. Besides, with the increase of β 2 , the failure time sharply decreases. Through the analysis in Figure 3 and Equation (7), as β 2 increases, the bond radius corrosion rate of the numerical model increases exponentially, and the failure time increases rapidly.
As shown in Figure 4, there is a significant correlation between σ a/ σ c and 0.8 correlation, while an insignificant 3 Geofluids correlation between σ a and damage time t f . Therefore, the relationship between 0.8-σ a/ σ c and the logarithmic of failure time t f is investigated. It can be seen from Figure 4 that the closer the value of σ a/ σ c to 0.8, the more sharply the failure time increases; when the value is far from 0.8 gradually, the failure time tends to be stable.
There is a nonlinear relationship between 0.8-σ a/ σ c and the logarithm of the failure time t f , which can be fitted by an exponential function. The fitting equation can be expressed as lg t f À Á = 0:290 × exp 0:322 According to the above sensitivity analysis, different values of β 1 , β 2 , and σ a were selected to assign to the PSC model. Based on the crack initiation time and multiple parameters adjustment, β 1 , β 2 , and σ a are determined to be 5 × 10 −14 , 21, and 20 MPa for the further simulation in this study. The values of these parameters are basically consistent with the data in the previous study on seepage damage behavior of coal measure sandstone [37].

Establishment of the Numerical Simulation Model.
Since there is no specific conversion formula between macroparameters and mesoparameters, the trial algorithm is widely used to calibrate the mesoparameters of the proposed models [38,39]. The calibration steps are as follows: the  To effectively calibrate the mesoparameters of the coal measure sandstone model, a cuboid model with a length of 50 mm, a width of 50 mm, and a height of 100 mm was established in this study. Table 1 shows the macroparameters for calibration. Table 2 shows the mesoparameters after calibration. The fluid simulation was performed by the CFD, while the rock mass simulation was performed by the particle discrete element. In order to prevent the particle model from exceeding the fluid boundary after deformation, the fluid domain was set larger than the boundary of the granular rock skeleton, and the rock skeleton was completely wrapped by the fluid domain. Figure 5 shows the schematic diagram of the coal measure sandstone model under the seepage-creep coupling effect by the PFC-CFD-PSC method.
The meshing of the fluid domain in Figure 5(a) should satisfy the following condition: where l is the shortest side length of the fluid and Δa is the side length of the fluid element. For this simulation test, the size of the fluid domain was designed as 60 mm × 60 mm × 100 mm, and the size of each fluid cell was 5 mm × 5 mm × 5 mm. The size of the coal measure sandstone particle model was 50 mm × 50 mm × 100 mm, and the fluid domain transparency in Figure 5 was 80%.

Fracture Extension Characteristics of Seepage-Creep
Numerical Model. The coal measure sandstone model in the numerical simulation was a 3D model, and the number of microcracks in the final failure model was large. Therefore, it was difficult to distinguish the macroscopic penetration cracks and the local damage in the entire 3D model. To intuitively observe the fracture characteristics of the model, 2D slices were obtained by cutting a 3D fracture map along the vertical direction. Figure 6 shows the cracks generated in coal measures sandstone under different seepage pressures (P in is the seepage pressure). Figure 6 shows the fracture morphology of coal measure sandstone under triaxial compression with a confining pressure of 4 MPa and residual strength of 70% peak strength. In the failure model, the shape and trend of the main cracks are drawn by red lines, and the thickness of the red line represents the size of the cracks. It can be found that under the same conditions, the larger the seepage pressure of the rock mass model, the more the penetration cracks, and the greater the degree of macroscopic damage. When the seepage pressure is 0.5 MPa, obvious inclined cracks appear inside the sample, and the failure mode is a typical shear failure (Figure 6(a)). When the seepage pressure is 1.5 MPa, there are two penetration cracks generated in the rock model. Specifically, one penetration crack was parallel to  the axial loading direction of the rock model, and its failure mode is a tensile failure, while the failure mode of the other penetration crack is a typical shear failure. Therefore, the failure mode of the seepage-creep rock model under the seepage pressure of 1.5 MPa is the composite shear-tensile failure. When the seepage pressure is 2.5 MPa, multiple

Pressure Gradient Distribution of Seepage Grid.
To display the internal hydraulic gradient of the 3D coal measure sandstone model, the planes along Y = 25 mm and Z = 9 mm were obtained, as shown in Figure 7. Figure 7 shows a schematic diagram of the hydraulic gradient in the fluid domain of the particle discrete element seepage-creep model under different seepage pressures. Before the calculation process of the seepage-creep numerical model, there are no particles in the edged fluid grid, that is, the coal measure sandstone only exists in the red box in Figure 7. Therefore, the rock mass has little influence on the hydraulic gradient at the edge grids of the seepage field.

Geofluids
As shown in Figure 7(a), when the seepage pressure is 0.5 MPa, the pressure gradient of the seepage field around and in the center of the rock is smaller; the larger pressure gradient is distributed in a circle at the top of the rock, and the value of pressure gradient ranges from 0.015 to 0.017. As shown in Figure 7(b), when the seepage pressure is 1.5 MPa, the pressure gradient of the seepage field around and in the center of the rock is small; the larger pressure gradient is also distributed within a circle at the top of the rock, and the value is between 0.085 and 0.094. As shown in Figure 7(c), when the seepage pressure is 2.5 MPa, there is a small pressure gradient of the seepage field around the rock top, ranging from 0.25 to 0.32; the value of pressure gradient at the top central part of the rock is between 0.35 and 0.51, while the value of the four corners is between 0.30 and 0.35. As shown in Figure 7(d), when the seepage pressure is 3.5 MPa, there is a small pressure gradient of the seepage field around the top of the rock, ranging from 0.35 to 0.6, while the pressure gradient of the middle top part of the rock is between 0.75 and 0.93. The high hydraulic gradient values of the model are concentrated at the top of the rock. Due to the obstructive effect of particles on the fluid flow, the hydraulic gradient sharply decreases at the beginning of the water flow. Therefore, the seepage effect has a greater impact on the top of the rock model. As shown in Figure 7, the seepage at the four top corners of the coal measure sandstone has a smaller influence on the rock particles.
Besides rock models under seepage pressures of 0.5 MPa, 1.5 MPa, 2.5 MPa, and 3.5 MPa, the rock models under seepage pressures of 1.0 MPa, 2.0 MPa, and 3.0 MPa were also simulated for the further analysis. The relationship between the highest-pressure gradient and the seepage pressure is obtained, as shown in Figure 8.
As shown in Figure 8, with the increase of the seepage pressure, the maximum pressure gradient also increases. The increasing rate was much smaller within the seepage pressure of 0.5-1.5 MPa, and the increasing value of seepage pressure was 0.078. When the seepage pressure is greater than 1.5 MPa, the maximum pressure gradient increases linearly, and the increasing rate is accelerated.

Conclusion
A seepage-creep coupling model based on the PFC-CFD-PSC method is proposed to explore the mutual influence between the seepage field and the creep field in this study, and the proposed model is applied to analyze the seepagecreep behavior of coal measure sandstone.
(1) The simulation results show that the seepage pressure has a significant effect on the damage degree of sandstone. The greater the seepage pressure on the coal measure sandstone model, the more penetration cracks are generated inside the rock and the greater the degree of macroscopic damage (2) When the seepage pressure is 0.5 MPa, the failure mode of the coal measure sandstone is a typical shear failure; when the seepage pressure is 1.5 MPa, the failure mode of coal measure sandstone is a composite tensile-shear failure; when the seepage pressure is 2.5 MPa, a shear failure can be observed; when the seepage pressure is 3.5 MPa, the obvious main crack mode presents a shear failure. When the residual strength is 70% of the peak strength,

Data Availability
All data, models, and code generated or used during the study appear in the submitted article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.

Acknowledgments
This work was supported by the National Natural Science Foundation of China (NSFC) (51974296, 52061135111).