1. Introduction
In recent years, hydraulic fracturing technology has been increasingly used in unconventional oil and gas production. Implementing the fracturing process to the formation can form artificial fractures to communicate existing natural fractures, laminar facies, and other weak natural facies. As a highly permeable channel for subsurface fluid flow, fracture conductivity is one of the evaluation indicators of fracturing effectiveness, and the level of conductivity is directly related to the oil and gas production capacity. The surface of the fracture is rough, compared with the classical flat plate assumption, the fluid flow in the real uneven fracture space will produce greater energy loss. Moreover, in the process of oil and gas production, the rough fracture closure law differs greatly from that of the flat plate fracture. So it is necessary to carry out researches related to the fluid flow law in the rough fracture in order to accurately predict the change of the flow conductivity of the real fracture.
Many scholars have conducted relevant studies using theoretical derivation, numerical experiments, laboratory experiments, and field tests to understand fracture flow behavior. The current breakdown of fluid flow in a single fracture has mainly focused on the influence of fracture surface roughness [
1,
2,
3], Reynolds number [
4], contact ratio [
5], and slip [
3] on the flow rate.
Gentier et al. (1989) [
6] studied the characteristics of intra-fracture flow at different closure stresses utilizing metal casts. They found that the direction of least resistance to flow is the main flow path. Ishibashi et al. (2015) [
7] characterized intra-fracture flow as a dominant flow channel and established a power exponential relationship between fracture opening, fracture permeability, and fracture length. Li B. et al. (2015) [
8] conducted experiments on fluid flow in a single fracture and found a dominant channel for nonlinear fluid flow in rough fractures, pointing out the effect of size on the flow pattern. Veatch et al. (1965) [
9] studied the seepage pattern of fluid in a turbulent rectangular channel using an experimental method using a rectangular body as a seepage flow channel. Deng (2010, 2011, 2012) [
10,
11,
12] conducted studies that established correlations between laboratory-scale and field large-scale fracture models to predict fracture conductivity under closure stress.
Zimmerman et al. (1992) [
13] proposed a correction factor for the contact area by equating the contact areas of the two fracture plates as circular, elliptical, and irregular shapes. Oron et al. (1998) [
14] used the Navier–Stokes equation to reexamine the fluid percolation pattern in a 2D rough fracture structure. They found that the contact area significantly affects the fluid percolation in the fracture. Gangi et al. (1978) [
15] investigated the relationship between fracture permeability and pressure. They derived a semi-empirical formula for fracture permeability and pressure by equating the rough fracture as a cylindrical projection. Indraratna et al. (2001) [
16] investigated the relationship between permeability and fracture face roughness based on Darcy’s law and verified it experimentally. Zou et al. (2014) [
17] used an experimental method to study the effect of fracture surface dislocation degree and roughness on the inflow capacity of self-supported shear fractures and found that roughness significantly affects the inflow capacity of such fractures only at ground closure stress. The greater the roughness, the better the infiltration capacity at low closure stresses.
Sarifzadeh et al. (2009) [
18] used the finite volume method to study the fluid seepage pattern in a constructed three-dimensional rough fracture structure. Taylor et al. (1999) [
19] used the finite element method to study the fluid seepage pattern in a fractured system under different ratios of equivalent permeability and matrix permeability.
To summarize the previous literature, the fluid flow law within rough fractures needs to focus on the preference of flow paths and the variation of flow velocities due to the unevenness of the walls. In this paper, based on indoor physical simulation experiments, full-diameter downhole cores are prepared as core inclusions into standard actual triaxial fracture specimens, and indoor experiments are carried out to obtain realistic fracture surfaces. The geometric reconstruction of the fracture surface is realized by using 3D scanning and digital core technology to abstract the geometric model of the fracture surface. A 3D fracture model is then established, and the velocity and stress fields of fluid flow are solved by the finite element method. At the same time, the effect of surface roughness parameters on fluid flow in closed fractures is evaluated by parametric analysis of four factors: size, uniformity, density, and shape of the micro-convex body. Flow velocity and flow curvature are used to evaluate the fluid flow pattern within the rough fractures quantitatively. Finally, the flow hindrance of the proppant population at different stages is also analyzed based on the proppant accumulation pattern in rough fractures found in the existing literature. This research provides a basis for a refined study of the closed flow conductivity of fractures, while measures to optimize flow in fractures can be proposed in a targeted manner.
2. Modeling of Liquid Flow through Rough Fractures
Some assumptions used during the model derivation process:
- (1)
The fluid flow in fractures can be described as a single phase, incompressible, steady laminar flow.
- (2)
The fracture has equal width. That is to say, the flow channel of the fluid is only affected by the surface roughness.
- (3)
It is assumed that the fluid has no force on the rock and no creep deformation occurs.
2.1. Flow Equations in Fractures
Different turbulence regions will be formed around the object when the fluid flows through a stationary object. In general, fluid flow in fractures can be described as a single phase, incompressible, steady laminar flow, which is governed by the Navier–Stokes equation as follows [
20,
21,
22]:
where
is the liquid density (kg/m
3),
is the velocity vector (m/s),
is the viscosity (mPa·s), and
is the reduced pressure (MPa). The left-hand item in Equation (1) is the advective acceleration terms, and the two items on the right-hand are the pressure gradient and the viscous force terms. For an incompressible fluid, the mass conservation equation is given below [
23,
24].
Equations (1) and (2) are the governing equations of fluid flow in a planar fracture. Moreover, they are solved using finite element software COMSOL Multiphysics [
25].
2.2. Geometric Treatment
The fracture surface mainly surrounds the fluid flow channel in the fracture, and the fracture morphology directly restricts the flow in the rough fracture. Therefore, clarifying the actual morphology of the fracture surface is the basis for investigating the fluid flow in the fracture.
A shale core column (100 mm × 150 mm) was taken from the well. The shale core column was wrapped I n concrete to form a 300 mm × 300 mm × 300 mm specimen. The core was kept in the middle of the specimen during casting (
Figure 1). After 48 h of resting, the sample was drilled down from the surface to the center of the core and cemented with epoxy resin [
26]. Physical simulations of hydraulic fracturing were conducted using an actual triaxial fracturing system. It should be noted here that this experiment is the research content of our team, and the specific experimental details are described in detail in the references, so we will not repeat them here.
After the experiments were completed, the main fracture surface near the wellbore was taken for 3D scanning to reconstruct and obtain the fracture surface profile (
Figure 2). As shown in
Figure 2, the fracture surface shows noticeable undulating changes. Because the rock is composed of mineral particles and cement of different sizes and strengths, there is apparent roughness on the fracture surface.
The contour and profile maps in
Figure 2 show that the micro-convex bodies on the rough fracture surface are morphologically diverse. The micro-convex bodies of rough fracture surfaces are formed by combining multiple morphologically irregular regions. Discriminated by shape factors, the micro-convex bodies of rough fracture surfaces can be approximated as multiple circular regions of different sizes (
Figure 2d).
2.3. Numerical Treatment
The general hydraulic fracturing process includes fracture initiation, expansion, and closure [
27,
28]. After the fracture is fully opened, the fracture aperture reaches its maximum value, and a certain degree of reduction in fracture aperture occurs under closure stress. Indoor experiments and field data show that the fractures can maintain stable aperture and connectivity, relying on self-supporting capability or proppants to provide a channel for oil and gas seepage [
29,
30].
This paper reconstructs the fracture surfaces formed after the real triaxial hydraulic fracturing test. The rough morphology of the actual fracture surface is entirely irregular, and each fracture surface is unique. There is no universal meaning to carrying out relevant studies on such independent samples that do not have repeatability, so this paper makes appropriate simplifications to the micro-convex body of the rough fracture surface and abstracts a reasonable mathematical model.
For different closure degrees of fracture, fluid flow channels are different, and the flow state is bound to differ. A three-dimensional rough fracture geometry model was developed to demonstrate the effect of the closure degree on the fluid flow in the fracture. In the model, a cube is set as the base fracture. Several cylindrical micro-convex bodies are placed on the lower fracture face to simulate the flow obstruction caused by the surface roughness in the self-proppant fracture. The fluid flows in the fracture under the pressure difference. All boundaries are non-permeable fixed boundaries except for the entrance and exit.
Figure 3 shows a sketch of the computational domain for the simulation of flow around cylinders.
Numerous laboratory tests have demonstrated that in situ cores produce rough fractures. The fracture aperture and roughness are well-known controlling characteristics for a fracture’s fluid conductivity [
31,
32,
33,
34,
35,
36]. Therefore, it is essential to simulate the effect of roughness properties of in situ cores in fluid flow processes. Due to the complex mineral composition and the varying sizes of the mineral grains, the resulting fracture surfaces are also morphologically diverse. Then fracture geometry models are established with the random distribution of micro-convex bodies, which are based on the fracture morphology reconstruction method in
Figure 1. Density (number per unit area), size (average radius), uniformity (standard deviation of radius), and shape of the micro-convex base are analyzed to evaluate the effect of surface roughness parameters on fluid flow in closed fractures.
The finite element software package COMSOL [
25] is adopted to numerically solve the flow (Equation (1)) and the mass conservation (Equation (2)). The numerical model uses a fully coupled algebraic multigrid (AMG) direct solver and Jacobian update every time step. The tetrahedral mesh is used in the model, and mesh refinement is performed around micro-convex bodies to ensure convergence of numerical results.
3. Sensitivity Analysis of Fluid Flow inside the Self-Supporting Fracture
First, the fluid flow in the fracture under different closure degrees is analyzed. The diameter and location of the micro-convex bodies in the model are constant, and only the height of the micro-convex bodies is changed to characterize the different closure degrees of fracture. Based on the substrate morphology of the actual fracture surface after reconstruction, the micro-convex body can be simplified to a cylinder with a radius of 5 mm or less. The radius of the micro-convex body in this model is 3 mm. The heights are 0 mm, 1 mm, 2 mm, 3 mm, 4 mm, and 5 mm, respectively (0 mm indicates a flat plate model without flow obstruction, and 5 mm indicates complete contact at all contact points, namely fracture closure).
Figure 4 shows the fluid flow at the root of a cylindrical micro-convex body in the fracture. Fluid flows under the pressure difference through the micro-convex body in the fracture, and the streamline bends. The streamlines bifurcate and turn after encountering the micro-convex body and gradually converge after bypassing the micro-convex body. It can be seen that the velocity is significantly reduced near the micro-convex body due to the obstruction. Comparing the above figures, the greater the degree of closure, the smaller the void area formed by the micro-convex body above and the upper surface of the fracture, and the more fluid can only flow in the fracture as a bypass flow. Therefore, the higher the micro-convex body, the greater the fracture closure and the more significant the flow rate reduction.
The average velocity of the flow field at the root of the micro-convex body is obtained for each model to compare the average flow velocity for different degrees of fracture closure (
Figure 5). The higher the micro-convex body, the higher the degree of fracture closure and the smaller the average flow velocity. When the height of the micro-convex body is close to or equal to the fracture aperture (Model 5 and Model 6), the fracture is considered close to the fully closed state, and the average flow velocity is minimum.
Figure 6 shows that the fluid is uniform and undisturbed before flowing close to the micro-convex body, and the streamlines are a set of parallel straight lines. When the flow approaches the micro-convex body, the streamline bends due to the obstruction of the micro-convex body. After bypassing the maximum lateral area of the micro-convex body, the streamlines turn together again. The streamlines return to a set of parallel straight lines at a certain distance after the micro-convex body. In
Figure 6a, the micro-convex bodies are distributed in a uniform linear array. The flow channels are mainly the gaps between two adjacent rows of columns. The streamlines can be kept chiefly straight and symmetrically distributed in clusters. In
Figure 6b, the micro-convex bodies are randomly distributed, and the flow channels vary in width. In the red area of the diagram, the fluid is subject to less resistance to flow through the wider channel and has a higher flow rate. Upstream from the micro-convex body along the flow direction, the streamline undergoes a sharp turn in very close proximity. However, there is a streamlined recovery zone downstream of the micro-convex body.
The fluid always flows along with a place of less resistance. That is to say, the liquid flows mainly between the areas with large voids. For an actual fracture surface, the distribution of micro-convex bodies is irregular, the probability that adjacent micro-convex bodies are in a straight line is small, and the formed pore area is usually where nonlinear zigzag flow occurs (as shown in
Figure 6b). The tortuosity increases the length of the flow channel, consuming more energy and thus reducing the velocity.
The fluid percolation tortuosity
in the fracture structure can be considered the result of the combined effect of the percolation streamline tortuosity
. Koponen et al. (1996) [
37] derived the calculation of tortuosity based on the Kozeny–Carman, Navier–Stokes equation, with a tiny cross-sectional body as the basic unit of integration and treating the seepage streamline as a capillary, and derived the equation:
where
is the tortuosity at the position
;
is the instantaneous velocity at the position
, m/s;
is the volume, m
3.
Based on the differential idea, Equation (3) can be simplified as follows:
Duda et al. (2011) [
38] argue Equation (3) can also be expressed as:
In the tortuosity equation, Koponen et al. (1996) [
37] and Duda et al. (2011) [
38] think that the tortuosity is closely related to the instantaneous velocity of the location. However, the previous analysis of fluid seepage in fractures shows that the instantaneous velocity is one of the crucial factors affecting the tortuosity at this location. Still, the ratio of the absolute velocity to the partial velocity in the direction of seepage directly determines the tortuosity at this location.
3.1. Effect of Micro-Convex Body Density on Fluid Flow
The number
of circular micro-convex bodies in a plane of the fixed area is used to differentiate the micro-convex bodies’ density. The flow phenomena were simulated for the number
of 10, 20, 30, 40, and 50, respectively. The velocity and pressure fields for a differential pressure of 1 Pa are shown in
Figure 7 and
Figure 8.
From
Figure 7, it can be seen that as the number of micro-convex bodies increases, the more obstructions the fluid encounters when flowing through the flat plate, and the smaller the average value of the full-field flow velocity. When the number of micro-convex bodies is less than 30, a clear dominant channel can be seen in the flow field (red area in
Figure 7), and the fluid can maintain a high flow velocity. Then, as the number of micro-convex bodies increases, the flow paths continue to bifurcate, and the flow velocity remains low.
For the pressure field (
Figure 8), the leading edge of the pressure gradient shows a clear correlation with the number and location of the micro-convex bodies. The pressure contours in the parallel plate model are approximately parallel to the inlet boundary. Micro-convex bodies accelerate the pressure drop. The higher the number of micro-convex bodies in the flow path, the faster the pressure drop, corresponding to the flow velocity change.
3.2. Effect of Micro-Convex Body Size on Fluid Flow
The flow phenomena were simulated for randomly distributed micro-convex bodies with a mean radius
of 0.5 mm, 1 mm, 2 mm, 3 mm, 4 mm, and 5 mm, respectively. The flow velocity field with a pressure difference of 1 Pa is shown in
Figure 9.
The mean value of the radius of the circular micro-convex body increases gradually from
Figure 9a–f. When the micro-convex body is small (
Figure 9a,b), the flow state is minimally affected, and low flow velocity occurs only near the micro-convex body (mainly concentrated downstream of the flow path). There is no vortex in the whole area, the streamline remains flat, and the fluid flows along the direction parallel to the fracture length toward the distal end. As the radius of the micro-convex body gradually increases, the flow velocity at each point in the flow field is highly uneven. An apparent dominant channel appears at the location where no micro-convex body is gathered, and the streamline gradually bends. Vortices exist in the region near the back stationary point under the micro-convex body, and local turbulence occurs in the flow field. When the radius of the micro-convex body increases to a large size (
Figure 9f), the overall flow velocity at each point is low, the maximum flow velocity in the region is 20 mm/s, and streamline tortuosity is more significant.
The streamline tortuosity and average flow velocity are calculated for different micro-convex sizes, as shown in
Figure 10. In the same flow region, the same number of circular micro-convex bodies are randomly distributed. As the size of micro-convex bodies gradually increases, the streamline tortuosity gradually increases, and the average velocity of the whole field linearly decreases.
3.3. Effect of Micro-Convex Body Uniformity on Fluid Flow
The standard deviation
of the micro-convex radius is used to define the uniformity of the randomly distributed micro-convex bodies. The flow phenomena are simulated for standard deviations of the radius of the micro-convex bodies from 0 m to 0.09 m (with an interval of 0.01 m), respectively. The flow velocity and pressure fields of the six models are shown in
Figure 11 and
Figure 12 for a pressure difference of 1 Pa, respectively.
In
Figure 11, the larger the standard deviation of the micro-convex radius, the more inhomogeneous the size distribution of the micro-convex body in the flow field. As the streamline meanders forward, the greater the obstacle the fluid needs to bypass, the faster the fluid pressure drop and the corresponding decrease in flow rate. In
Figure 12, the significant pressure reduction phenomenon always occurs near the considerable micro-convex bodies, such as the sudden pressure drop in the region where the micro-convex bodies gather in the middle of
Figure 12d and near the largest micro-convex bodies in
Figure 12f.
The streamline tortuosity and average velocity are calculated for different micro-convex body uniformity, as shown in
Figure 13. The streamline meander curvature varies sinusoidally with the homogeneity of the micro-convex body for reasons to be further discussed. With the increase of the standard deviation of the micro-convex body radius, the average flow velocity of the whole field tends to fluctuate and decrease. Significantly, in the two models with a standard deviation of 0.08 m and 0.09 m, the flow velocity decreases, which corresponds to the variation of the pressure field in the above figure.
3.4. Effect of Micro-Convex Body Shape on Fluid Flow
Treating the actual rough fracture surface substrate shape as circular is a simplification. The actual fracture may have an uneven edge of the micro-convex body due to the splitting of particles or the presence of specific flaky mineral components. Therefore, a flow model is developed with a randomly distributed micro-convex body with a square base shape. After superimposing the respective pressure field and streamline diagrams, the flow fields of the circular base model and the square base model under 1 Pa differential pressure are plotted as shown in
Figure 14.
Comparing the two figures in
Figure 14, the streamlines bend after encountering the obstruction. The deflection path of the streamline in
Figure 14a is rounded and slow. In
Figure 14b, the fluid turns sharply after encountering the obstruction, and on the side where the fluid flows in, near the two corners of the square, the pressure drops abruptly, and the streamline close to the micro-convex body shows a bend close to 90°, and the flow becomes less stable.
A randomly distributed geometric model with the number of 10, 20, 30, 40, and 50 micro-convex bodies are generated for the circular and square micro-convex base modes. The variation of streamline tortuosity and flow velocity with the number of micro-convex bodies is calculated in
Figure 15.
From
Figure 15, it can be seen that the streamline tortuosity gradually increases, and the flow velocity gradually decreases as the number of micro-convex bodies increases. The tortuosity of the streamlines of the circular micro-convex model is always more extensive than that of the square micro-convex model. After the curvature of the streamlines, the streamlines are close to the circular form at all positions in the near micro-convex body region for the circular micro-convex body. In this process, there is always a velocity component in the y-direction, and the streamline returns to the flat state more slowly. However, for the square micro-convex body, the streamline has gradually approached a balanced form as it approaches the upper and lower boundaries of the micro-convex body. Therefore, the streamlined tortuosity of the square micro-convex model is smaller than that of the circular micro-convex model. The same reason is given for the lower flow velocity of the circular micro-convex model than that of the square micro-convex model. On rough surfaces, with more curvature of locations that are not zero, the nonlinear flow zones increase.
4. Discussion
In this paper, we study the flow state of the fluid when it encounters micro-convex bodies in a fracture, and it is found that the streamline bends to form a nonlinear zone after facing the micro-convex body.
According to boundary layer theory, the fluid surrounding the flow cylinder can be divided into the boundary layer region and the potential flow region [
39]. The flow within the potential flow zone can be regarded as the flow of an ideal fluid, as shown in
Figure 16. The fluid flows along the
x-axis towards the cylinder and “collides” with the column at stagnation point A, where the fluid velocity is 0. Then the fluid reverses direction and flows against the surface of the cylinder. Thus, a circle of radius r is a streamline. Finally, the fluid flows along the
x-axis again. Point B is also a stagnation point because the fluid has to turn sharply in the horizontal direction, and the fluid velocity can only be equal to zero. Point A is the front stagnation point, and point B is the back stationary point. The streamline immediately adjacent to the cylinder is composed of a circle of radius r with the front and rear
x-axis.
When the fluid flows from point A to point M, the flow velocity increases, and the pressure decreases. The decompression accelerates; the flow velocity from point M to point B falls, and the pressure rises with the ascending deceleration. In other words, the pressure at point B is lower than at point A. The flow from point A to point M is the flow with the pressure gradient, and the flow from point M to point B is the flow against the pressure gradient. Observing the role of the pressure gradient, the force caused by the pressure gradient will push the fluid mass forward with an accelerating effect. Against the impact of the pressure gradient, the fluid close to the cylinder flows to the mainstream of the opposite direction of motion, resulting in the separation of the boundary layer. The boundary layer separation will produce a vortex and, constantly carried away by the mainstream, forms a tail vortex area behind the cylinder. Due to the presence of the vortex, the fluid in the tail vortex area generates a significant friction loss, consumes energy, and generates a large drag loss. Therefore, reducing the inverse pressure gradient can reduce pressure loss and velocity reduction.
Proppants are usually injected into the fracture to maintain the effective opening of a hydraulic fracture. The proppant particles will settle and accumulate during transport. Regarding the accumulation morphology of proppant particles in rough fractures, some scholars have conducted experimental studies, and the experimental results are shown in
Figure 17.
Figure 17 shows that in rough fractures, the layer of proppant packing is separated during the proppant piling-up. The proppant accumulates on the fracture surface in non-uniform agglomerates, leading to particle cavities and proppant pillars forming. Li et al. (2020) [
41] divided the formation of proppant pillars at the fracture surface into three stages, as shown in
Figure 18. The pillars start as randomly sized, then gradually stabilize, and eventually the particles aggregate in large numbers to form more enormous proppant pillars.
Based on the proppant accumulation pattern in the fracture at different stages in
Figure 18, the geometric model of the flow in the fracture is established, and the flow field in the fracture is shown in
Figure 19.
The fluid streamlines bend after encountering the proppant pillar, and there are apparent high-speed channels in the fractures. When the proppant particles start to aggregate, the pillars are of different sizes, and the small pillars have less influence on the flow at this time. As the particles continue to gather, the pillar gradually stabilizes to a uniform size, at which point the interference of pillars to the flow is no longer negligible. When the particles continue to accumulate and form a large proppant pillar, the fluid flow rate decreases significantly after bypassing the proppant pillar. When proppants aggregate to a certain extent, a distinct tail vortex appears downstream of the column. The larger the aggregation area is, the larger the vortex area is. It can be seen that the center of the vortex is a low-flow velocity region, which will cause a more significant flow energy loss. The larger the vortex region, the lower the average flow velocity of the whole field.
Thus, if only flow resistance is considered for the fracturing construction process, it is recommended to use pulsed pumping of proppant to form a uniform clumped proppant pillar as much as possible, which can effectively avoid the non-uniformity of the flow field and slow down the reduction of flow velocity. At the same time, considering the effective support performance of the fracture, balancing the flow resistance and optimal support, is the key to establishing the efficient flow conductivity of the fracture.
5. Conclusions
In this paper, based on the 3D reconstruction results of the rough fracture surface after the fracturing of real cores, the flow model of the rough fracture surface with a random distribution of micro-convex bodies was established. The effects of the density, size, uniformity, and shape of the micro-convex base on the fluid bypass flow were analyzed, and the following conclusions obtained.
Fluid flow velocity varies in fractures with different degrees of closure. When the height of micro-convex bodies is close to or equal to the width of the fracture, the fracture can be considered close to a fully closed state. Almost all the fluid in the fracture has to pass through the obstruction of the micro-convex body in the form of tortuous flow, and the average flow velocity in the fracture is the smallest. The flow tortuosity and the number of vortexes increases with the number (density) and size of micro-convex bodies. The considerable pressure reduction always occurs in the vicinity of the large area of micro-bumps, where the corresponding flow velocity reduction is also more pronounced.
The curvature of the micro-convex body edge affects the extent of the nonlinear region in the flow field, which also indicates that an increase in the roughness of the fracture surface increases the turbulent region within the fracture. The more locations where the curvature is not zero, the more nonlinear zones of flow, and the existence of nonlinear zones will significantly reduce the flow rate and affect the output of oil and gas.
With regard to fracture closures, fracture surface roughness is instead beneficial to maintain fracture width during production. Therefore, the actual oil and gas production process usually increases the number of micro-convex bodies on the fracture face by pumping in proppants. Optimizing the proppants placement pattern to make the tail of the flow-obstructing micro-convex body as close to streamlined as possible can effectively avoid non-uniformity and slow down the flow rate reduction.