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

Authors: Yongjian Pan; Huajun Wang; Yanlin Zhao; Qiang Liu; Shilin Luo
Journal: Geofluids
 (2021)

---

## 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.


---

## Body

## 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

Figure1 shows a fixed control volume, where ΔV=Δx×Δy×Δz, and the number of particles in the fluid unit is np.Figure 1
Fluid through the particle assemble.Assuming that the fluid flows only in thex-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–25]:
(1)fdsum=∑i=1npfdix=−fintxΔV−dpdxπ6∑i=1npdpi3,where fintx is the interaction force between the fluid and solid phases in the unit volume, dpi 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,⋯,np, j=x,y,z.Considering the porosity, it can be expressed as:(2)n=1−1ΔVπ6∑i=1npdpi3.Substituting Eq. (2) into Eq. (1) yields:
(3)fdsum=∑i=1npfdix=−fintx1−n+dpdxπ6∑i=1npdpi3.Thus, the total force acting on the particle is(4)fdij=−fintj1−n+∇pjπ6dpi3.The interaction force between the particle and fluid under fluid-solid coupling can be expressed as:(5)fintj=n⋅∇pj.Substituting Eq. (5) into Eq. (4), the total force on the particle under the action of fluid-solid coupling can be obtained as:
(6)fdij=−∇pj1−nπ6dpi3.


#### 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.(7)dpdx=−ρfgKux0=−νfρfkux0,where νf,ρf is the kinematic viscosity and density of the fluid, respectively, ux0is 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:(8)k=c⋅dp¯2n31−n2,where c is the Kozeny-Carman constant, generally taken between 0.003 and 0.0055; dp¯ 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.(9)dpdx=−150⋅1−n2n3dp¯2νfρfux0+1.75⋅1−nn3dp¯ρfux02.By substitutingnux=ux0 into Eq. (9) yields:
(10)dpdx=−150⋅1−n2n2dp¯2νfρfux+1.75⋅1−nndp¯ρfux2.Since the particle forces are generated by relative fluid flow, the relative velocityurx=υ¯x−ux (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:
(11)dpdx=150⋅1−n2n2dp¯2νfρfux+1.75⋅1−nndp¯ρfυ¯x−ux⋅υ¯x−ux.Equation (11) can be expressed as:
(12)∇pj=150⋅1−n2n2dp¯2νfρfux+1.75⋅1−nndp¯ρfυ¯j−uj⋅υ¯j−uj.For large porosity (n≥0.8), the pressure gradient can be calculated by the following equation.
(13)∇pj=34ρfCd1−n⋅n−2.7d¯pυ¯j−uj⋅υ¯j−uj,where Cd is the drag coefficient of the sphere, which is a function related to the Reynolds coefficient Rep:
(14)Cd=24Rep1+0.15⋅Rep0.687Rep≤1000,0.44Rep≥1000,Rep=nυ¯j−ujd¯pυf.


#### 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
(15)∂n∂t=−∇⋅nu,(16)∂nu∂t=−∇⋅nuu−nρf∇p−nρf∇⋅τ+ng+fintρf,where u is the flow velocity vector, τ is the viscous stress tensor, g is the acceleration of gravity, fint 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:
(17)τ=μ⋅e•d,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:
(18)∂n∂t+∂nux∂x+∂nuy∂y+∂nuz∂z=0,(19)∂nux∂t+∂nux2∂x+∂nuxuy∂y+∂nuxuz∂z=−nρf∂p∂x−1ρf∂nτxx∂x+∂nτyx∂y+∂nτzx∂z+ngx+fintxρf,(20)∂nuy∂t+∂nuyux∂x+∂nuy2∂y+∂nuyuz∂z=−nρf∂p∂y−1ρf∂nτxy∂x+∂nτyy∂y+∂nτzy∂z+ngy+fintyρf,(21)∂nuz∂t+∂nuzux∂x+∂nuz2∂y+∂nuz2∂z=−nρf∂p∂z−1ρf∂nτxz∂x+∂nτyz∂y+∂nτzz∂z+ngz+fintzρf.The unit volume solid-liquid interaction force can be expressed as:(22)fintj=βintjυ¯j−uj,where the solid-liquid friction coefficient can be calculated by the following equation
(23)βintj=150⋅1−n2n⋅dp¯2uf+1.75⋅1−ndp¯⋅ρf⋅υ¯j−ujn≤0.8,34⋅ρf⋅CD⋅1−n⋅n−1.7dp¯⋅υ¯j−ujn>0.8,where βintj friction coefficient is defined as fluid, combining Eqs. (21) and (22), we can obtain that:
(24)fdij=−βintj1−nυ¯j−uj+∇pjπ6dpi3.


## 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

Figure1 shows a fixed control volume, where ΔV=Δx×Δy×Δz, and the number of particles in the fluid unit is np.Figure 1
Fluid through the particle assemble.Assuming that the fluid flows only in thex-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–25]:
(1)fdsum=∑i=1npfdix=−fintxΔV−dpdxπ6∑i=1npdpi3,where fintx is the interaction force between the fluid and solid phases in the unit volume, dpi 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,⋯,np, j=x,y,z.Considering the porosity, it can be expressed as:(2)n=1−1ΔVπ6∑i=1npdpi3.Substituting Eq. (2) into Eq. (1) yields:
(3)fdsum=∑i=1npfdix=−fintx1−n+dpdxπ6∑i=1npdpi3.Thus, the total force acting on the particle is(4)fdij=−fintj1−n+∇pjπ6dpi3.The interaction force between the particle and fluid under fluid-solid coupling can be expressed as:(5)fintj=n⋅∇pj.Substituting Eq. (5) into Eq. (4), the total force on the particle under the action of fluid-solid coupling can be obtained as:
(6)fdij=−∇pj1−nπ6dpi3.


### 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.(7)dpdx=−ρfgKux0=−νfρfkux0,where νf,ρf is the kinematic viscosity and density of the fluid, respectively, ux0is 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:(8)k=c⋅dp¯2n31−n2,where c is the Kozeny-Carman constant, generally taken between 0.003 and 0.0055; dp¯ 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.(9)dpdx=−150⋅1−n2n3dp¯2νfρfux0+1.75⋅1−nn3dp¯ρfux02.By substitutingnux=ux0 into Eq. (9) yields:
(10)dpdx=−150⋅1−n2n2dp¯2νfρfux+1.75⋅1−nndp¯ρfux2.Since the particle forces are generated by relative fluid flow, the relative velocityurx=υ¯x−ux (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:
(11)dpdx=150⋅1−n2n2dp¯2νfρfux+1.75⋅1−nndp¯ρfυ¯x−ux⋅υ¯x−ux.Equation (11) can be expressed as:
(12)∇pj=150⋅1−n2n2dp¯2νfρfux+1.75⋅1−nndp¯ρfυ¯j−uj⋅υ¯j−uj.For large porosity (n≥0.8), the pressure gradient can be calculated by the following equation.
(13)∇pj=34ρfCd1−n⋅n−2.7d¯pυ¯j−uj⋅υ¯j−uj,where Cd is the drag coefficient of the sphere, which is a function related to the Reynolds coefficient Rep:
(14)Cd=24Rep1+0.15⋅Rep0.687Rep≤1000,0.44Rep≥1000,Rep=nυ¯j−ujd¯pυf.


### 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
(15)∂n∂t=−∇⋅nu,(16)∂nu∂t=−∇⋅nuu−nρf∇p−nρf∇⋅τ+ng+fintρf,where u is the flow velocity vector, τ is the viscous stress tensor, g is the acceleration of gravity, fint 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:
(17)τ=μ⋅e•d,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:
(18)∂n∂t+∂nux∂x+∂nuy∂y+∂nuz∂z=0,(19)∂nux∂t+∂nux2∂x+∂nuxuy∂y+∂nuxuz∂z=−nρf∂p∂x−1ρf∂nτxx∂x+∂nτyx∂y+∂nτzx∂z+ngx+fintxρf,(20)∂nuy∂t+∂nuyux∂x+∂nuy2∂y+∂nuyuz∂z=−nρf∂p∂y−1ρf∂nτxy∂x+∂nτyy∂y+∂nτzy∂z+ngy+fintyρf,(21)∂nuz∂t+∂nuzux∂x+∂nuz2∂y+∂nuz2∂z=−nρf∂p∂z−1ρf∂nτxz∂x+∂nτyz∂y+∂nτzz∂z+ngz+fintzρf.The unit volume solid-liquid interaction force can be expressed as:(22)fintj=βintjυ¯j−uj,where the solid-liquid friction coefficient can be calculated by the following equation
(23)βintj=150⋅1−n2n⋅dp¯2uf+1.75⋅1−ndp¯⋅ρf⋅υ¯j−ujn≤0.8,34⋅ρf⋅CD⋅1−n⋅n−1.7dp¯⋅υ¯j−ujn>0.8,where βintj friction coefficient is defined as fluid, combining Eqs. (21) and (22), we can obtain that:
(24)fdij=−βintj1−nυ¯j−uj+∇pjπ6dpi3.


## 2.2.1. The Force on the Particle

Figure1 shows a fixed control volume, where ΔV=Δx×Δy×Δz, and the number of particles in the fluid unit is np.Figure 1
Fluid through the particle assemble.Assuming that the fluid flows only in thex-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–25]:
(1)fdsum=∑i=1npfdix=−fintxΔV−dpdxπ6∑i=1npdpi3,where fintx is the interaction force between the fluid and solid phases in the unit volume, dpi 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,⋯,np, j=x,y,z.Considering the porosity, it can be expressed as:(2)n=1−1ΔVπ6∑i=1npdpi3.Substituting Eq. (2) into Eq. (1) yields:
(3)fdsum=∑i=1npfdix=−fintx1−n+dpdxπ6∑i=1npdpi3.Thus, the total force acting on the particle is(4)fdij=−fintj1−n+∇pjπ6dpi3.The interaction force between the particle and fluid under fluid-solid coupling can be expressed as:(5)fintj=n⋅∇pj.Substituting Eq. (5) into Eq. (4), the total force on the particle under the action of fluid-solid coupling can be obtained as:
(6)fdij=−∇pj1−nπ6dpi3.


## 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.(7)dpdx=−ρfgKux0=−νfρfkux0,where νf,ρf is the kinematic viscosity and density of the fluid, respectively, ux0is 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:(8)k=c⋅dp¯2n31−n2,where c is the Kozeny-Carman constant, generally taken between 0.003 and 0.0055; dp¯ 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.(9)dpdx=−150⋅1−n2n3dp¯2νfρfux0+1.75⋅1−nn3dp¯ρfux02.By substitutingnux=ux0 into Eq. (9) yields:
(10)dpdx=−150⋅1−n2n2dp¯2νfρfux+1.75⋅1−nndp¯ρfux2.Since the particle forces are generated by relative fluid flow, the relative velocityurx=υ¯x−ux (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:
(11)dpdx=150⋅1−n2n2dp¯2νfρfux+1.75⋅1−nndp¯ρfυ¯x−ux⋅υ¯x−ux.Equation (11) can be expressed as:
(12)∇pj=150⋅1−n2n2dp¯2νfρfux+1.75⋅1−nndp¯ρfυ¯j−uj⋅υ¯j−uj.For large porosity (n≥0.8), the pressure gradient can be calculated by the following equation.
(13)∇pj=34ρfCd1−n⋅n−2.7d¯pυ¯j−uj⋅υ¯j−uj,where Cd is the drag coefficient of the sphere, which is a function related to the Reynolds coefficient Rep:
(14)Cd=24Rep1+0.15⋅Rep0.687Rep≤1000,0.44Rep≥1000,Rep=nυ¯j−ujd¯pυf.


## 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
(15)∂n∂t=−∇⋅nu,(16)∂nu∂t=−∇⋅nuu−nρf∇p−nρf∇⋅τ+ng+fintρf,where u is the flow velocity vector, τ is the viscous stress tensor, g is the acceleration of gravity, fint 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:
(17)τ=μ⋅e•d,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:
(18)∂n∂t+∂nux∂x+∂nuy∂y+∂nuz∂z=0,(19)∂nux∂t+∂nux2∂x+∂nuxuy∂y+∂nuxuz∂z=−nρf∂p∂x−1ρf∂nτxx∂x+∂nτyx∂y+∂nτzx∂z+ngx+fintxρf,(20)∂nuy∂t+∂nuyux∂x+∂nuy2∂y+∂nuyuz∂z=−nρf∂p∂y−1ρf∂nτxy∂x+∂nτyy∂y+∂nτzy∂z+ngy+fintyρf,(21)∂nuz∂t+∂nuzux∂x+∂nuz2∂y+∂nuz2∂z=−nρf∂p∂z−1ρf∂nτxz∂x+∂nτyz∂y+∂nτzz∂z+ngz+fintzρf.The unit volume solid-liquid interaction force can be expressed as:(22)fintj=βintjυ¯j−uj,where the solid-liquid friction coefficient can be calculated by the following equation
(23)βintj=150⋅1−n2n⋅dp¯2uf+1.75⋅1−ndp¯⋅ρf⋅υ¯j−ujn≤0.8,34⋅ρf⋅CD⋅1−n⋅n−1.7dp¯⋅υ¯j−ujn>0.8,where βintj friction coefficient is defined as fluid, combining Eqs. (21) and (22), we can obtain that:
(24)fdij=−βintj1−nυ¯j−uj+∇pjπ6dpi3.


## 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 Table1. 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, 100×100×200mm, 100×100×250mm, and 100×100×300mm, as shown in Table 2.Table 1
Mesoscopic parameters of numerical simulation materials.
KindsParticle size (mm)Density (kg·m-3)Normal stiffness (N·m-1)Tangential stiffness (N·m-1)Coefficient of frictionViscosity (MPa·s)Ball0.1 ~ 2.023001.0×1041.0×1040.5Fluid1000Functions related to runtimeWall1.0×1071.0×1070.3Table 2
Design of numerical experiment.
Water pressure (MPa)Crack width (mm)Number particlesModel size (mm)0.58116586100×100×200155448100×100×250174879100×100×3000.88116586100×100×200155448100×100×250174879100×100×30018116586100×100×200155448100×100×250174879100×100×3001.28116586100×100×200155448100×100×250174879100×100×300


#### 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 is0.3×0.1×0.1m (Length×Width×Height), as shown in Figure 2.Figure 2
A hollow cuboid.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 radiusR adopts from Rmin to Rmax uniform distribution, and the average radius of particles is R¯=Rmin+Rmax/2. Considering the calculation accuracy and calculation time, the minimum radius of particles Rmin=0.05mm, the maximum radius Rmax=1mm, and the number N of particles is determined by the porosity n.
(25)N=1−n⋅V4/3⋅π⋅R¯3,where n is the porosity; V is the cuboid volume, m3; 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 Figure3.Figure 3
Particle assemble.


#### 3.1.2. Generate Fractures

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


#### 3.1.3. Set up Fluid Unit

A total of 250 fluid units are set in this model, and the unit size is2×1×1mm. 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−4s, as shown in Figure 5.Figure 5
Schematic diagram of fluid unit model.


### 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 Table1. 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.Figure 6
Model diagram of a moment of particle equilibrium disruption.


### 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 Figure7. The conditions for convergence are set as follows.
(26)ratiosum=sumofaverageunbalancedforcesofallparticlessumofaveragecontactforcesofallparticles,ratiomax=maximumunbalancedforceofparticlesmaximumcontactforceofparticles,Figure 7
Schematic diagram of PFC iterative process.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 (Figure8). 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, 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.Figure 8
The vortex phenomenon and particle movement trend in the initial stage.Figure 9
The change of the particle flow direction.Figure 10
The overall increase in the process of particle flow.According to Figure11, it can be found that the seepage flow of particles under different water pressure all increase with time, such as model size 100×100×300mm; with the increase of pressure gradient, the percolation volume of particles after 0.45 s are 4.4×10−4m3, 6.9×10−4m3, 9.2×10−4m3, and 1.2×10−3m3, 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×250mm; and the percolation volume of particles after 0.35 s was 4.7×10−4m3, 6.9×10−4m3, 8.5×10−4m3, and 1.4×10−3m3 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−4m3, 6.9×10−4m3, 9.2×10−4m3, and 1.3×10−3m3, 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.Figure 11
The relationship curve between seepage flow of particles and time under different water pressures for particle size of100×100×300mm at (a) P=0.5MPa, (b) P=0.8MPa, (c) P=1.0MPa, and (d) P=1.2MPa; and particle size of 100×100×250mm at (e) P=0.5MPa, (f) P=0.8MPa, (g) P=1.0MPa, and (h) P=1.2MPa; and particle size of 100×100×200mm at (i) P=0.5MPa, (j) P=0.8MPa, (k) P=1.0MPa, and (l) P=1.2MPa.
(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 Figure12.Figure 12
Simulated fitting curves for particle size is100×100×300mm (a) P=0.5MPa, (b) P=0.8MPa, (c) P=1.0MPa, and (d) P=1.2MPa; and particle size is 100×100×250mm (e) P=0.5MPa, (f) P=0.8MPa, (g) P=1.0MPa, and (h) P=1.2MPa; and particle size is 100×100×200mm (i) P=0.5MPa, (j) P=0.8MPa, (k) P=1.0MPa, and (l) P=1.2MPa.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure12 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 coefficientR2is 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 coefficientR2is 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 Figure13 shows the fitted curves of pressure gradient and particle speed.Figure 13
The relationship between pressure gradient and the particle speed.Figure13 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 R2 equal to 0.9312, and the fitting formula is as follows:
(27)v=1.328×10−2−1.156×10−21+expΔp−4.42×106/3.54×105.


## 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 Table1. 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, 100×100×200mm, 100×100×250mm, and 100×100×300mm, as shown in Table 2.Table 1
Mesoscopic parameters of numerical simulation materials.
KindsParticle size (mm)Density (kg·m-3)Normal stiffness (N·m-1)Tangential stiffness (N·m-1)Coefficient of frictionViscosity (MPa·s)Ball0.1 ~ 2.023001.0×1041.0×1040.5Fluid1000Functions related to runtimeWall1.0×1071.0×1070.3Table 2
Design of numerical experiment.
Water pressure (MPa)Crack width (mm)Number particlesModel size (mm)0.58116586100×100×200155448100×100×250174879100×100×3000.88116586100×100×200155448100×100×250174879100×100×30018116586100×100×200155448100×100×250174879100×100×3001.28116586100×100×200155448100×100×250174879100×100×300


### 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 is0.3×0.1×0.1m (Length×Width×Height), as shown in Figure 2.Figure 2
A hollow cuboid.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 radiusR adopts from Rmin to Rmax uniform distribution, and the average radius of particles is R¯=Rmin+Rmax/2. Considering the calculation accuracy and calculation time, the minimum radius of particles Rmin=0.05mm, the maximum radius Rmax=1mm, and the number N of particles is determined by the porosity n.
(25)N=1−n⋅V4/3⋅π⋅R¯3,where n is the porosity; V is the cuboid volume, m3; 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 Figure3.Figure 3
Particle assemble.


### 3.1.2. Generate Fractures

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


### 3.1.3. Set up Fluid Unit

A total of 250 fluid units are set in this model, and the unit size is2×1×1mm. 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−4s, as shown in Figure 5.Figure 5
Schematic diagram of fluid unit model.


## 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 is0.3×0.1×0.1m (Length×Width×Height), as shown in Figure 2.Figure 2
A hollow cuboid.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 radiusR adopts from Rmin to Rmax uniform distribution, and the average radius of particles is R¯=Rmin+Rmax/2. Considering the calculation accuracy and calculation time, the minimum radius of particles Rmin=0.05mm, the maximum radius Rmax=1mm, and the number N of particles is determined by the porosity n.
(25)N=1−n⋅V4/3⋅π⋅R¯3,where n is the porosity; V is the cuboid volume, m3; 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 Figure3.Figure 3
Particle assemble.


## 3.1.2. Generate Fractures

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


## 3.1.3. Set up Fluid Unit

A total of 250 fluid units are set in this model, and the unit size is2×1×1mm. 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−4s, as shown in Figure 5.Figure 5
Schematic diagram of fluid unit model.


## 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 Table1. 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.Figure 6
Model diagram of a moment of particle equilibrium disruption.


## 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 Figure7. The conditions for convergence are set as follows.
(26)ratiosum=sumofaverageunbalancedforcesofallparticlessumofaveragecontactforcesofallparticles,ratiomax=maximumunbalancedforceofparticlesmaximumcontactforceofparticles,Figure 7
Schematic diagram of PFC iterative process.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 (Figure8). 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, 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.Figure 8
The vortex phenomenon and particle movement trend in the initial stage.Figure 9
The change of the particle flow direction.Figure 10
The overall increase in the process of particle flow.According to Figure11, it can be found that the seepage flow of particles under different water pressure all increase with time, such as model size 100×100×300mm; with the increase of pressure gradient, the percolation volume of particles after 0.45 s are 4.4×10−4m3, 6.9×10−4m3, 9.2×10−4m3, and 1.2×10−3m3, 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×250mm; and the percolation volume of particles after 0.35 s was 4.7×10−4m3, 6.9×10−4m3, 8.5×10−4m3, and 1.4×10−3m3 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−4m3, 6.9×10−4m3, 9.2×10−4m3, and 1.3×10−3m3, 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.Figure 11
The relationship curve between seepage flow of particles and time under different water pressures for particle size of100×100×300mm at (a) P=0.5MPa, (b) P=0.8MPa, (c) P=1.0MPa, and (d) P=1.2MPa; and particle size of 100×100×250mm at (e) P=0.5MPa, (f) P=0.8MPa, (g) P=1.0MPa, and (h) P=1.2MPa; and particle size of 100×100×200mm at (i) P=0.5MPa, (j) P=0.8MPa, (k) P=1.0MPa, and (l) P=1.2MPa.
(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 Figure12.Figure 12
Simulated fitting curves for particle size is100×100×300mm (a) P=0.5MPa, (b) P=0.8MPa, (c) P=1.0MPa, and (d) P=1.2MPa; and particle size is 100×100×250mm (e) P=0.5MPa, (f) P=0.8MPa, (g) P=1.0MPa, and (h) P=1.2MPa; and particle size is 100×100×200mm (i) P=0.5MPa, (j) P=0.8MPa, (k) P=1.0MPa, and (l) P=1.2MPa.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure12 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 coefficientR2is 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 coefficientR2is 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 Figure13 shows the fitted curves of pressure gradient and particle speed.Figure 13
The relationship between pressure gradient and the particle speed.Figure13 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 R2 equal to 0.9312, and the fitting formula is as follows:
(27)v=1.328×10−2−1.156×10−21+expΔp−4.42×106/3.54×105.


## 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 (Figure8). 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, 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.Figure 8
The vortex phenomenon and particle movement trend in the initial stage.Figure 9
The change of the particle flow direction.Figure 10
The overall increase in the process of particle flow.According to Figure11, it can be found that the seepage flow of particles under different water pressure all increase with time, such as model size 100×100×300mm; with the increase of pressure gradient, the percolation volume of particles after 0.45 s are 4.4×10−4m3, 6.9×10−4m3, 9.2×10−4m3, and 1.2×10−3m3, 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×250mm; and the percolation volume of particles after 0.35 s was 4.7×10−4m3, 6.9×10−4m3, 8.5×10−4m3, and 1.4×10−3m3 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−4m3, 6.9×10−4m3, 9.2×10−4m3, and 1.3×10−3m3, 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.Figure 11
The relationship curve between seepage flow of particles and time under different water pressures for particle size of100×100×300mm at (a) P=0.5MPa, (b) P=0.8MPa, (c) P=1.0MPa, and (d) P=1.2MPa; and particle size of 100×100×250mm at (e) P=0.5MPa, (f) P=0.8MPa, (g) P=1.0MPa, and (h) P=1.2MPa; and particle size of 100×100×200mm at (i) P=0.5MPa, (j) P=0.8MPa, (k) P=1.0MPa, and (l) P=1.2MPa.
(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 Figure12.Figure 12
Simulated fitting curves for particle size is100×100×300mm (a) P=0.5MPa, (b) P=0.8MPa, (c) P=1.0MPa, and (d) P=1.2MPa; and particle size is 100×100×250mm (e) P=0.5MPa, (f) P=0.8MPa, (g) P=1.0MPa, and (h) P=1.2MPa; and particle size is 100×100×200mm (i) P=0.5MPa, (j) P=0.8MPa, (k) P=1.0MPa, and (l) P=1.2MPa.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure12 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 coefficientR2is 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 coefficientR2is 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 Figure13 shows the fitted curves of pressure gradient and particle speed.Figure 13
The relationship between pressure gradient and the particle speed.Figure13 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 R2 equal to 0.9312, and the fitting formula is as follows:
(27)v=1.328×10−2−1.156×10−21+expΔp−4.42×106/3.54×105.


## 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



---

*Source: 5599748-2021-04-09.xml*