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.
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 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 Figure 1, water and sand flowed into the area formed by two parallel smooth fracture surfaces. The distance between the two fracture surfaces is , and the length of the fracture is . The projection of the upper and lower fracture surfaces on the section is a straight line. The boundary of the computational domain is the inlet section , the upper fracture section , the outlet section , and the lower fracture surface .
2.2.2. The Computational Domain of the Rough Fractures
As shown in Figure 2, 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 is a broken line which divides the surface into 50 equal broken line segments marked by and , . The computational domain is formed by the inlet section , the upper fracture surface , the outlet section , and the lower fracture section .
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:
where is the Reynolds number with the hydraulic diameter of . 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 .
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 the direction. The Stress-Omega RSM was used in this paper. 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 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 that 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 , 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.
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. The 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.
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 energy , 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.
Table 1 shows the material properties. Table 2 shows the interaction parameters of sand particles and the wall. Table 3 shows the parameters among sand particles.
The 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 was 2650 kg/m3, the sand particle size 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 Section 2.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:
where is the absolute value of the pressure gradient in the experiment. 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].
where is the fracture aperture and is the fracture width. Equation (2) could be obtained by substituting of 70 mm and of 1.8 mm into Equation (3). Equation (2) was used as a two-dimensional model of a rough fracture to approximately get the equation of the three-dimensional model.
Figure 6 shows the change curves of the absolute value of the pressure gradient 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 was larger than the actual value. Similarly, the value of by the numerical simulation was less than the measured value. As shown in Figure 6, the values of obtained by the RNG k-ε model and Realizable k-ε model were much larger than the measured values, while the values of 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 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 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.
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 and , respectively. In Figure 7, when , 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.
(a) LES
(b) Stress-Omega RSM
(c) SST k-ω model
(a) LES
(b) Stress-Omega RSM
(c) SST k-ω model
In Figures 7 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 and , 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
Figure 9 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 to was the continuous-phase flow field pressure curve, and the interval from to was the water-sediment two-phase flow field pressure curve after the sand particles were injected.
In Figure 9, 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.
In Figure 11, with the inlet velocity of 0.869 m/s, the distributions of velocities, turbulence energy, and pressure on the cross-section of (midline position) could be observed when , 0.22 s, and 0.27 s. In Figures 11(a) and 11(b), in the interval of , the flow velocity and turbulent kinetic energy changed continuously along . In the interval of , they remained stable, indicating that the flow was fully developed. In Figure 11(c), the pressure was approximately linearly distributed along . The velocities, turbulent kinetic energy, and pressure curves were basically consistent at different times on the cross-section of , indicating that they did not change over time on this section.
(a) Velocity
(b) Turbulent kinetic energy
(c) Pressure
Figure 12 shows the distributions of velocities, turbulent kinetic energy, and pressure at different times on the cross-section of . In Figures 13(a) and 13(b), the velocity and turbulent kinetic energy curves basically coincide on the cross-section of , indicating that they were stable and did not change over time. In Figure 14(c), the pressure on the cross-section of changed with time in a small change amplitude.
(a) Velocity
(b) Turbulent kinetic energy
(c) Pressure
(a)
(b)
(a)
(b)
Based on the above analysis, the physical quantities did not change along the flow direction in the interval of . Therefore, the interval of was chosen to analyze the particle distribution and discrete-phase momentum source term distribution.
Figure 13 shows the sand distribution in the interval of , 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.
Figure 14 shows the distribution of momentum source terms in the discrete-phase model in the fracture interval of . In Figure 14(a), the absolute value of the momentum source term in the direction component was relatively large at the inlet segment, and the direction was along the negative direction of . 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 , the flow was fully developed, and the momentum source terms in the direction component were evenly distributed along the direction. In Figure 14(b), the absolute value of the momentum source term in the direction component was relatively large in the inlet segment. It was positive at the upper part of 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
Figure 15 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.
(a)
(b)
In Figure 15(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 Figure 15(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
Figure 16 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 , and it was the pressure curve of the water-sediment two-phase flow field during the interval .
In Figure 16, 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.
Figure 17 shows distributions of velocities, turbulent kinetic energy, and pressure on the cross-section of at 0.17 s, 0.22 s, and 0.27 s, with the inlet velocity of 0.869 m/s.
(a) Velocity
(b) Turbulent kinetic energy
(c) Pressure
It could be found that when 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 , 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 Figure 17(b), when , the turbulent kinetic energy fluctuated violently between 0.03 m2/s2 and 0.07 m2/s2, which reflected the basic properties of turbulence. When , 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 Figure 17(c), the pressure fluctuates locally and the overall trend decreases linearly with .
Figure 18 shows the distributions of velocity, turbulent kinetic energy, and pressure distribution on the cross-section of , 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 and . 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 and changed greatly with time.
(a)
(b)
(c)
3.3.3. The Spatial Distribution Law of Physical Quantities in the Flow Fields
Figure 19 shows the pressure nephogram when , 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].
(a)
(b)
(c)
(d)
(e)
Figure 20 shows the flow-line diagram of the fracture segment from to when , 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 21 shows the velocity nephogram at , 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 and , 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].
(a)
(b)
(c)
(d)
(e)
Figure 22 shows the turbulent kinetic energy nephogram at , 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 . It was the lowest near the wall of the concave corner of the fracture. It was also randomly distributed along the flow.
(a)
(b)
(c)
(d)
(e)
Figure 23 shows the distribution of sand particles in the rough fracture at , 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 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 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.
(a)
(b)
(c)
(d)
(e)
Figure 24 shows the distribution of the momentum source term of the discrete-phase model in the direction component at , with the inlet velocity of 0.869 m/s. Within the segment of , the momentum source term was relatively large in the direction component, indicating that the force between the phases was large along the flow direction. Within the segment of , the momentum source term was relatively large in the direction component, indicating that the momentum exchange between phases was large.
(a)
(b)
(c)
(d)
(e)
Figure 25 shows the momentum source term of the discrete-phase model in the direction component at , with the inlet velocity of 0.869 m/s. It could be found that it was relatively large when . 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.
(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 was 2650 kg/m3. The sand particle size was 0.04 mm, and the volume fraction was . 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
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.