Abstract
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 seepage-creep 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.5 MPa and increases sharply under the seepage pressure of greater than 1.5 MPa.
1. Introduction
Deep coal mining has been widely used in the mining areas of China [1–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–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–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–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–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–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 parallel-bonded 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 parallel-bounded 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.
2. 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.
2.1. 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 is the overall force of the fluid on the particle, is the gravity of the particle, and and 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 on the particle can be calculated as follows [33, 34]: where is the dragging coefficient, and ; is the density of the fluid; is the radius of the particle; is the fluid viscosity coefficient; and are the velocity vectors of the fluid and particles; and is a function related to the porosity.
The buoyancy of the particle in fluids can be calculated as follows: where is the acceleration of gravity.
2.1.2. 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 is the movement velocity of the particle.
2.1.3. Calculation of Seepage Coefficient
To ensure the accuracy and speed of calculation, a cube method was used to 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 covering different fluid elements can be calculated by the cube method, while the volume of particle in one fluid element can be directly calculated. The mainly step of the cube method is converting the particle with radius which is not at the same fluid element into a cube with the same position and the side length , as shown in Figure 2. Then, the volume of the particle in the calculated fluid element can be expressed as follows: where is the volume of the particle in the calculated fluid element, is the volume of the overlap part of the transformed cube and the calculated fluid element, is the volume of the calculated particle, and is the volume of the transformed cube.
The permeability coefficient of a certain fluid element 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).
2.2. Creep Simulation in the PFC-CFD-PSC Model
2.2.1. Realization of Particle Creep Behaviors
With the passage of time, the seepage damage of the coal measure 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.
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 is the bonded radius, is the stress of the calculated contact bond, is the tensile strength of the bond, is the threshold stress, and and are constants of the material and change with the temperature and chemical environment.
2.2.2. Sensitivity Analysis of Parameters in Creep Equation
The sensitivity analysis of parameters , , and 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 , , and were finally determined as , 21, and 20 MPa, respectively. Figures 3 and 4 show the relationship between , , and and the failure time of the model.


As shown in Figure 3, and the logarithm of damage time can be fitted by a polynomial function , and the determination coefficient is 0.98. With the increase of , the failure time of the model decreases sharply. Through the analysis in Figure 3 and Equation (7), when is smaller than , 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 and damage time can be fitted by the linear relation , and the determination coefficient is . Besides, with the increase of , the failure time sharply decreases. Through the analysis in Figure 3 and Equation (7), as 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 and 0.8 correlation, while an insignificant correlation between and damage time . Therefore, the relationship between 0.8- and the logarithmic of failure time is investigated. It can be seen from Figure 4 that the closer the value of 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- and the logarithm of the failure time , which can be fitted by an exponential function. The fitting equation can be expressed as
According to the above sensitivity analysis, different values of , , and were selected to assign to the PSC model. Based on the crack initiation time and multiple parameters adjustment, , , and are determined to be , 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].
2.3. 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 macroscopic law of the numerical test is compared with the indoor experimental results, and then, the mesoparameters are adjusted to ensure that the macroscopic law of the numerical simulation is basically consistent with the experimental value.
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.

(a) Fluid domain

(b) Coal measure sandstone in the fluid
The meshing of the fluid domain in Figure 5(a) should satisfy the following condition: where is the shortest side length of the fluid and is the side length of the fluid element.
For this simulation test, the size of the fluid domain was designed as mm, and the size of each fluid cell was . The size of the coal measure sandstone particle model was , and the fluid domain transparency in Figure 5 was 80%.
3. Results and Discussion
3.1. 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 ( is the seepage pressure).

(a) MPa

(b) MPa

(c) MPa

(d) MPa
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 cracks are generated in the failure model, among which there were three parallel inclined main cracks. The shear failure is the main failure mode of the rock model under the seepage pressure of 2.5 MPa. The final damage degree of the rock model under the seepage pressure of 3.5 MPa is much larger than that under the seepage pressure of 2.5 MPa. Besides, the number of main cracks in the coal measure sandstone under the seepage pressure of 3.5 MPa is the largest, and the shear failure is the main failure mode of the coal measure sandstone under the seepage pressure of 3.5 MPa. It can be concluded that the seepage pressure has a significant effect on the failure mode of the sandstone. With the increase of seepage pressure, the number of penetration cracks in the failure model also increases. When the seepage pressures are 0.5 MPa, 1.5 MPa, 2.5 MPa, and 3.5 MPa, the number of the broken bond of the rock model is 11868, 18372, 19041, and 20377, respectively. It indicates that the increase of seepage pressure can increase the damage degree of the seepage-creep model.
3.2. Pressure Gradient Distribution of Seepage Grid
To display the internal hydraulic gradient of the 3D coal measure sandstone model, the planes along mm and mm were obtained, as shown in Figure 7.

(a) MPa

(b) MPa

(c) MPa

(d) MPa
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. 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.
4. 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 seepage-creep 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, the numbers of bond failures in the rock model under seepage pressures of 0.5 MPa, 1.5 MPa, 2.5 MPa, and 3.5 MPa are 11868, 18372, 19041, and 20377, respectively(3)The simulation results show that the high hydraulic gradient is mainly concentrated at the top of the rock. The seepage at the top rock has a greater impact on the rock model, while there is a smaller fluid force at the top of the four corners of the rock. The maximum value of the fluid grid pressure gradient increases with the increase of the seepage pressure, and it increases slowly under the seepage pressure of less than 1.5 MPa and increases sharply under the seepage pressure of greater than 1.5 MPa
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).