# Study on Characteristics of the Water-Sediment Two-Phase Flow in Fractures Based on Numerical Simulation

Authors: Lili Cao; Chuan Wu; Hanhua Xu; Hao Yu; Ming Li; Yu Han; Zhanqing Chen
Journal: Geofluids
 (2022)

---

## Abstract

The instability of water-sediment flow in fractures can easily induce water-sediment disasters. Therefore, it is of great significance for the prevention and control of water-sediment inrush to study the water-sediment two-phase flow in fractures. Based on the water-sediment two-phase flow theory, a model of the water-sediment two-phase flow system was established. The Ansys Fluent software was used to study the characteristics of the water-sediment two-phase flow in smooth and rough fractures. The spatial-temporal evolution laws of the water-sediment two-phase flow were studied; the results indicated that they did not change with time in the smooth fractured flow fields, while changing continuously with time in the rough fractured flow fields and in a dynamic steady state. The research results can provide references for the water-sediment two-phase flow in fractures and rock mass.


---

## Body

## 1. Introduction

Water and sand inrush is one of the main disasters in the mining of shallow coal seams in Western China [1]. Fractures or small faults often connect the water-rich layer with unbalanced sand bodies due to weathering [2, 3]. Then, under the action of gravity, the water-sediment mixture will submerge the equipment underground and even cause casualties. Water-sediment two-phase flow is one of the important incentives of water-sediment inrush disasters [4]. Therefore, it is of great significance to study the water-sediment two-phase flow in fractures to ensure the safe and efficient production of mines [5, 6].Scholars have carried out numerous researches on water-sediment two-phase mixture seepage characteristics and their influencing factors [7, 8]. Liu et al. analyzed the variation laws of permeability parameters of rock-fractured rock mass fracture roughness, sand particle size, and sediment concentration [9]. Yang et al. revealed the influence of sediment-water interactions on movement characteristics and rheological characteristics of water-sediment mixture [10]. Qi and Bo studied the failure mechanism and evolution characteristics of water-sediment inrush disasters caused by the instability of the filling medium of karst caves [11]. Yang et al. conducted a nonlinear flow experiment on the fracture network and revealed the mechanism of water inrush in fractures [12]. Zhang et al. analyzed the influence of the aeolian sand’s particle size distribution, void ratio, and initial mass on the flow characteristics [13]. Xu et al. simulated the process of water-sediment inrush in the goaf and proposed the prevention and control techniques [14]. However, due to the limitation of test conditions, there were still some shortcomings in the study of the flow characteristics of the water-sand mixture by experimental means. For example, it was difficult to find the test equipment that met the requirements of high performance when the water-sand mixture needed to be transported at a stable flow rate. In addition, it was difficult to describe the vortex structure efficiently and accurately using the existing experimental means, which made it difficult to conduct a more in-depth study on the flow field of water-sediment two-phase flow.Considering the complex geological conditions of undergrounding mining [15, 16], some scholars have used the finite element analysis and discrete element method to reveal the whole process of water-sediment two-phase flow and have achieved remarkable results [17]. Du et al. used the FORTRAN programming language for numerical calculation and determined the main factors affecting the characteristic parameters of water-sediment two-phase mixture seepage in fractured rock mass [18]. Pu analyzed the laws of water-sediment two-phase flow in fractures using the lattice Boltzmann method. They also simulated the water-sand inrush by numerical simulation and analyzed the impacts of particle size and fracture width on the water-sand inrush velocity [19]. Guo et al. simulated the whole process of water-sand inrush process during the mining of the working face by using the self-developed simulation test system [20]. Lei et al. used the PFC3D to simulate the development process of overburden fissures and the water hydraulic of Paleogene aquifers [21]. In summary, scholars have conducted studies on the characteristics of water-sediment two-phase seepage in fractures using experimental methods, but there are still some deficiencies. The research only focuses on the seepage of a single fluid and rarely involves the liquid-solid two-phase movement in fractures. There are few studies on seepage of the particle phase and continuous phase simultaneously in fractures. Meanwhile, there is currently no universally accepted standard for the description of fracture surface morphology [2, 22].In this paper, Ansys Fluent 17.0 will be used to simulate water-sediment two-phase flow in fractures with smooth and rough surfaces. Firstly, the mechanical model and numerical model of the water-sediment two-phase flow will be established. Then, different turbulence simulation methods will be used to calculate the flow field of the continuous phase. The calculated results will be compared with the experimental results, and a preferable method will be selected for fractures with rough and complex surfaces. Finally, the dispersed phase particles will be injected into the flow field of the continuous phase in dynamic equilibrium, and the coupling calculation will be performed. The calculated results will be compared with the experimental results to further study the influence factors of the water-sediment seepage field. The research results are aimed at providing a theoretical foundation for the mechanism of water-sediment inrush during coal mining.


## 2. The Mechanical Model of Water-Sediment Two-Phase Flow in Fractures

### 2.1. The Model Hypothesis and Selection of Water-Sediment Flow in Fractures

The following assumptions were made during the numerical simulation of water-sediment flow in fractures [23, 24]. (I) Water is incompressible. That is, the density of water was a constant. (II) Sand particles were spheres with the same particle sizes. (III) The rigid entities were spherical sand particle units and would not obtain obvious damage. (IV) The water-sediment mixture flowed into fractures from the flow guide transition plate with a larger aperture. The transition was smooth at the inlet of the fracture, and the flow velocity was uniform on the cross-section of the inlet of the fracture. (V) The disperse phase sand particle velocity was the same with that of the continuous-phase fluid at the inlet of the fracture. (VI) When sand particles escaped from the exit of the fracture, the tracking would be stopped. (VII) The physical quantity did not change along the X3 direction, and the flow was taken as a two-dimensional flow.The model of water-sediment flow in fractures was composed of the continuous-phase model and the discrete-phase model. The Reynolds averaging method and the large eddy simulation method were used to simulate the turbulent flow for the continuous phase [25]. According to different methods of stress treatment, the Reynolds averaging model could be divided into the eddy viscosity model and Reynolds stress model [26]. The former is to introduce the turbulent viscosity to treat the Reynolds stress instead of treating the stress directly. The Reynolds stress is expressed by a function of the turbulent viscosity. The commonly used eddy viscosity models in engineering include the Spalart-Allmaras model, k-ε model, and k-ω model. In this paper, the selected turbulence models include the RNG k-ε model in the k-ε model, Realizable k-ε model, and SST k-ω model in the k-ω model. In the Reynolds stress model (RSM), the equations of the Reynolds stress are established directly. The isotropic viscosity hypotheses can be avoided. Meanwhile, the influences of fluid rotation, streamline bending, and sharp change in the strain rate will be considered. Compared with the eddy viscosity models, RSM is more suitable for accurate prediction of complex flows. The Stress-Omega RSM was adopted in this paper. The Stress-Omega model of the wall-modeled large eddy simulation (WMLES) was used for the mathematical model of subgrid-scale stress [27].It is necessary to simulate the rotation and translation of sand particles in the study of water-sediment flow in fractures. In this paper, the random orbit model was used to simulate the particle diffusion caused by the turbulent motion of the continuous phase. Meanwhile, based on the distinct element method and the program BALL proposed by Cundall and Strack [28], the Hertzian-Dashpot model and the Rolling Friction model were used to simulate the normal contact and the tangential contact force among particles, respectively.


### 2.2. Computational Domain of the Water-Sediment Flow in Fractures

The computational domainΩ was a set of particles surrounded by two fracture surfaces, the inlet section and the outlet section. The set of particles was time varying. That is, water and sand entered continuously from the inlet section and flowed out from the outlet section, and the control volume was constant. The computational domains of the smooth fractures and the rough fractures were descried as follows.


#### 2.2.1. The Computational Domain of the Smooth Fractures

As shown in Figure1, water and sand flowed into the area formed by two parallel smooth fracture surfaces. The distance between the two fracture surfaces is h, and the length of the fracture is L. The projection of the upper and lower fracture surfaces on the OX1X2 section is a straight line. The boundary of the computational domain is the inlet section OD1¯, the upper fracture section D1D2¯, the outlet section D2D3¯, and the lower fracture surface D3O¯.Figure 1
The computational domain of the smooth fracture.


#### 2.2.2. The Computational Domain of the Rough Fractures

As shown in Figure2, water and sand flowed in the area formed by two anastomotic fracture surfaces. The curve of the upper and lower fracture surfaces on the cross-section OX1X2 is a broken line which divides the surface into 50 equal broken line segments marked by FLiu and FLid, i=1,2,3,⋯50. The computational domain Ω is formed by the inlet section OC1¯, the upper fracture surface ∪i=150FLiu, the outlet section C2C3¯, and the lower fracture section ∪i=150FLid.Figure 2
The computational domain of the rough fracture.


### 2.3. Boundary Conditions and the Initial Conditions

The boundary and initial conditions needed to be specified to solve the equations for the water-sediment two-phase flow fields. The boundary and initial conditions of the continuous phase and the discrete phase were given, respectively.


#### 2.3.1. Boundary Conditions of the Continuous Phase

The inlet boundary was set as the velocity inlet. The four-level seepage velocity was used in the experiment. The hydraulic diameter was twice the aperture of the fracture inlet. The turbulence intensity could be obtained by the following equation:(1)I=0.16ReDH−1/8,whereReDH is the Reynolds number with the hydraulic diameter of DH. The outlet boundary was set as the pressure outlet. The outlet connected with the air, so the outlet pressure was the standard atmospheric pressure of 1. The wall boundary was set as a no-slip wall boundary. The gravity direction was along the negative direction of X2.


#### 2.3.2. Boundary Conditions of the Dispersed Phase

The discrete-phase sand particle velocity was the same as the continuous-phase fluid velocity at the inlet of the fracture. The volume concentration of sand particles was set to 4.06%. The outlet of the fracture was set as the escape boundary of the discrete phase, where the tracing would stop. The sand particles and the wall boundary were set as the rebound boundary. The values of the initial boundary conditions were fixed and did not change with time.


### 2.4. Numerical Calculation Methods and Parameter Settings

The grid division is a key step in the numerical calculation by using the finite volume method. Its quality directly affects the accuracy of the numerical calculation results. In this paper, the structured grids were used, because the boundary of the smooth fracture model was relatively regular. The grid nodes were evenly distributed at the inlet to make the sand particles evenly distributed along theX2 direction. The Stress-Omega RSM was used in this paper. y+ was the dimensionless distance between the first layer of grid nodes and the wall, which was approximately 1. Therefore, the fine grids were arranged in the near-wall region. After calculation, it was necessary to check whether y+ could meet the requirements.In order to ensure the accuracy of the calculation results, the grid independence test was carried out. The size of the initial grid in the boundary layer was 0.005 mm, and the global grid size was 0.08 mm. After the grid independence test, the size of the grid in the boundary layer was 0.003 mm, and the global grid size was 0.02 mm. Through calculation, it was found that both grids could satisfy thaty+ was less than 1. The inlet pressure difference of two grids was less than 1%. The velocity difference was less than 0.3% on the cross-section of X1=20mm, indicating that the initial grids could meet the requirements of the calculation accuracy. In view of the calculation accuracy and efficiency, the initial grids were selected for calculation. Figure 3 shows the grid division.Figure 3
The grid of the smooth fracture.The hybrid grid was used, due to sharp bending of the wall of the rough fracture model. Multilayer structured grids were set in the boundary layer, and the remaining computational domain was unstructured quadrilateral grids. They+ was tested, and the grid independence test was conducted. The total number of grids was about 145,000, as shown in Figure 4. The grids in the boundary layer are shown in Figure 5.Figure 4
Grids of the rough fracture.Figure 5
Grids in the boundary layer of the rough fracture.After the grid division of the computational domain, it was necessary to discretize the governing equations and definite conditions (boundary conditions and initial conditions) in the space domain and the time domain. The Least Squares Cell-Based method was used. The pressure interpolation method was the pressure staggering option (PRESTO). The momentum conservation equation, the equation for the turbulence kinetic energyk, and the equation for the specific dissipation rate ω adopted the third-order MUSCL scheme to reduce the numerical diffusion. The Reynolds stress equation used the QUICK discrete scheme (the quadrilateral grids used the QUICK discrete scheme, and other grids used the second-order discrete scheme). The second-order implicit time integration scheme was used to solve the continuous-phase equations, and the Runge-Kutta method was used to calculate the discrete-phase equations.Table1 shows the material properties. Table 2 shows the interaction parameters of sand particles and the wall. Table 3 shows the parameters among sand particles.Table 1
Material properties [29].
WaterSand particlesDensity (kg·m-3)998.22650Dynamic viscosity (kg·m-1·s-1)1.003×10−3/Elastic modulus (GPa)/55.9Poisson’s ratio/0.13Particle size (mm)/0.04Table 2
Interaction parameters between the sand particles and the wall.
PropertiesBetween sand particles and wallFriction factor0.45Normal recovery factor0.2Tangential recovery factor0.9Table 3
Interaction parameters among sand particles [29].
PropertiesAmong sand particlesCoefficient of static friction0.3Coefficient of sliding friction0.2Coefficient of restitution0.05The pressure-based solver was used for calculation. The pressure-velocity coupling PSIO (pressure implicit with splitting of operator) algorithm was used. The interphase coupling calculation method was used, and interactions between the continuous phase and the discrete phase were considered. After adjustment, the time step was set as 10-5 s. The tracking time step of the sand particles was one percent of the continuous-phase calculation time step. The velocity and inlet pressure changes were monitored during the calculation process. Parallel computing was used to enhance the calculation speed.


### 2.5. The Experimental Scheme

In this paper, some water-sediment seepage experimental data were from the results of Zhanqing Chen’s Research Group, State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology. Restricted by the test conditions, only the water-sediment two-phase flow in regular and rough fractures was studied without the test of water-sediment two-phase flow in smooth fractures.Numerical simulation is a necessary supplement and extension of experimental research and theoretical analysis. It can simulate complex flow problems and obtain the data that is difficult to obtain in experiments [30]. In this paper, the numerical simulation was used to study the characteristics of the water-sediment two-phase flow in fractures.To better describe the flows in fractures, it was necessary to compare the simulation results and the experimental data of the unidirectional flow in fractures before performing the simulation of two-phase flow in fractures. In this paper, several different turbulence simulation calculation methods were used to calculate the continuous-phase flow fields, and the calculation results were compared with the experimental results. Based on the comparison results, the appropriate turbulence simulation method was chosen.The Ansys Fluent was used to simulate the water-sediment two-phase flow in the smooth and rough fractures. The geometric description of the fracture surface was given before the simulation. Then, the grid division was performed in the computational domain. The experimental parameters were given. The sand densityρp was 2650 kg/m3, the sand particle size Dp was 0.04 mm, and the volume concentration Ф of the sand particles was 4.06%. On this basis, the numerical model was established. Finally, the inlet velocities of the fractures were set as 0.349 m/s, 0.532 m/s, 0.697 m/s, and 0.869 m/s. Then, the program was debugged.After the simulation, the spatial-temporal evolution laws and influence factors of the water-sediment two-phase flow were analyzed according to the simulation results.


## 2.1. The Model Hypothesis and Selection of Water-Sediment Flow in Fractures

The following assumptions were made during the numerical simulation of water-sediment flow in fractures [23, 24]. (I) Water is incompressible. That is, the density of water was a constant. (II) Sand particles were spheres with the same particle sizes. (III) The rigid entities were spherical sand particle units and would not obtain obvious damage. (IV) The water-sediment mixture flowed into fractures from the flow guide transition plate with a larger aperture. The transition was smooth at the inlet of the fracture, and the flow velocity was uniform on the cross-section of the inlet of the fracture. (V) The disperse phase sand particle velocity was the same with that of the continuous-phase fluid at the inlet of the fracture. (VI) When sand particles escaped from the exit of the fracture, the tracking would be stopped. (VII) The physical quantity did not change along the X3 direction, and the flow was taken as a two-dimensional flow.The model of water-sediment flow in fractures was composed of the continuous-phase model and the discrete-phase model. The Reynolds averaging method and the large eddy simulation method were used to simulate the turbulent flow for the continuous phase [25]. According to different methods of stress treatment, the Reynolds averaging model could be divided into the eddy viscosity model and Reynolds stress model [26]. The former is to introduce the turbulent viscosity to treat the Reynolds stress instead of treating the stress directly. The Reynolds stress is expressed by a function of the turbulent viscosity. The commonly used eddy viscosity models in engineering include the Spalart-Allmaras model, k-ε model, and k-ω model. In this paper, the selected turbulence models include the RNG k-ε model in the k-ε model, Realizable k-ε model, and SST k-ω model in the k-ω model. In the Reynolds stress model (RSM), the equations of the Reynolds stress are established directly. The isotropic viscosity hypotheses can be avoided. Meanwhile, the influences of fluid rotation, streamline bending, and sharp change in the strain rate will be considered. Compared with the eddy viscosity models, RSM is more suitable for accurate prediction of complex flows. The Stress-Omega RSM was adopted in this paper. The Stress-Omega model of the wall-modeled large eddy simulation (WMLES) was used for the mathematical model of subgrid-scale stress [27].It is necessary to simulate the rotation and translation of sand particles in the study of water-sediment flow in fractures. In this paper, the random orbit model was used to simulate the particle diffusion caused by the turbulent motion of the continuous phase. Meanwhile, based on the distinct element method and the program BALL proposed by Cundall and Strack [28], the Hertzian-Dashpot model and the Rolling Friction model were used to simulate the normal contact and the tangential contact force among particles, respectively.


## 2.2. Computational Domain of the Water-Sediment Flow in Fractures

The computational domainΩ was a set of particles surrounded by two fracture surfaces, the inlet section and the outlet section. The set of particles was time varying. That is, water and sand entered continuously from the inlet section and flowed out from the outlet section, and the control volume was constant. The computational domains of the smooth fractures and the rough fractures were descried as follows.


### 2.2.1. The Computational Domain of the Smooth Fractures

As shown in Figure1, water and sand flowed into the area formed by two parallel smooth fracture surfaces. The distance between the two fracture surfaces is h, and the length of the fracture is L. The projection of the upper and lower fracture surfaces on the OX1X2 section is a straight line. The boundary of the computational domain is the inlet section OD1¯, the upper fracture section D1D2¯, the outlet section D2D3¯, and the lower fracture surface D3O¯.Figure 1
The computational domain of the smooth fracture.


### 2.2.2. The Computational Domain of the Rough Fractures

As shown in Figure2, water and sand flowed in the area formed by two anastomotic fracture surfaces. The curve of the upper and lower fracture surfaces on the cross-section OX1X2 is a broken line which divides the surface into 50 equal broken line segments marked by FLiu and FLid, i=1,2,3,⋯50. The computational domain Ω is formed by the inlet section OC1¯, the upper fracture surface ∪i=150FLiu, the outlet section C2C3¯, and the lower fracture section ∪i=150FLid.Figure 2
The computational domain of the rough fracture.


## 2.2.1. The Computational Domain of the Smooth Fractures

As shown in Figure1, water and sand flowed into the area formed by two parallel smooth fracture surfaces. The distance between the two fracture surfaces is h, and the length of the fracture is L. The projection of the upper and lower fracture surfaces on the OX1X2 section is a straight line. The boundary of the computational domain is the inlet section OD1¯, the upper fracture section D1D2¯, the outlet section D2D3¯, and the lower fracture surface D3O¯.Figure 1
The computational domain of the smooth fracture.


## 2.2.2. The Computational Domain of the Rough Fractures

As shown in Figure2, water and sand flowed in the area formed by two anastomotic fracture surfaces. The curve of the upper and lower fracture surfaces on the cross-section OX1X2 is a broken line which divides the surface into 50 equal broken line segments marked by FLiu and FLid, i=1,2,3,⋯50. The computational domain Ω is formed by the inlet section OC1¯, the upper fracture surface ∪i=150FLiu, the outlet section C2C3¯, and the lower fracture section ∪i=150FLid.Figure 2
The computational domain of the rough fracture.


## 2.3. Boundary Conditions and the Initial Conditions

The boundary and initial conditions needed to be specified to solve the equations for the water-sediment two-phase flow fields. The boundary and initial conditions of the continuous phase and the discrete phase were given, respectively.


### 2.3.1. Boundary Conditions of the Continuous Phase

The inlet boundary was set as the velocity inlet. The four-level seepage velocity was used in the experiment. The hydraulic diameter was twice the aperture of the fracture inlet. The turbulence intensity could be obtained by the following equation:(1)I=0.16ReDH−1/8,whereReDH is the Reynolds number with the hydraulic diameter of DH. The outlet boundary was set as the pressure outlet. The outlet connected with the air, so the outlet pressure was the standard atmospheric pressure of 1. The wall boundary was set as a no-slip wall boundary. The gravity direction was along the negative direction of X2.


### 2.3.2. Boundary Conditions of the Dispersed Phase

The discrete-phase sand particle velocity was the same as the continuous-phase fluid velocity at the inlet of the fracture. The volume concentration of sand particles was set to 4.06%. The outlet of the fracture was set as the escape boundary of the discrete phase, where the tracing would stop. The sand particles and the wall boundary were set as the rebound boundary. The values of the initial boundary conditions were fixed and did not change with time.


## 2.3.1. Boundary Conditions of the Continuous Phase

The inlet boundary was set as the velocity inlet. The four-level seepage velocity was used in the experiment. The hydraulic diameter was twice the aperture of the fracture inlet. The turbulence intensity could be obtained by the following equation:(1)I=0.16ReDH−1/8,whereReDH is the Reynolds number with the hydraulic diameter of DH. The outlet boundary was set as the pressure outlet. The outlet connected with the air, so the outlet pressure was the standard atmospheric pressure of 1. The wall boundary was set as a no-slip wall boundary. The gravity direction was along the negative direction of X2.


## 2.3.2. Boundary Conditions of the Dispersed Phase

The discrete-phase sand particle velocity was the same as the continuous-phase fluid velocity at the inlet of the fracture. The volume concentration of sand particles was set to 4.06%. The outlet of the fracture was set as the escape boundary of the discrete phase, where the tracing would stop. The sand particles and the wall boundary were set as the rebound boundary. The values of the initial boundary conditions were fixed and did not change with time.


## 2.4. Numerical Calculation Methods and Parameter Settings

The grid division is a key step in the numerical calculation by using the finite volume method. Its quality directly affects the accuracy of the numerical calculation results. In this paper, the structured grids were used, because the boundary of the smooth fracture model was relatively regular. The grid nodes were evenly distributed at the inlet to make the sand particles evenly distributed along theX2 direction. The Stress-Omega RSM was used in this paper. y+ was the dimensionless distance between the first layer of grid nodes and the wall, which was approximately 1. Therefore, the fine grids were arranged in the near-wall region. After calculation, it was necessary to check whether y+ could meet the requirements.In order to ensure the accuracy of the calculation results, the grid independence test was carried out. The size of the initial grid in the boundary layer was 0.005 mm, and the global grid size was 0.08 mm. After the grid independence test, the size of the grid in the boundary layer was 0.003 mm, and the global grid size was 0.02 mm. Through calculation, it was found that both grids could satisfy thaty+ was less than 1. The inlet pressure difference of two grids was less than 1%. The velocity difference was less than 0.3% on the cross-section of X1=20mm, indicating that the initial grids could meet the requirements of the calculation accuracy. In view of the calculation accuracy and efficiency, the initial grids were selected for calculation. Figure 3 shows the grid division.Figure 3
The grid of the smooth fracture.The hybrid grid was used, due to sharp bending of the wall of the rough fracture model. Multilayer structured grids were set in the boundary layer, and the remaining computational domain was unstructured quadrilateral grids. They+ was tested, and the grid independence test was conducted. The total number of grids was about 145,000, as shown in Figure 4. The grids in the boundary layer are shown in Figure 5.Figure 4
Grids of the rough fracture.Figure 5
Grids in the boundary layer of the rough fracture.After the grid division of the computational domain, it was necessary to discretize the governing equations and definite conditions (boundary conditions and initial conditions) in the space domain and the time domain. The Least Squares Cell-Based method was used. The pressure interpolation method was the pressure staggering option (PRESTO). The momentum conservation equation, the equation for the turbulence kinetic energyk, and the equation for the specific dissipation rate ω adopted the third-order MUSCL scheme to reduce the numerical diffusion. The Reynolds stress equation used the QUICK discrete scheme (the quadrilateral grids used the QUICK discrete scheme, and other grids used the second-order discrete scheme). The second-order implicit time integration scheme was used to solve the continuous-phase equations, and the Runge-Kutta method was used to calculate the discrete-phase equations.Table1 shows the material properties. Table 2 shows the interaction parameters of sand particles and the wall. Table 3 shows the parameters among sand particles.Table 1
Material properties [29].
WaterSand particlesDensity (kg·m-3)998.22650Dynamic viscosity (kg·m-1·s-1)1.003×10−3/Elastic modulus (GPa)/55.9Poisson’s ratio/0.13Particle size (mm)/0.04Table 2
Interaction parameters between the sand particles and the wall.
PropertiesBetween sand particles and wallFriction factor0.45Normal recovery factor0.2Tangential recovery factor0.9Table 3
Interaction parameters among sand particles [29].
PropertiesAmong sand particlesCoefficient of static friction0.3Coefficient of sliding friction0.2Coefficient of restitution0.05The pressure-based solver was used for calculation. The pressure-velocity coupling PSIO (pressure implicit with splitting of operator) algorithm was used. The interphase coupling calculation method was used, and interactions between the continuous phase and the discrete phase were considered. After adjustment, the time step was set as 10-5 s. The tracking time step of the sand particles was one percent of the continuous-phase calculation time step. The velocity and inlet pressure changes were monitored during the calculation process. Parallel computing was used to enhance the calculation speed.


## 2.5. The Experimental Scheme

In this paper, some water-sediment seepage experimental data were from the results of Zhanqing Chen’s Research Group, State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology. Restricted by the test conditions, only the water-sediment two-phase flow in regular and rough fractures was studied without the test of water-sediment two-phase flow in smooth fractures.Numerical simulation is a necessary supplement and extension of experimental research and theoretical analysis. It can simulate complex flow problems and obtain the data that is difficult to obtain in experiments [30]. In this paper, the numerical simulation was used to study the characteristics of the water-sediment two-phase flow in fractures.To better describe the flows in fractures, it was necessary to compare the simulation results and the experimental data of the unidirectional flow in fractures before performing the simulation of two-phase flow in fractures. In this paper, several different turbulence simulation calculation methods were used to calculate the continuous-phase flow fields, and the calculation results were compared with the experimental results. Based on the comparison results, the appropriate turbulence simulation method was chosen.The Ansys Fluent was used to simulate the water-sediment two-phase flow in the smooth and rough fractures. The geometric description of the fracture surface was given before the simulation. Then, the grid division was performed in the computational domain. The experimental parameters were given. The sand densityρp was 2650 kg/m3, the sand particle size Dp was 0.04 mm, and the volume concentration Ф of the sand particles was 4.06%. On this basis, the numerical model was established. Finally, the inlet velocities of the fractures were set as 0.349 m/s, 0.532 m/s, 0.697 m/s, and 0.869 m/s. Then, the program was debugged.After the simulation, the spatial-temporal evolution laws and influence factors of the water-sediment two-phase flow were analyzed according to the simulation results.


## 3. Results of the Numerical Simulation

### 3.1. Selection of the Turbulence Model

According to the assumptions made in Section2.1, the physical quantities did not change in the width direction during the simulation. However, the fracture specimens used in the experiment had certain widths, which would affect the fluid flow. Thus, it was necessary to clarify the influence of the fracture width before choosing the turbulence model. The relationship between the absolute values of the pressure gradient between the real experiment considering the width and the numerical model without considering the width was expressed by the following equation:
(2)Gp=2.01Gp′,whereGp is the absolute value of the pressure gradient in the experiment. Gp′ is the absolute value of the pressure gradient in the numerical model. Equation (3) could be used to get the value of α, that is, 2.01 [31].
(3)a=11−h/D192/π2∑n=1∞tanh2n−1πD/2h/2n−15,whereh is the fracture aperture and D is the fracture width. Equation (2) could be obtained by substituting D of 70 mm and h of 1.8 mm into Equation (3). Equation (2) was used as a two-dimensional model Gp′ of a rough fracture to approximately get the equation of the three-dimensional model.Figure6 shows the change curves of the absolute value of the pressure gradient Gp in the continuous-phase flow in rough fractures in different turbulence models. There was linear loss and partial loss during the seepage test, so the value of Gp was larger than the actual value. Similarly, the value of Gp by the numerical simulation was less than the measured value. As shown in Figure 6, the values of Gp obtained by the RNG k-ε model and Realizable k-ε model were much larger than the measured values, while the values of Gp obtained by using the Shear-Stress Transport (SST) k-ω model, Stress-Omega RSM, and LES were less than the measured values. It could be seen that the measured Gp value changed nonlinearly with the seepage velocity, which was consistent with the change tendencies of simulation results obtained by the SST k-ω model and Stress-Omega RSM. The value of Gp obtained by LES showed an approximately linear relationship with the seepage velocity. Therefore, the SST k-ω model and Stress-Omega RSM could better simulate the pressure loss of the flow in rough fractures. In order to determine the optimal turbulence model, it was necessary to further analyze the structures of their flow fields.Figure 6
Comparison of simulation and test results.LES could directly calculate the turbulent motion larger than the grid scale through instantaneous Navier-Stokes equations, which could predict the formation and distribution of large-scale eddy currents. In this paper, LES, Stress-Omega RSM, and SST k-ω model were used to simulate the water-sediment flow in rough fractures with the fracture inlet velocity of 0.869 m/s. The vortex structure obtained by LES was used as a reference to compare with the results of the other two models. Figures 7 and 8 show the streamline diagrams obtained by the three simulation methods when t=0.06s and t=0.1s, respectively. In Figure 7, when t=0.06s, the streamline diagrams were basically the same obtained by RSM and LES. Specifically, two kinds of vortices were formed at the concave corner of the fracture. One was a three-vortex structure composed of one large and two small vortices, and the large vortex was located downstream of the small vortex. The other was a double-vortex structure composed of one large and one small vortex. The small vortex was located at the bottom of the concave corner of the fracture. However, the form of the vertex obtained by the SST k-ω model was single and was composed of a large vortex and a small vortex at the bottom of the concave corner of the fracture.Figure 7
Streamline diagrams whent=0.06s.
(a)
LES(b)
Stress-Omega RSM(c)
SST k-ω modelFigure 8
Streamline diagrams whent=0.1s.
(a)
LES(b)
Stress-Omega RSM(c)
SST k-ω modelIn Figures7 and 8, the streamline diagrams are different, showing the randomness of the two forms of vortices along the flow direction, while the streamline diagrams obtained by the SST k-ω model were nearly the same when t=0.06s and t=0.1s, which could not reflect the evolution of the vortex over time. In summary, Stress-Omega RSM could accurately calculate the pressure loss and simulate the structure of the flow field in the rough fractures, so it was chosen to simulate the continuous-phase flow in fractures.


### 3.2. Simulation Results of the Water-Sediment Two-Phase Flow in Smooth Fractures

Figure9 shows the change curve of the average pressure at the inlet of the fracture with time when the inlet velocity was 0.869 m/s. The interval from t=0 to t=0.015s was the continuous-phase flow field pressure curve, and the interval from t=0.015s to t=0.5s was the water-sediment two-phase flow field pressure curve after the sand particles were injected.Figure 9
Average pressure-time curve at the inlet of the fracture.In Figure9, the average pressure at the inlet of the fracture increased gradually after the injection of sand particles. With the continuous increase of injected sand particles, it increased to the extreme value and then began to reduce and finally showed small amplitude oscillations within a certain range. It could be seen that the sand particles significantly decreased the average pressure at the inlet. Figure 10 shows the pressure gradient absolute value-velocity curve, and an approximately linear relationship existed.Figure 10
The pressure gradient absolute value-velocity curve.In Figure11, with the inlet velocity of 0.869 m/s, the distributions of velocities, turbulence energy, and pressure on the cross-section of X2=0.9mm (midline position) could be observed when t=0.17s, 0.22 s, and 0.27 s. In Figures 11(a) and 11(b), in the interval of X1≤30mm, the flow velocity and turbulent kinetic energy changed continuously along X2=0.9mm. In the interval of X1≥30mm, they remained stable, indicating that the flow was fully developed. In Figure 11(c), the pressure was approximately linearly distributed along X2=0.9. The velocities, turbulent kinetic energy, and pressure curves were basically consistent at different times on the cross-section of X2=0.9mm, indicating that they did not change over time on this section.Figure 11
Distribution of physical quantities on the cross-section ofX2=0.9mm.
(a)
Velocity(b)
Turbulent kinetic energy(c)
PressureFigure12 shows the distributions of velocities, turbulent kinetic energy, and pressure at different times on the cross-section of X1=40mm. In Figures 13(a) and 13(b), the velocity and turbulent kinetic energy curves basically coincide on the cross-section of X1=40mm, indicating that they were stable and did not change over time. In Figure 14(c), the pressure on the cross-section of X1=40mm changed with time in a small change amplitude.Figure 12
Distribution of physical quantities on the cross-section ofX1=40mm.
(a)
Velocity(b)
Turbulent kinetic energy(c)
PressureFigure 13
Distributions of sand particles: (a)0mm≤X1≤20mm fracture segment; (b) 20mm≤X1≤40mm fracture segment.
(a)(b)Figure 14
Distribution of the momentum source terms with0mm≤X1≤40mm in the discrete-phase model: (a) X1 direction component; (b)X2direction component.
(a)(b)Based on the above analysis, the physical quantities did not change along the flow direction in the interval ofX1≥25mm. Therefore, the interval of 0≤X1≤40mm was chosen to analyze the particle distribution and discrete-phase momentum source term distribution.Figure13 shows the sand distribution in the interval of 0mm≤X1≤40mm, which was divided into two segments. It could be found that there was a layer of particles with lower velocities on the lower wall of the fracture, and there were no particles on the upper wall. It was indicated that particles were deflected downward due to the action of gravity, and a small part of particles was deposited on the lower wall of the fracture.Figure14 shows the distribution of momentum source terms in the discrete-phase model in the fracture interval of 0mm≤X1≤40mm. In Figure 14(a), the absolute value of the momentum source term in the X1 direction component was relatively large at the inlet segment, and the direction was along the negative direction of X1. It was indicated that the force of the continuous-phase fluid on particles was opposite to the flow direction at the inlet segment of the fracture. In the interval of X1≥40mm, the flow was fully developed, and the momentum source terms in the X1 direction component were evenly distributed along the X2 direction. In Figure 14(b), the absolute value of the momentum source term in the X2 direction component was relatively large in the inlet segment. It was positive at the upper part of X2=0.9mm and was negative at the lower part.


### 3.3. Simulation Results of the Water-Sediment Two-Phase Flow in Rough Fractures

#### 3.3.1. Comparison of Numerical Simulation Results and Test Results

Figure15 shows the comparison between the simulated results and the experimental results of the pressure gradient absolute value-seepage velocity curves. Figure 15(a) shows the comparison result of the single phase flow, and Figure 15(b) shows the comparison result of the water-sediment two-phase flow.Figure 15
Comparison of the absolute value of the pressure gradient-seepage velocity curves: (a) single-phase flow; (b) water-sediment two-phase flow.
(a)(b)In Figure15(a), the pressure gradient absolute value-seepage velocity curve of the single phase in fractures obtained by the numerical simulation is basically consistent with that in the experiment, showing a nonlinear relationship. The numerical simulation result was smaller than the experimental result, and the relative error was between 30.1% and 33.4%. This is because Equation (2) describes the relationship between the pressure gradient absolute value of the three-dimensional model of the smooth fracture and the pressure gradient absolute value of the two-dimensional model. The flow was laminar, while in this paper, the flow was turbulent through rough fractures [32]. It is worth noting that the relative errors between the numerical simulation results and the experimental results were relatively close at different flow rates.In Figure15(b), the absolute value of the pressure gradient-seepage velocity of the water-sediment two-phase flow in fractures obtained by the numerical simulation is basically consistent with the curve obtained by the experiment, showing a nonlinear relationship. The numerical simulation result was smaller than the experimental result, and the relative error was between 18.5% and approximately 46.7%. This is because Equation (2) describes the relationship between the absolute value of the pressure gradient in the three-dimensional model of the smooth fracture in laminar flow and that in the two-dimensional model. The flow was a single-phase laminar flow, while this paper focuses on the two-phase turbulence in rough fractures. It should be noted that the absolute error of the numerical simulation results and the test results were relatively close at different flow velocities, and the relative error decreased with the increase of the flow velocity [33].


#### 3.3.2. Change Laws of Physical Quantities in Flow Fields with Time

Figure16 gives the change curve of the average pressure at the inlet with time, when the inlet velocity was 0.869 m/s. It was the pressure curve of the continuous-phase flow field during the time interval 0≤t≤0.12s, and it was the pressure curve of the water-sediment two-phase flow field during the interval 0.12≤t≤0.20s.Figure 16
Average pressure-time curve at fracture inlet.In Figure16, the inlet pressure of the fracture fluctuated violently around 20 kPa, without an obvious decreasing trend. The average value of adjacent peaks and valleys basically reached dynamic stability. This indicated that sand particles had no significant effects on the inlet pressure.Figure17 shows distributions of velocities, turbulent kinetic energy, and pressure on the cross-section of X2=1.4mm at 0.17 s, 0.22 s, and 0.27 s, with the inlet velocity of 0.869 m/s.Figure 17
Distributions of physical quantities on the cross-section ofX2=1.4mm.
(a)
Velocity(b)
Turbulent kinetic energy(c)
PressureIt could be found that whenX1≥2mm and the flow velocity fluctuated around 2.7 m/s, the difference between peak and valley values was small. The velocity curves do not match at different times, and there is no significant difference between the peak value and the valley value. When X2≤2mm, the flow velocity increased sharply. This is because the inlet velocity was uniformly distributed. The fluid velocity recombined from the inlet and changes constantly during the flow through each section. It was decreased due to the viscosity at the walls. The fluid in the middle part outside the boundary layer accelerated.In Figure17(b), when X1≥20mm, the turbulent kinetic energy fluctuated violently between 0.03 m2/s2 and 0.07 m2/s2, which reflected the basic properties of turbulence. When X1≤20mm, it increased sharply, indicating that the fluid velocity recombined from the inlet and increased the turbulence intensity. There was no obvious difference between the peak values and valley values.In Figure17(c), the pressure fluctuates locally and the overall trend decreases linearly with X1.Figure18 shows the distributions of velocity, turbulent kinetic energy, and pressure distribution on the cross-section of X1=40mm, at 0.17 s, 0.22 s, and 0.27 s, with the inlet velocity of 0.869 m/s. In Figure 18(a), the maximum water velocity was between 2.5 m/s and 3.0 m/s, and the corresponding position was between X2=2.2mm and X2=2.4mm. In Figure 18(b), the turbulent kinetic energy was very small near the wall, and the peak values and positions changed dramatically over time. In Figure 18(c), the pressure fluctuated smoothly along X2 and changed greatly with time.Figure 18
Distributions of physical quantities on the cross-section ofX1=40mm: (a) velocity; (b) turbulent kinetic energy; (c) pressure.
(a)(b)(c)


#### 3.3.3. The Spatial Distribution Law of Physical Quantities in the Flow Fields

Figure19 shows the pressure nephogram when t=0.27s, with the inlet velocity of 0.869 m/s. Many lamellate low-pressure and high-pressure areas could be observed in the flow fields. They mainly existed near the fracture tip, indicating that the roughness of the fracture surface increased the pressure loss of the flow field [34].Figure 19
Pressure nephogram: (a)X1=0mm to X1=20mm fracture segment; (b) X1=20mm to X1=40mm fracture segment; (c) X1=40mm to X1=60mm fracture segment; (d) X1=60mm to X1=80mm fracture segment; (e) X1=80mm to X1=100mm fracture segment.
(a)(b)(c)(d)(e)Figure20 shows the flow-line diagram of the fracture segment from X1=30mm to X1=40mm when t=0.27s, with the inlet velocity of 0.869 m/s. The flow separated at the concave corners of the fracture. One or two clockwise-rotating vortices and one counterclockwise-rotating vortex were formed at the concave corners of the lower fracture surface, and one or two counterclockwise-rotating vortices and one clockwise-rotating vortex were formed at the upper fracture surface. Take the lower fracture surface as an example. The causes of the vortices were analyzed. The fluid flowed forward through the fracture. After flowing through the fracture tip, the flow cross-section expanded suddenly. It was impossible for the fluid to suddenly change the direction along the fracture surface due to the inertial force. At this point, smooth transitions of fluid occurred, manifesting that the main flow lines bent near the wall and the flow expanded. Some fluids did not flow forward with the main flow between the outer surface and the wall surface of the expanded part of the main flow. The expansion of the main flow cross-section decreased the flow velocity. At this point, the pressure gradually increased along the flow direction. Thus, the adverse pressure gradient was formed. Some fluids between the main flow and the wall flowed countercurrently along the wall, forming a clockwise vortex. The main flow continued to move forward. When it encountered the next fracture tip, the flow section suddenly shrank, and the main flow section also gradually shrank, forming a clockwise vortex at the wall. It could be found that the downstream vortex range was larger and the flow lines were denser, indicating that the downstream vortex was stronger. Two vortices merged at some concave corners of the fracture. There was a corner vortex rotating counterclockwise at the bottom of concave recessed corner of the fracture, which was formed by two vortices rotating clockwise. There were two types of vortex structures in the flow field. One was the three-vortex structure composed of one large and two small vortices. The other was the double-vortex structure composed of one large and one small vortex. They were randomly distributed along the flow direction.Figure 20
Flow-line diagram.Figure21 shows the velocity nephogram at t=0.27s, with the inlet velocity of 0.869 m/s. When the fluid entered the fracture, the velocity distribution changed drastically. The flow in the fracture was roughly divided into two parts. The main flow was between X2=1mm and X2=1.8mm, and the velocity was between 2.21 m/s and 3.16 m/s. There were many discontinuously distributed high-velocity areas. The main flow was bent at the fracture tip, which was unevenly distributed along the flow direction. The other part of the flow was in the vortex area at the concave corner of the fracture. The fluid velocity was relatively low in the center of the vortex and the boundary layer of the wall [33].Figure 21
The velocity nephogram: (a)X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)Figure22 shows the turbulent kinetic energy nephogram at t=0.27s, with the inlet velocity of 0.869 m/s. It was relatively small at the inlet, because the inlet velocity was evenly distributed, and the flow was not fully developed. As the fluid flowed in the fracture, the turbulent kinetic energy gradually increased. There was a high turbulent kinetic energy region close to the wall at the fracture tip, due to strong collision between the fluid and the fracture tip and high fluctuating velocity. The turbulent kinetic energy was slightly lower near the center line X2=1.4mm. It was the lowest near the wall of the concave corner of the fracture. It was also randomly distributed along the flow.Figure 22
Turbulent kinetic energy nephogram: (a)X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)Figure23 shows the distribution of sand particles in the rough fracture at t=0.27s, with the inlet velocity of 0.869 m/s. In Figures 5–15, the sand particles were densely distributed and the velocity was higher with the 1mm≤X2≤1.8mm segment. The sand particle velocity was higher than the continuous-phase fluid velocity between 2.3 m/s and 3.29 m/s. The sand particle flow direction bent at the fracture tip. The volume concentration of the sand particles was relatively low from the main flow segment to the middle area of the fracture wall. The sand particles were mainly distributed one upstream side of the wall, while the sand concentration was very low on the downstream side of the wall. Within the 0mm≤X1≤20mm segment, the volume fraction of sand particles was relatively high in the lower fracture surface. This was because the flow at the inlet segment was not fully developed, and sand particles were prone to sedimentation due to force of gravity. With the development of flow, the distribution of sand particles tended to be even on the upper and lower fracture surfaces.Figure 23
Distribution of sand particles: (a)X1=0mm to X1=20mm fracture segment; (b) X1=20mm to X1=40mm fracture segment; (c) X1=40mm to X1=60mm fracture segment; (d) X1=60mm to X1=80mm fracture segment; (e) X1=80mm to X1=100mm fracture segment.
(a)(b)(c)(d)(e)Figure24 shows the distribution of the momentum source term of the discrete-phase model in the X1 direction component at t=0.27s, with the inlet velocity of 0.869 m/s. Within the segment of 1mm≤X2≤1.8mm, the momentum source term was relatively large in the X1 direction component, indicating that the force between the phases was large along the flow direction. Within the segment of 0mm≤X1≤20mm, the momentum source term was relatively large in the X1 direction component, indicating that the momentum exchange between phases was large.Figure 24
The momentum source term of the discrete-phase model in theX1 direction component: (a) X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)Figure25 shows the momentum source term of the discrete-phase model in the X2 direction component at t=0.27s, with the inlet velocity of 0.869 m/s. It could be found that it was relatively large when X2=1.4mm. The area with higher absolute value was in a discontinuous and sheet distribution, and positive and negative values appeared alternatively. When the value was positive, the continuous-phase fluid had an upward force on discrete-phase sand particles. When the value was negative, there was a downward force.Figure 25
The momentum source term of the discrete-phase model in theX2 direction component: (a) X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)


## 3.1. Selection of the Turbulence Model

According to the assumptions made in Section2.1, the physical quantities did not change in the width direction during the simulation. However, the fracture specimens used in the experiment had certain widths, which would affect the fluid flow. Thus, it was necessary to clarify the influence of the fracture width before choosing the turbulence model. The relationship between the absolute values of the pressure gradient between the real experiment considering the width and the numerical model without considering the width was expressed by the following equation:
(2)Gp=2.01Gp′,whereGp is the absolute value of the pressure gradient in the experiment. Gp′ is the absolute value of the pressure gradient in the numerical model. Equation (3) could be used to get the value of α, that is, 2.01 [31].
(3)a=11−h/D192/π2∑n=1∞tanh2n−1πD/2h/2n−15,whereh is the fracture aperture and D is the fracture width. Equation (2) could be obtained by substituting D of 70 mm and h of 1.8 mm into Equation (3). Equation (2) was used as a two-dimensional model Gp′ of a rough fracture to approximately get the equation of the three-dimensional model.Figure6 shows the change curves of the absolute value of the pressure gradient Gp in the continuous-phase flow in rough fractures in different turbulence models. There was linear loss and partial loss during the seepage test, so the value of Gp was larger than the actual value. Similarly, the value of Gp by the numerical simulation was less than the measured value. As shown in Figure 6, the values of Gp obtained by the RNG k-ε model and Realizable k-ε model were much larger than the measured values, while the values of Gp obtained by using the Shear-Stress Transport (SST) k-ω model, Stress-Omega RSM, and LES were less than the measured values. It could be seen that the measured Gp value changed nonlinearly with the seepage velocity, which was consistent with the change tendencies of simulation results obtained by the SST k-ω model and Stress-Omega RSM. The value of Gp obtained by LES showed an approximately linear relationship with the seepage velocity. Therefore, the SST k-ω model and Stress-Omega RSM could better simulate the pressure loss of the flow in rough fractures. In order to determine the optimal turbulence model, it was necessary to further analyze the structures of their flow fields.Figure 6
Comparison of simulation and test results.LES could directly calculate the turbulent motion larger than the grid scale through instantaneous Navier-Stokes equations, which could predict the formation and distribution of large-scale eddy currents. In this paper, LES, Stress-Omega RSM, and SST k-ω model were used to simulate the water-sediment flow in rough fractures with the fracture inlet velocity of 0.869 m/s. The vortex structure obtained by LES was used as a reference to compare with the results of the other two models. Figures 7 and 8 show the streamline diagrams obtained by the three simulation methods when t=0.06s and t=0.1s, respectively. In Figure 7, when t=0.06s, the streamline diagrams were basically the same obtained by RSM and LES. Specifically, two kinds of vortices were formed at the concave corner of the fracture. One was a three-vortex structure composed of one large and two small vortices, and the large vortex was located downstream of the small vortex. The other was a double-vortex structure composed of one large and one small vortex. The small vortex was located at the bottom of the concave corner of the fracture. However, the form of the vertex obtained by the SST k-ω model was single and was composed of a large vortex and a small vortex at the bottom of the concave corner of the fracture.Figure 7
Streamline diagrams whent=0.06s.
(a)
LES(b)
Stress-Omega RSM(c)
SST k-ω modelFigure 8
Streamline diagrams whent=0.1s.
(a)
LES(b)
Stress-Omega RSM(c)
SST k-ω modelIn Figures7 and 8, the streamline diagrams are different, showing the randomness of the two forms of vortices along the flow direction, while the streamline diagrams obtained by the SST k-ω model were nearly the same when t=0.06s and t=0.1s, which could not reflect the evolution of the vortex over time. In summary, Stress-Omega RSM could accurately calculate the pressure loss and simulate the structure of the flow field in the rough fractures, so it was chosen to simulate the continuous-phase flow in fractures.


## 3.2. Simulation Results of the Water-Sediment Two-Phase Flow in Smooth Fractures

Figure9 shows the change curve of the average pressure at the inlet of the fracture with time when the inlet velocity was 0.869 m/s. The interval from t=0 to t=0.015s was the continuous-phase flow field pressure curve, and the interval from t=0.015s to t=0.5s was the water-sediment two-phase flow field pressure curve after the sand particles were injected.Figure 9
Average pressure-time curve at the inlet of the fracture.In Figure9, the average pressure at the inlet of the fracture increased gradually after the injection of sand particles. With the continuous increase of injected sand particles, it increased to the extreme value and then began to reduce and finally showed small amplitude oscillations within a certain range. It could be seen that the sand particles significantly decreased the average pressure at the inlet. Figure 10 shows the pressure gradient absolute value-velocity curve, and an approximately linear relationship existed.Figure 10
The pressure gradient absolute value-velocity curve.In Figure11, with the inlet velocity of 0.869 m/s, the distributions of velocities, turbulence energy, and pressure on the cross-section of X2=0.9mm (midline position) could be observed when t=0.17s, 0.22 s, and 0.27 s. In Figures 11(a) and 11(b), in the interval of X1≤30mm, the flow velocity and turbulent kinetic energy changed continuously along X2=0.9mm. In the interval of X1≥30mm, they remained stable, indicating that the flow was fully developed. In Figure 11(c), the pressure was approximately linearly distributed along X2=0.9. The velocities, turbulent kinetic energy, and pressure curves were basically consistent at different times on the cross-section of X2=0.9mm, indicating that they did not change over time on this section.Figure 11
Distribution of physical quantities on the cross-section ofX2=0.9mm.
(a)
Velocity(b)
Turbulent kinetic energy(c)
PressureFigure12 shows the distributions of velocities, turbulent kinetic energy, and pressure at different times on the cross-section of X1=40mm. In Figures 13(a) and 13(b), the velocity and turbulent kinetic energy curves basically coincide on the cross-section of X1=40mm, indicating that they were stable and did not change over time. In Figure 14(c), the pressure on the cross-section of X1=40mm changed with time in a small change amplitude.Figure 12
Distribution of physical quantities on the cross-section ofX1=40mm.
(a)
Velocity(b)
Turbulent kinetic energy(c)
PressureFigure 13
Distributions of sand particles: (a)0mm≤X1≤20mm fracture segment; (b) 20mm≤X1≤40mm fracture segment.
(a)(b)Figure 14
Distribution of the momentum source terms with0mm≤X1≤40mm in the discrete-phase model: (a) X1 direction component; (b)X2direction component.
(a)(b)Based on the above analysis, the physical quantities did not change along the flow direction in the interval ofX1≥25mm. Therefore, the interval of 0≤X1≤40mm was chosen to analyze the particle distribution and discrete-phase momentum source term distribution.Figure13 shows the sand distribution in the interval of 0mm≤X1≤40mm, which was divided into two segments. It could be found that there was a layer of particles with lower velocities on the lower wall of the fracture, and there were no particles on the upper wall. It was indicated that particles were deflected downward due to the action of gravity, and a small part of particles was deposited on the lower wall of the fracture.Figure14 shows the distribution of momentum source terms in the discrete-phase model in the fracture interval of 0mm≤X1≤40mm. In Figure 14(a), the absolute value of the momentum source term in the X1 direction component was relatively large at the inlet segment, and the direction was along the negative direction of X1. It was indicated that the force of the continuous-phase fluid on particles was opposite to the flow direction at the inlet segment of the fracture. In the interval of X1≥40mm, the flow was fully developed, and the momentum source terms in the X1 direction component were evenly distributed along the X2 direction. In Figure 14(b), the absolute value of the momentum source term in the X2 direction component was relatively large in the inlet segment. It was positive at the upper part of X2=0.9mm and was negative at the lower part.


## 3.3. Simulation Results of the Water-Sediment Two-Phase Flow in Rough Fractures

### 3.3.1. Comparison of Numerical Simulation Results and Test Results

Figure15 shows the comparison between the simulated results and the experimental results of the pressure gradient absolute value-seepage velocity curves. Figure 15(a) shows the comparison result of the single phase flow, and Figure 15(b) shows the comparison result of the water-sediment two-phase flow.Figure 15
Comparison of the absolute value of the pressure gradient-seepage velocity curves: (a) single-phase flow; (b) water-sediment two-phase flow.
(a)(b)In Figure15(a), the pressure gradient absolute value-seepage velocity curve of the single phase in fractures obtained by the numerical simulation is basically consistent with that in the experiment, showing a nonlinear relationship. The numerical simulation result was smaller than the experimental result, and the relative error was between 30.1% and 33.4%. This is because Equation (2) describes the relationship between the pressure gradient absolute value of the three-dimensional model of the smooth fracture and the pressure gradient absolute value of the two-dimensional model. The flow was laminar, while in this paper, the flow was turbulent through rough fractures [32]. It is worth noting that the relative errors between the numerical simulation results and the experimental results were relatively close at different flow rates.In Figure15(b), the absolute value of the pressure gradient-seepage velocity of the water-sediment two-phase flow in fractures obtained by the numerical simulation is basically consistent with the curve obtained by the experiment, showing a nonlinear relationship. The numerical simulation result was smaller than the experimental result, and the relative error was between 18.5% and approximately 46.7%. This is because Equation (2) describes the relationship between the absolute value of the pressure gradient in the three-dimensional model of the smooth fracture in laminar flow and that in the two-dimensional model. The flow was a single-phase laminar flow, while this paper focuses on the two-phase turbulence in rough fractures. It should be noted that the absolute error of the numerical simulation results and the test results were relatively close at different flow velocities, and the relative error decreased with the increase of the flow velocity [33].


### 3.3.2. Change Laws of Physical Quantities in Flow Fields with Time

Figure16 gives the change curve of the average pressure at the inlet with time, when the inlet velocity was 0.869 m/s. It was the pressure curve of the continuous-phase flow field during the time interval 0≤t≤0.12s, and it was the pressure curve of the water-sediment two-phase flow field during the interval 0.12≤t≤0.20s.Figure 16
Average pressure-time curve at fracture inlet.In Figure16, the inlet pressure of the fracture fluctuated violently around 20 kPa, without an obvious decreasing trend. The average value of adjacent peaks and valleys basically reached dynamic stability. This indicated that sand particles had no significant effects on the inlet pressure.Figure17 shows distributions of velocities, turbulent kinetic energy, and pressure on the cross-section of X2=1.4mm at 0.17 s, 0.22 s, and 0.27 s, with the inlet velocity of 0.869 m/s.Figure 17
Distributions of physical quantities on the cross-section ofX2=1.4mm.
(a)
Velocity(b)
Turbulent kinetic energy(c)
PressureIt could be found that whenX1≥2mm and the flow velocity fluctuated around 2.7 m/s, the difference between peak and valley values was small. The velocity curves do not match at different times, and there is no significant difference between the peak value and the valley value. When X2≤2mm, the flow velocity increased sharply. This is because the inlet velocity was uniformly distributed. The fluid velocity recombined from the inlet and changes constantly during the flow through each section. It was decreased due to the viscosity at the walls. The fluid in the middle part outside the boundary layer accelerated.In Figure17(b), when X1≥20mm, the turbulent kinetic energy fluctuated violently between 0.03 m2/s2 and 0.07 m2/s2, which reflected the basic properties of turbulence. When X1≤20mm, it increased sharply, indicating that the fluid velocity recombined from the inlet and increased the turbulence intensity. There was no obvious difference between the peak values and valley values.In Figure17(c), the pressure fluctuates locally and the overall trend decreases linearly with X1.Figure18 shows the distributions of velocity, turbulent kinetic energy, and pressure distribution on the cross-section of X1=40mm, at 0.17 s, 0.22 s, and 0.27 s, with the inlet velocity of 0.869 m/s. In Figure 18(a), the maximum water velocity was between 2.5 m/s and 3.0 m/s, and the corresponding position was between X2=2.2mm and X2=2.4mm. In Figure 18(b), the turbulent kinetic energy was very small near the wall, and the peak values and positions changed dramatically over time. In Figure 18(c), the pressure fluctuated smoothly along X2 and changed greatly with time.Figure 18
Distributions of physical quantities on the cross-section ofX1=40mm: (a) velocity; (b) turbulent kinetic energy; (c) pressure.
(a)(b)(c)


### 3.3.3. The Spatial Distribution Law of Physical Quantities in the Flow Fields

Figure19 shows the pressure nephogram when t=0.27s, with the inlet velocity of 0.869 m/s. Many lamellate low-pressure and high-pressure areas could be observed in the flow fields. They mainly existed near the fracture tip, indicating that the roughness of the fracture surface increased the pressure loss of the flow field [34].Figure 19
Pressure nephogram: (a)X1=0mm to X1=20mm fracture segment; (b) X1=20mm to X1=40mm fracture segment; (c) X1=40mm to X1=60mm fracture segment; (d) X1=60mm to X1=80mm fracture segment; (e) X1=80mm to X1=100mm fracture segment.
(a)(b)(c)(d)(e)Figure20 shows the flow-line diagram of the fracture segment from X1=30mm to X1=40mm when t=0.27s, with the inlet velocity of 0.869 m/s. The flow separated at the concave corners of the fracture. One or two clockwise-rotating vortices and one counterclockwise-rotating vortex were formed at the concave corners of the lower fracture surface, and one or two counterclockwise-rotating vortices and one clockwise-rotating vortex were formed at the upper fracture surface. Take the lower fracture surface as an example. The causes of the vortices were analyzed. The fluid flowed forward through the fracture. After flowing through the fracture tip, the flow cross-section expanded suddenly. It was impossible for the fluid to suddenly change the direction along the fracture surface due to the inertial force. At this point, smooth transitions of fluid occurred, manifesting that the main flow lines bent near the wall and the flow expanded. Some fluids did not flow forward with the main flow between the outer surface and the wall surface of the expanded part of the main flow. The expansion of the main flow cross-section decreased the flow velocity. At this point, the pressure gradually increased along the flow direction. Thus, the adverse pressure gradient was formed. Some fluids between the main flow and the wall flowed countercurrently along the wall, forming a clockwise vortex. The main flow continued to move forward. When it encountered the next fracture tip, the flow section suddenly shrank, and the main flow section also gradually shrank, forming a clockwise vortex at the wall. It could be found that the downstream vortex range was larger and the flow lines were denser, indicating that the downstream vortex was stronger. Two vortices merged at some concave corners of the fracture. There was a corner vortex rotating counterclockwise at the bottom of concave recessed corner of the fracture, which was formed by two vortices rotating clockwise. There were two types of vortex structures in the flow field. One was the three-vortex structure composed of one large and two small vortices. The other was the double-vortex structure composed of one large and one small vortex. They were randomly distributed along the flow direction.Figure 20
Flow-line diagram.Figure21 shows the velocity nephogram at t=0.27s, with the inlet velocity of 0.869 m/s. When the fluid entered the fracture, the velocity distribution changed drastically. The flow in the fracture was roughly divided into two parts. The main flow was between X2=1mm and X2=1.8mm, and the velocity was between 2.21 m/s and 3.16 m/s. There were many discontinuously distributed high-velocity areas. The main flow was bent at the fracture tip, which was unevenly distributed along the flow direction. The other part of the flow was in the vortex area at the concave corner of the fracture. The fluid velocity was relatively low in the center of the vortex and the boundary layer of the wall [33].Figure 21
The velocity nephogram: (a)X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)Figure22 shows the turbulent kinetic energy nephogram at t=0.27s, with the inlet velocity of 0.869 m/s. It was relatively small at the inlet, because the inlet velocity was evenly distributed, and the flow was not fully developed. As the fluid flowed in the fracture, the turbulent kinetic energy gradually increased. There was a high turbulent kinetic energy region close to the wall at the fracture tip, due to strong collision between the fluid and the fracture tip and high fluctuating velocity. The turbulent kinetic energy was slightly lower near the center line X2=1.4mm. It was the lowest near the wall of the concave corner of the fracture. It was also randomly distributed along the flow.Figure 22
Turbulent kinetic energy nephogram: (a)X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)Figure23 shows the distribution of sand particles in the rough fracture at t=0.27s, with the inlet velocity of 0.869 m/s. In Figures 5–15, the sand particles were densely distributed and the velocity was higher with the 1mm≤X2≤1.8mm segment. The sand particle velocity was higher than the continuous-phase fluid velocity between 2.3 m/s and 3.29 m/s. The sand particle flow direction bent at the fracture tip. The volume concentration of the sand particles was relatively low from the main flow segment to the middle area of the fracture wall. The sand particles were mainly distributed one upstream side of the wall, while the sand concentration was very low on the downstream side of the wall. Within the 0mm≤X1≤20mm segment, the volume fraction of sand particles was relatively high in the lower fracture surface. This was because the flow at the inlet segment was not fully developed, and sand particles were prone to sedimentation due to force of gravity. With the development of flow, the distribution of sand particles tended to be even on the upper and lower fracture surfaces.Figure 23
Distribution of sand particles: (a)X1=0mm to X1=20mm fracture segment; (b) X1=20mm to X1=40mm fracture segment; (c) X1=40mm to X1=60mm fracture segment; (d) X1=60mm to X1=80mm fracture segment; (e) X1=80mm to X1=100mm fracture segment.
(a)(b)(c)(d)(e)Figure24 shows the distribution of the momentum source term of the discrete-phase model in the X1 direction component at t=0.27s, with the inlet velocity of 0.869 m/s. Within the segment of 1mm≤X2≤1.8mm, the momentum source term was relatively large in the X1 direction component, indicating that the force between the phases was large along the flow direction. Within the segment of 0mm≤X1≤20mm, the momentum source term was relatively large in the X1 direction component, indicating that the momentum exchange between phases was large.Figure 24
The momentum source term of the discrete-phase model in theX1 direction component: (a) X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)Figure25 shows the momentum source term of the discrete-phase model in the X2 direction component at t=0.27s, with the inlet velocity of 0.869 m/s. It could be found that it was relatively large when X2=1.4mm. The area with higher absolute value was in a discontinuous and sheet distribution, and positive and negative values appeared alternatively. When the value was positive, the continuous-phase fluid had an upward force on discrete-phase sand particles. When the value was negative, there was a downward force.Figure 25
The momentum source term of the discrete-phase model in theX2 direction component: (a) X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)


## 3.3.1. Comparison of Numerical Simulation Results and Test Results

Figure15 shows the comparison between the simulated results and the experimental results of the pressure gradient absolute value-seepage velocity curves. Figure 15(a) shows the comparison result of the single phase flow, and Figure 15(b) shows the comparison result of the water-sediment two-phase flow.Figure 15
Comparison of the absolute value of the pressure gradient-seepage velocity curves: (a) single-phase flow; (b) water-sediment two-phase flow.
(a)(b)In Figure15(a), the pressure gradient absolute value-seepage velocity curve of the single phase in fractures obtained by the numerical simulation is basically consistent with that in the experiment, showing a nonlinear relationship. The numerical simulation result was smaller than the experimental result, and the relative error was between 30.1% and 33.4%. This is because Equation (2) describes the relationship between the pressure gradient absolute value of the three-dimensional model of the smooth fracture and the pressure gradient absolute value of the two-dimensional model. The flow was laminar, while in this paper, the flow was turbulent through rough fractures [32]. It is worth noting that the relative errors between the numerical simulation results and the experimental results were relatively close at different flow rates.In Figure15(b), the absolute value of the pressure gradient-seepage velocity of the water-sediment two-phase flow in fractures obtained by the numerical simulation is basically consistent with the curve obtained by the experiment, showing a nonlinear relationship. The numerical simulation result was smaller than the experimental result, and the relative error was between 18.5% and approximately 46.7%. This is because Equation (2) describes the relationship between the absolute value of the pressure gradient in the three-dimensional model of the smooth fracture in laminar flow and that in the two-dimensional model. The flow was a single-phase laminar flow, while this paper focuses on the two-phase turbulence in rough fractures. It should be noted that the absolute error of the numerical simulation results and the test results were relatively close at different flow velocities, and the relative error decreased with the increase of the flow velocity [33].


## 3.3.2. Change Laws of Physical Quantities in Flow Fields with Time

Figure16 gives the change curve of the average pressure at the inlet with time, when the inlet velocity was 0.869 m/s. It was the pressure curve of the continuous-phase flow field during the time interval 0≤t≤0.12s, and it was the pressure curve of the water-sediment two-phase flow field during the interval 0.12≤t≤0.20s.Figure 16
Average pressure-time curve at fracture inlet.In Figure16, the inlet pressure of the fracture fluctuated violently around 20 kPa, without an obvious decreasing trend. The average value of adjacent peaks and valleys basically reached dynamic stability. This indicated that sand particles had no significant effects on the inlet pressure.Figure17 shows distributions of velocities, turbulent kinetic energy, and pressure on the cross-section of X2=1.4mm at 0.17 s, 0.22 s, and 0.27 s, with the inlet velocity of 0.869 m/s.Figure 17
Distributions of physical quantities on the cross-section ofX2=1.4mm.
(a)
Velocity(b)
Turbulent kinetic energy(c)
PressureIt could be found that whenX1≥2mm and the flow velocity fluctuated around 2.7 m/s, the difference between peak and valley values was small. The velocity curves do not match at different times, and there is no significant difference between the peak value and the valley value. When X2≤2mm, the flow velocity increased sharply. This is because the inlet velocity was uniformly distributed. The fluid velocity recombined from the inlet and changes constantly during the flow through each section. It was decreased due to the viscosity at the walls. The fluid in the middle part outside the boundary layer accelerated.In Figure17(b), when X1≥20mm, the turbulent kinetic energy fluctuated violently between 0.03 m2/s2 and 0.07 m2/s2, which reflected the basic properties of turbulence. When X1≤20mm, it increased sharply, indicating that the fluid velocity recombined from the inlet and increased the turbulence intensity. There was no obvious difference between the peak values and valley values.In Figure17(c), the pressure fluctuates locally and the overall trend decreases linearly with X1.Figure18 shows the distributions of velocity, turbulent kinetic energy, and pressure distribution on the cross-section of X1=40mm, at 0.17 s, 0.22 s, and 0.27 s, with the inlet velocity of 0.869 m/s. In Figure 18(a), the maximum water velocity was between 2.5 m/s and 3.0 m/s, and the corresponding position was between X2=2.2mm and X2=2.4mm. In Figure 18(b), the turbulent kinetic energy was very small near the wall, and the peak values and positions changed dramatically over time. In Figure 18(c), the pressure fluctuated smoothly along X2 and changed greatly with time.Figure 18
Distributions of physical quantities on the cross-section ofX1=40mm: (a) velocity; (b) turbulent kinetic energy; (c) pressure.
(a)(b)(c)


## 3.3.3. The Spatial Distribution Law of Physical Quantities in the Flow Fields

Figure19 shows the pressure nephogram when t=0.27s, with the inlet velocity of 0.869 m/s. Many lamellate low-pressure and high-pressure areas could be observed in the flow fields. They mainly existed near the fracture tip, indicating that the roughness of the fracture surface increased the pressure loss of the flow field [34].Figure 19
Pressure nephogram: (a)X1=0mm to X1=20mm fracture segment; (b) X1=20mm to X1=40mm fracture segment; (c) X1=40mm to X1=60mm fracture segment; (d) X1=60mm to X1=80mm fracture segment; (e) X1=80mm to X1=100mm fracture segment.
(a)(b)(c)(d)(e)Figure20 shows the flow-line diagram of the fracture segment from X1=30mm to X1=40mm when t=0.27s, with the inlet velocity of 0.869 m/s. The flow separated at the concave corners of the fracture. One or two clockwise-rotating vortices and one counterclockwise-rotating vortex were formed at the concave corners of the lower fracture surface, and one or two counterclockwise-rotating vortices and one clockwise-rotating vortex were formed at the upper fracture surface. Take the lower fracture surface as an example. The causes of the vortices were analyzed. The fluid flowed forward through the fracture. After flowing through the fracture tip, the flow cross-section expanded suddenly. It was impossible for the fluid to suddenly change the direction along the fracture surface due to the inertial force. At this point, smooth transitions of fluid occurred, manifesting that the main flow lines bent near the wall and the flow expanded. Some fluids did not flow forward with the main flow between the outer surface and the wall surface of the expanded part of the main flow. The expansion of the main flow cross-section decreased the flow velocity. At this point, the pressure gradually increased along the flow direction. Thus, the adverse pressure gradient was formed. Some fluids between the main flow and the wall flowed countercurrently along the wall, forming a clockwise vortex. The main flow continued to move forward. When it encountered the next fracture tip, the flow section suddenly shrank, and the main flow section also gradually shrank, forming a clockwise vortex at the wall. It could be found that the downstream vortex range was larger and the flow lines were denser, indicating that the downstream vortex was stronger. Two vortices merged at some concave corners of the fracture. There was a corner vortex rotating counterclockwise at the bottom of concave recessed corner of the fracture, which was formed by two vortices rotating clockwise. There were two types of vortex structures in the flow field. One was the three-vortex structure composed of one large and two small vortices. The other was the double-vortex structure composed of one large and one small vortex. They were randomly distributed along the flow direction.Figure 20
Flow-line diagram.Figure21 shows the velocity nephogram at t=0.27s, with the inlet velocity of 0.869 m/s. When the fluid entered the fracture, the velocity distribution changed drastically. The flow in the fracture was roughly divided into two parts. The main flow was between X2=1mm and X2=1.8mm, and the velocity was between 2.21 m/s and 3.16 m/s. There were many discontinuously distributed high-velocity areas. The main flow was bent at the fracture tip, which was unevenly distributed along the flow direction. The other part of the flow was in the vortex area at the concave corner of the fracture. The fluid velocity was relatively low in the center of the vortex and the boundary layer of the wall [33].Figure 21
The velocity nephogram: (a)X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)Figure22 shows the turbulent kinetic energy nephogram at t=0.27s, with the inlet velocity of 0.869 m/s. It was relatively small at the inlet, because the inlet velocity was evenly distributed, and the flow was not fully developed. As the fluid flowed in the fracture, the turbulent kinetic energy gradually increased. There was a high turbulent kinetic energy region close to the wall at the fracture tip, due to strong collision between the fluid and the fracture tip and high fluctuating velocity. The turbulent kinetic energy was slightly lower near the center line X2=1.4mm. It was the lowest near the wall of the concave corner of the fracture. It was also randomly distributed along the flow.Figure 22
Turbulent kinetic energy nephogram: (a)X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)Figure23 shows the distribution of sand particles in the rough fracture at t=0.27s, with the inlet velocity of 0.869 m/s. In Figures 5–15, the sand particles were densely distributed and the velocity was higher with the 1mm≤X2≤1.8mm segment. The sand particle velocity was higher than the continuous-phase fluid velocity between 2.3 m/s and 3.29 m/s. The sand particle flow direction bent at the fracture tip. The volume concentration of the sand particles was relatively low from the main flow segment to the middle area of the fracture wall. The sand particles were mainly distributed one upstream side of the wall, while the sand concentration was very low on the downstream side of the wall. Within the 0mm≤X1≤20mm segment, the volume fraction of sand particles was relatively high in the lower fracture surface. This was because the flow at the inlet segment was not fully developed, and sand particles were prone to sedimentation due to force of gravity. With the development of flow, the distribution of sand particles tended to be even on the upper and lower fracture surfaces.Figure 23
Distribution of sand particles: (a)X1=0mm to X1=20mm fracture segment; (b) X1=20mm to X1=40mm fracture segment; (c) X1=40mm to X1=60mm fracture segment; (d) X1=60mm to X1=80mm fracture segment; (e) X1=80mm to X1=100mm fracture segment.
(a)(b)(c)(d)(e)Figure24 shows the distribution of the momentum source term of the discrete-phase model in the X1 direction component at t=0.27s, with the inlet velocity of 0.869 m/s. Within the segment of 1mm≤X2≤1.8mm, the momentum source term was relatively large in the X1 direction component, indicating that the force between the phases was large along the flow direction. Within the segment of 0mm≤X1≤20mm, the momentum source term was relatively large in the X1 direction component, indicating that the momentum exchange between phases was large.Figure 24
The momentum source term of the discrete-phase model in theX1 direction component: (a) X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)Figure25 shows the momentum source term of the discrete-phase model in the X2 direction component at t=0.27s, with the inlet velocity of 0.869 m/s. It could be found that it was relatively large when X2=1.4mm. The area with higher absolute value was in a discontinuous and sheet distribution, and positive and negative values appeared alternatively. When the value was positive, the continuous-phase fluid had an upward force on discrete-phase sand particles. When the value was negative, there was a downward force.Figure 25
The momentum source term of the discrete-phase model in theX2 direction component: (a) X1=0mm to X1=20mm; (b) X1=20mm to X1=40mm; (c) X1=40mm to X1=60mm; (d) X1=60mm to X1=80mm; (e) X1=80mm to X1=100mm.
(a)(b)(c)(d)(e)


## 4. Conclusions

In this paper, the water-sediment two-phase flow theory was used to establish the mechanical model of the water-sediment flow in the fractures. The numerical simulation software was used to simulate the water-sediment two-phase flow in smooth fractures and rough fractures. The conclusions are as follows:(I)
The RNG k-ε model, Realizable k-ε model, SST k-ω model, and Stress-Omega RSM and LES were used to simulate the continuous-phase flow in rough fractures. The simulated results were compared with the experimental data. The comparison results showed that Stress-Omega RSM could accurately calculate the pressure loss and simulate the flow field structures of rough fractures. Thus, it was chosen to simulate the continuous-phase flow in fractures(II)
The water-sediment two-phase flow theory was used to simulate the water-sediment two-phase flow in smooth and rough fractures. The sand densityρp was 2650 kg/m3. The sand particle size Dp was 0.04 mm, and the volume fraction was Φ=4.06%. The physical quantities such as pressure, velocity, and turbulence kinetic energy and the momentum source item, flow diagrams, and particle distributions were given with fracture inlet velocities of 0.349 m/s, 0.532 m/s, 0.697 m/s, and 0.869 m/s, respectively. The simulation results showed that the water-sediment two-phase flow pressure loss was basically the same with the experimental data, verifying the correctness of the numerical simulation method(III)
In the flow fields, the change law of the physical quantities with the time and the law of spatial distribution were analyzed in detail. The results indicated that they did not change with time in the smooth fractured flow fields, while changing continuously with time in the rough fractured flow fields and in a dynamic and stable state. In smooth fractures, there was an asymmetric spatial distribution of physical quantities due to the influence of gravity, while the flow field shows the randomness of spatial distribution in a rough fracture, because water-sediment flow was restrained by the walls



---

*Source: 6235044-2022-05-06.xml*