Numerical Analysis of the Mud Inflow Model of Fractured Rock Mass Based on Particle Flow

Zhejiang Engineering Survey and Design Institute Group Co. Ltd, Ningbo, 315010 Zhejiang Province, China Work Safety Key Lab on Prevention and Control of Gas and Roof Disasters for Southern Coal Mines, Hunan Provincial Key Laboratory of Safe Mining Techniques of Coal Mines, Hunan University of Science and Technology, Xiangtan, 411201 Hunan Province, China School of Resource Environment and Safety Engineering, Hunan University of Science and Technology, Xiangtan, 411201 Hunan Province, China School of Civil Engineering, Chongqing University, Chongqing 400045, China


Introduction
Engineering construction under karst geological conditions often causes various disasters, such as water inrush, water gushing, and collapse. Among them, the problems of water inrush and mud outburst are particularly significant in fractured rock masses [1,2]. Due to the complex dynamic process of water inrush and water gushing formation and time, the audience is influenced by many factors, such as formation lithology, hydrogeological conditions, construction methods, and change of mechanical behavior of rock mass [3]. The erosion of fluid reduces the strength and stability of engineering rock masses and brings serious threats to project construction and operation.
For underground engineering, the rock mass excavation not only involves the formation of a new equilibrium state after the original equilibrium state is broken [4] but also the rock masses in a certain range around the cavern will become loose under the influence of joints, fractures, and other structural surfaces slip and collapse. For this case, the particle discrete element PFC can be better simulated [5][6][7][8][9]. Yang et al. [10] performed a discrete element simulation of the fractured red sandstone under uniaxial compression and revealed the failure mechanism of the fractured red sandstone; Fan and Cao [11] established a numerical calculation model containing two fractures based on the particle contact adhesion model in the discrete element numerical analysis software PFC3D and discussed the effect of rock bridge angle on the mechanical properties of defective rock; Liu et al. [12] studied the crack expansion of indented rock specimens numerically; Potyondy and Cundall [13] proposed a rock bonded particles numerical model (BPM) and used the discrete element numerical analysis software PFC to simulate; Zhu and Xu [14] considered the effects of blasting disturbance, geological conditions of surrounding rock, and tunnel cross-section form on tunnel collapse and carried out a series of numerical simulation calculations, and a series of numerical simulation calculations were carried out, the calculated tunnel collapse amount was compared with the commonly used loose load calculation results, and the closest calculation value was obtained. Furthermore, many scholars [15][16][17] conducted a great deal of research using PFC particle discrete elements in fluid-solid coupling calculations; Zhou et al. [18,19] based on particle flow theory using the FISHTANK function library of the PFC software calculation program built-in FISH language successfully simulated the seepage and piping of sand and cohesive soil materials and obtained the variation law of pressure and flow velocity in the seepage process; Luo et al. [20] used the PFC particle flow platform and the builtin FISH language to compile a program to simulate the law of particles falling freely in the fluid and the change law of seepage gradient and flow rate under different pressure differences, and the results basically conform to Darcy's law; Bai [21] simulates the whole process of foundation pit water inrush through PFC2D; Wang et al. [22] used PFC3D numerical software to study and explore the influence of factors such as fault water pressure and rock fracture properties on the water inrush from tunnel and mud inrush.
Due to the complex nature of mud inflow seepage and the uncertainty of boundary conditions often make the numerical calculation of seepage impact far from the actual test results, and different numerical software has different scope of application, such as finite difference software for continuous nonlinear multifield coupling problems (FLAC3D), discrete element software has unique advantages in solving structural surface control problems (3DEC), and particle flow software is particularly suitable for brittle material fracture development and bulk flow deformation problems (PFC3D). Considering that the rock mass is a discontinuous medium to simulate the problem of large deformation of rock mass in underground engineering, PFC numerical software has the advantages of reasonable treatment of permeability boundary and effective application of multiple flow-solid coupling models. This paper uses the discrete element method, using particle flow software (PFC3D), combined with fluid dynamics theory, establishes a fractured rock mass mud inflow model, simulates the whole process of fractured rock mud inflow, and discusses the impact of particle size and flow rate on the change of pressure gradient, which provides positive guidance significance for underground engineering construction in karst areas.

Basic Assumptions.
Since there is no real fluid in the PFC3D program, a fluid unit containing a certain number of particles is used instead of a fluid for the calculation. The flow process and law of the particles in the fracture and the indicators of pressure, velocity vector, and flow velocity on the fluid unit are obtained by the action of the particles rolling or flowing in the fracture. To facilitate the modeling of mud inflow in a single fracture rock, the following assumptions are proposed: (1) The force generated by fluid flow acts on the particle in the form of physical force and is applied to the fluid in the same magnitude (2) Particle units in the model are considered as rigid bodies, and the contact between particles occurs only in a small area, approximating point contact The contact behavior is flexible, allowing a certain amount of superposition between the rigid particles at the contact, but the value is much smaller than the radius of the particles, and the contact force is linked to the superposition between the particles through the force-displacement law  Figure 1 shows a fixed control volume, where ΔV = Δx × Δy × Δz, and the number of particles in the fluid unit is n p . Assuming that the fluid flows only in the x-direction and there is a pressure gradient dp/dx in that direction, so that the forces on the particles inside the fluid element in the seepage direction are balanced, the total permeability on the particles is expressed as [23][24][25]: where f int x is the interaction force between the fluid and solid phases in the unit volume, d pi is the particle diameter, p is the seepage pressure, the negative sign of the first term on the right side of the equation indicates that the force acting on the fluid is positive, and the negative sign of the second term indicates that the pressure is gradually decreasing along the x-positive direction. i = 1, 2, 3, ⋯, n p , j = x, y, z.
Considering the porosity, it can be expressed as: Substituting Eq. (2) into Eq. (1) yields: Thus, the total force acting on the particle is The interaction force between the particle and fluid under fluid-solid coupling can be expressed as: Substituting Eq. (5) into Eq. (4), the total force on the particle under the action of fluid-solid coupling can be obtained as: 2.2.2. Pressure Gradient in Pore Media. The pressure gradient is used in pore media to represent the fluid-solid coupling interaction, which can usually be calculated by empirical equations. Darcy's law is widely used in the fluid phase of porous media, and the pressure gradient in the expression of Darcy's law is proportional to the apparent relative velocity as follows.
where ν f , ρ f is the kinematic viscosity and density of the fluid, respectively, u x0 is the apparent relative velocity, K is the permeability coefficient, and k is the permeability. According to the Kozeny-Carman equation, the permeability of the pore medium can be expressed as: where c is the Kozeny-Carman constant, generally taken between 0.003 and 0.0055; d p is the average diameter of the particles. Darcy's law applies to laminar flows with Reynolds coefficients of 1~10. For relatively large Reynolds coefficients and nonlinear equations, the following equation can be used.
By substituting nu x = u x0 into Eq. (9) yields: Since the particle forces are generated by relative fluid flow, the relative velocity u rx = υ x − u x (where υ x is the average velocity of particles in a given control area) is replaced by the absolute fluid velocity, which following formula can be obtained: Equation (11) can be expressed as: For large porosity (n ≥ 0:8), the pressure gradient can be calculated by the following equation.
where C d is the drag coefficient of the sphere, which is a function related to the Reynolds coefficient R ep :

Navier-Stokes Equation and Continuity Equation.
The Navier-stokes equation for a fluid in liquid-phase [24], nonturbulent, viscous incompressible fluid-solid two-phase, can be expressed as 3 Geofluids where u is the flow velocity vector, τ is the viscous stress tensor, g is the acceleration of gravity, f int is the interaction force between the particle and the fluid in a unit volume, and ρ f is the fluid density, where the viscous stress tensor τ can be expressed as: where μ and e •ðdÞ are the viscosity and the stress deflection tensor, respectively. Expanding Eqs. (15) and (16) in the Cartesian coordinate system, we obtain: The unit volume solid-liquid interaction force can be expressed as: where the solid-liquid friction coefficient can be calculated by the following equation where β int j friction coefficient is defined as fluid, combining Eqs. (21) and (22), we can obtain that:    Table 2.
3.1.1. Generated Particles. Since the model generation sequence and method are the same, only one example is used for illustration here. First, 6 walls are generated as faces, and the right-hand rule is used in PFC3D to determine the effective face of the wall (four fingers surround the 4 coordinate points of the wall, and the direction pointed by the thumb is the effective face); the 6 walls generated that are composed of dimensions is 0:3 × 0:1 × 0:1 m (Length × Width × Height), as shown in Figure 2.
Then, the particles are generated in the hollow cuboid, the radius distribution of particles in PFC3D has normal distribution, uniform distribution, Gaussian distribution, etc. In order to better simulate the inhomogeneity of particle distribution, this numerical simulation adopts uniform distribution in random distribution, the distribution of particle radius R adopts from R min to R max uniform distribution, and the average radius of particles is R = ðR min + R max Þ/2. Considering the calculation accuracy and calculation time, the minimum radius of particles R min = 0:05 mm, the maximum radius R max = 1 mm, and the number N of particles is determined by the porosity n.
where n is the porosity; V is the cuboid volume, m 3 ; and R is the average particle radius, m.
In order to make the particles randomly generated in the cuboid with a given porosity, ensure the uniformity of the particle distribution, and make the generated particles reach the initial equilibrium state in a short time, the particles are randomly generated by the expanding radius method, and the gravity method was used to make the particles reach equilibrium, as shown in Figure 3.

Generate
Fractures. Delete one wall in the x-direction, while establishing multiple walls on that side to form a group, the gap width between the wall, and the wall d = 6R max , and the gap width is not less than 6 mm that a fracture; according to the fracture width d = 6R max to establish the wall perpendicular to the wall group, the fracture width used in this paper are 8 mm, as shown in Figure 4.

Set up Fluid Unit.
A total of 250 fluid units are set in this model, and the unit size is 2 × 1 × 1 mm. The boundary condition of the side wall is smooth, and the fluid flow velocity parallel to the wall surface is not equal to 0. The water pressure P is applied at x = 0, and the water pressure P is equal to 0 at x = 0:3, forming a water pressure difference △ P. For computational simplicity and rapid convergence, the fluid time step is set to 5:0 × 10 −4 s, as shown in Figure 5.

Geofluids
3.2. Initial Model Balance. The viscous damping coefficient is set at the solid connection for the energy consumption at the connection between the particles, and the ratio of the normal and tangential damping constants to the critical damping constant is 0.1. The friction coefficient between the particles and the particles is 0.5, and the friction coefficient between the particles and the wall is 0.3. Numerically simulated material fine views parameters are shown in Table 1. The particles reach the equilibrium state under the action of gravity, and then the fracture is generated by deleting the wall and establishing the wall, the width of the fracture should exceed 6 times the radius of the largest particle, and water and particles are discharged from this fracture. The changes of water velocity, porosity, particle loss, and other parameters are also tracked and recorded. The specimen at a moment of particle equilibrium disruption is shown in Figure 6, where yellow represents clay, green represents the magnitude of the particle normal contact force, and red arrows represent the direction of fluid flow.

Solution Method.
During the calculation of fluid circulation, Newton's second law and force-displacement law are The particles move in the middle of the model Figure 8: The vortex phenomenon and particle movement trend in the initial stage.
Fracture Figure 9: The change of the particle flow direction. 6 Geofluids applied repeatedly to the particles. Newton's second law is used to update the position of the particles and the wall to readjust the contact relationship between the particles, while the force-displacement law is used to update the contact forces in the contacting parts. The two alternate, iterating, and traversing the entire particle set in time steps until equilibrium is reached or damage occurs, and the pressure and velocity vectors within each fluid unit are calculated by the semi-implicit PLE algorithm. The calculation process is shown in Figure 7. The conditions for convergence are set as follows.
when ratio sum and ratio max are both less than 0.0001, stop calculation.

Analysis of Flow Law of Particle Seepage.
According to the numerical simulation results, it can be found that most of the particles move to the middle of the model when the water pressure is first applied, and a very small number of particles immediately enter the fracture, meanwhile, the movement of particles at the corner of the wall will produce a phenomenon similar to vortex flow ( Figure 8). As the numerical simulation calculation continues to advance, the particles gradually enter the fracture, the fluid viscosity keeps decreasing, the flow direction of the particles will change, the vortex flow phenomenon disappears, the particles at the corner of the wall flow toward the fracture, and the rest of the particles flow along the vertical direction ( Figure 9). In addition, the outflow of particles, at the position x = 0 at the bottom of the model, is found to rise as a whole, which is the result of activating the buoyancy factor command during the simulation and setting the density of the fluid to be comparable to the density of the particles. Using the PFC3D program, com-piled in the FISH language, the volume of flowing particles in the fracture was recorded in real time at different pressure, and the results of the numerical simulation are shown in Figure 10.
According to Figure 11, it can be found that the seepage flow of particles under different water pressure all increase with time, such as model size 100 × 100 × 300 mm; with the increase of pressure gradient, the percolation volume of particles after 0.45 s are 4:4 × 10 −4 m 3 , 6:9 × 10 −4 m 3 , 9:2 × 10 −4 m 3 , and 1:2 × 10 −3 m 3 , and its adjacent two pressure gradient the percolation volume of particles under the difference increased by 56%, 33%, and 30%, respectively, with the model size of 100 × 100 × 250 mm; and the percolation volume of particles after 0.35 s was 4:7 × 10 −4 m 3 , 6:9 × 10 −4 m 3 , 8:5 × 10 −4 m 3 , and 1:4 × 10 −3 m 3 with the increase of pressure gradient, and the percolation volume of particles under the difference of its two adjacent pressure gradients were increased by 47%, 23%, 64%, and the model size of 100 × 100 × 200 mm; with the increase of pressure gradient, the percolation volume of particles after 0.25 s was 3:7 × 10 −4 m 3 , 6:9 × 10 −4 m 3 , 9:2 × 10 −4 m 3 , and 1:3 × 10 −3 m 3 , and the percolation volume of particles under the difference of its two adjacent pressure gradients increased by 86%, 33%, and 42%. In addition, the model size and running time will increase the seepage volume of particles; in the case of the model size is large enough and the running time is relatively short, the seepage volume of particles increases and decreases with the increase of pressure gradient, which is due to the fact that when the particles enter the fracture, the width of the fracture remains the same, and the increase of the number of particles will make the crowding and friction between the particles and the friction between the particles and the wall, causing a certain degree of fracture. In the case of small model size and long running time, most of the particles enter the fracture in that running time, and the increase of seepage volume of particles will decrease first and then increase with the increase of pressure gradient, which indicates that the blocked particles in the fracture decrease, the porosity of the model increases, and the mobility of particles is enhanced, so the faster the speed of entering the fracture, the greater the increase of seepage volume of particles.  13 Geofluids 3.4.2. Analysis of Particle Seepage Rate. The particle percolation problem in the fracture is not really a water percolation problem in the fracture, and water is a continuous medium flowing in the fracture that almost satisfies Darcy's law, while the motion of particles is approximately regarded as the motion between small balls, which collide with each other and there is friction; in addition to the characteristics of the flow, the small ball itself also has the nature of rotation. Therefore, the most fundamental difference between water and particle flow in a single fracture is that water is a continuum and particle is a discrete body, by compiling the FISH language program and using the Hist command to record the amount of particle percolation in the fracture at different running times, and using Origin software for fitting, the fitting results are shown in Figure 12. Figure 12 shows that in the initial stage, the particle flow rate increases with the increase of running time, and the slope of the curve increases; consider fitting with a quadratic function, the fitting correlation coefficientR 2 is above 0.92, and the slope of the curve remains unchanged after the percolation rate of particles reaches stability; consider fitting with a primary function, the fitting coefficientR 2 is mostly above 0.90; and when the running time is long enough, the model after all the particles in the model are seeped out by the fracture, the particle flow rate is kept constant, and the slope is 0, as shown in Figure 12(l).
To analyze the relationship between pressure gradient and seepage velocity, the slope of the primary function (seepage velocity) was fitted for each model with different pressure gradients, and Figure 13 shows the fitted curves of pressure gradient and particle speed. Figure 13 shows that the flow velocity is not exactly proportional to the pressure gradient (data points), and it does not conform to Darcy's law. There are various reasons for this, such as changes in fluid viscosity, changes in model porosity, and collisions between particles and between particles and fracture walls or even the clogging effect of particles on fractures. When the pressure gradient is small, the flow velocity increases slowly with the increase of pressure gradient, after the pressure gradient exceeds a certain value, the flow velocity increases rapidly with the increase of gradient, and finally the region increases slowly, and the whole shows an "S" curve of growth. By further fitting the data (see the red curve in Figure 13), the pressure gradient and flow velocity approximately satisfy an exponential function relationship, and the fitted correlation coefficient R 2 equal to 0.9312, and the fitting formula is as follows: 2E6 3E6 4E6 5E6 6E6 7E6 Figure 13: The relationship between pressure gradient and the particle speed.

Conclusion
(1) When the water pressure is first applied, the flow direction of the particles will show different, a small number of particles immediately into the fracture, the vast majority of particles are flowing toward the middle of the model, the movement of particles at the corners of the wall occurs similar to the phenomenon of vortex, with the increase in running time, the vortex phenomenon disappears, the flow direction of the particles at the corners point to the fracture, the rest of the particles flow along the vertical direction (2) In the initial stage, the flow rate of particles increases with time, the slope of the curve increases, and the quadratic function is used for fitting. After the percolation rate of particles reaches stability, the slope of the curve remains the same, and the primary function is used for fitting (3) In the beginning stage, the particle flow velocity increases with the increase of the pressure gradient. In addition, the particle flow rate and pressure gradient are influenced by many factors, such as the change of fluid viscosity, the change of model porosity, and the collision between particles, particles and fracture wall, so that they approximately satisfy an exponential function of "S" curve

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

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