Abstract
Water inrush and mud outburst are one of the crucial engineering disasters commonly encountered during the construction of many railways and tunnels in karst areas. In this paper, based on fluid dynamics theory and discrete element method, we established a fractured rock mass mud inflow model using particle flow PFC3D numerical software, simulated the whole process of fractured rock mass mud inflow, and discussed the effect of particle size and flow velocity on the change of pressure gradient. The numerical simulation results show that the movement of particles at the corner of the wall when the water pressure is first applied occurs similar to the vortex phenomenon, with the running time increases, the flow direction of particles changes, the vortex phenomenon disappears, and the flow direction of particles at the corner points to the fracture; in the initial stage, the slope of the particle flows rate curves increases in time, and the quadratic function is used for fitting. After the percolation velocity of particles reaches stability, the slope of the curve remains constant, and the primary function is used for fitting; the particle flow rate and pressure gradient are influenced by a variety of factors, and they approximately satisfy the exponential function of an “S” curve.
1. 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–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–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 built-in 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.
2. Mud Inflow Theory
2.1. 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(3)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(4)The walls are parallel to each other, smooth and the width of the fracture can be adjusted by changing the distance between thin walls
2.2. Mud Inflow Theory and Its Implementation
2.2.1. The Force on the Particle
Figure 1 shows a fixed control volume, where , and the number of particles in the fluid unit is .
Assuming that the fluid flows only in the -direction and there is a pressure gradient 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–25]: where is the interaction force between the fluid and solid phases in the unit volume, is the particle diameter, 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 -positive direction. , .
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 is the kinematic viscosity and density of the fluid, respectively, is the apparent relative velocity, is the permeability coefficient, and is the permeability.
According to the Kozeny-Carman equation, the permeability of the pore medium can be expressed as: where is the Kozeny-Carman constant, generally taken between 0.003 and 0.0055; 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 into Eq. (9) yields:
Since the particle forces are generated by relative fluid flow, the relative velocity (where 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 (), the pressure gradient can be calculated by the following equation. where is the drag coefficient of the sphere, which is a function related to the Reynolds coefficient :
2.2.3. 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 where is the flow velocity vector, is the viscous stress tensor, is the acceleration of gravity, is the interaction force between the particle and the fluid in a unit volume, and is the fluid density, where the viscous stress tensor can be expressed as: where and 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 friction coefficient is defined as fluid, combining Eqs. (21) and (22), we can obtain that:
3. Numerical Analysis of Mud Inflow in Fractured Rock Masses
3.1. Establishment of a Numerical Model of Mud Inflow
The mud inflow model consists of walls and balls, where the particles are simulated by small balls with certain dimensions, the fractures are realized by setting the distance between the walls, and the fluid is simulated by fluid units. Since the particles in the model are not directly assigned to the macroscopic parameters of the material, as shown in Table 1. It is necessary to calibrate the parameters of the particles to simulate four different groups of water pressure, with water pressure of 0.5 MPa, 0.8 MPa, 1.0 MPa, 1.2 MPa, and three different model sizes, , , and , as shown in 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 (), 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 adopts from to uniform distribution, and the average radius of particles is . Considering the calculation accuracy and calculation time, the minimum radius of particles , the maximum radius , and the number of particles is determined by the porosity . where is the porosity; is the cuboid volume, m3; and 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.
3.1.2. Generate Fractures
Delete one wall in the -direction, while establishing multiple walls on that side to form a group, the gap width between the wall, and the wall , and the gap width is not less than 6 mm that a fracture; according to the fracture width to establish the wall perpendicular to the wall group, the fracture width used in this paper are 8 mm, as shown in Figure 4.
3.1.3. Set up Fluid Unit
A total of 250 fluid units are set in this model, and the unit size is . 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 is applied at , and the water pressure is equal to 0 at , forming a water pressure difference . For computational simplicity and rapid convergence, the fluid time step is set to , as shown in Figure 5.
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.
3.3. Solution Method
During the calculation of fluid circulation, Newton’s second law and force-displacement law are 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 ratiosum and ratiomax are both less than 0.0001, stop calculation.
3.4. Analysis of Numerical Simulation Results
3.4.1. 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 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, compiled 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 ; with the increase of pressure gradient, the percolation volume of particles after 0.45 s are , , , and , 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 ; and the percolation volume of particles after 0.35 s was , , , and 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 , , , and , 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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
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 coefficientis 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 coefficientis 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 equal to 0.9312, and the fitting formula is as follows:
4. 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.
Acknowledgments
This research is supported by the National Natural Science Foundation of China (Nos. 51774131 and 51274097). The authors are thankful for all of the support for this basic research.