1. Introduction
Reducing carbon emissions has become a global consensus and the demand for clean energy is increasing worldwide [
1,
2]. As an important component of clean fossil energy, coalbed methane has gradually become the focus of global exploration and development [
3], and abandoned mine gas extraction is gradually becoming a new way of coalbed methane extraction. Since the permeability of fractures in abandoned mines is much higher than that of the matrix, fractured gas reservoirs occupy an important position in the coalbed methane extraction process [
4]. The complexity of two-phase fluid flow arises mainly from the complex flow interface between the wetting and non-wetting phases [
5]. Different flow structures imply different flow characteristics. Currently, common flow structure classifications include bubble flow, slug flow, stratified flow, etc. [
6,
7]. Among them, the main feature of slug flow is a series of long bubbles separated by slugs [
8]. It is precisely for this reason that the pressure in the flow direction exhibits significant fluctuations. During the coalbed methane extraction process in abandoned mines, the presence of slug flow is a serious limitation to stable extraction [
9]. Therefore, the investigation of the gas–liquid slug flow in rock fractures and its controlling factors is of great importance to instruct reservoir applications.
In a study of two-phase fluid flow, the flow structure map consisting of the flow structure under the action of different two-phase fluid apparent flow velocities is usually used to reflect the evolution law of the two-phase fluid flow structure [
10,
11]. However, the flow structure map only reflects the flow interface where the flow structure is significantly transformed and fails to reflect the variability of the same flow structure under the action of different two-phase fluid flow rates. For the slug structure, there are also significant differences in the two-phase fluid flow characteristics represented by different slug lengths and the number of slugs in a rock fracture. Therefore, it is necessary to study the evolution of the slug flow structure under the action of different flow ratios to improve the understanding of the two-phase fluid slug flow structure.
Fracture tortuosity variation affects two-phase fluid flow behavior [
12,
13]. The two-phase fluid flow structure in rough fractures is more complex and tortuous [
14,
15].Through visualization tests, it was found that the two-phase fluid has an intermittent flow structure in the rough fracture and the fluid channels in the fracture become more concentrated as the fracture tortuosity and inhomogeneity increase [
16]. Meanwhile, similar investigation results were also obtained using a numerical simulation [
17]. However, the evolution pattern of the slug structure under the influence of fracture tortuosity has not been adequately studied.
The wetting angle is the angle between the liquid–solid interface and the tangent line of the liquid surface at the contact point between the liquid phase and the solid phase [
18]. The variation in it changes the flow structure of the two-phase fluid. The research shows significant differences in flow structures in fractures with different wettability [
19]. For droplet and bubble wettability, related reports have appeared in recent years [
20]. However, little has been reported on the evolution of the slug flow structure in fractures under the action of surface wettability. Meanwhile, due to the non-homogeneity of rock mineral components, the wetting angle of the fracture surface shows obvious non-homogeneity [
21]. However, most of the studies on the wettability of fracture surfaces are based on the average wettability angle of the fracture surface, so the current studies on the two-phase fluid flow structure in rough fractures differ greatly from the two-phase fluid distribution in real rock fractures [
19]. At present, there are still very few research reports on the evolution of two-phase fluid flow structures in rough fractures under the influence of surface wetting angle inhomogeneity.
To study the evolution of slug flow structures in rough rock fractures, the transformation of slug length and slug number with variation in the two-phase fluid flow ratio, fracture roughness and the surface contact angle are investigated by performing computational fluid dynamics modeling. A quantitative characterization method for the evolution of the fracture slug flow structure is proposed. Finally, a simple prediction method for the two-phase fluid flow structure in rough fractures is proposed by comprehensively analyzing the range of values of the characteristic parameters of the slug flow structure.
2. Numerical Model of Two-Phase Flow
2.1. Governing Equations of Two-Phase Fluid Flow
The diffusion effect at the interface is equal for an incompressible and insoluble two-phase Newtonian fluid flow via rock fractures, and the convection velocity complies with the requirement that the divergence be zero [
22,
23,
24]. The Navier–Stokes equations with time dependence and the continuity equation are used as follows:
where
u is the velocity vector,
ρ is the fluid density,
p is the pressure,
I is the identity matrix,
Fg is the gravitational force,
F commprises the other external forces in the system, and
Fst is the surface tension force which can be expressed as [
25]
where
n is the unit normal to the interface,
τ is the surface tension, and
δ is the Dirac delta function.
The following are definitions for fluid density and dynamic viscosity:
where
φ is a smooth continuous level set function,
ρg is the gas density,
ρw is the water density,
μg is the dynamic viscosity of the gas, and
μw is the dynamic viscosity of water.
Due to the current lack of experimental and theoretical studies on gas–water two-phase flow under high-temperature and high-pressure conditions, numerical simulations based on these conditions face challenges in being fully validated through relevant experiments or theoretical models. As a result, many researchers have used the gas–water physical and mechanical parameters under ideal conditions to ensure the feasibility of experiments and improve the reliability of simulation results. For computational fluid dynamics modeling, the ambient temperature is 20 °C, the pressure differential is simulated to be low, and the gas dissolution in water is neglected. Moreover, the atmospheric pressure is much higher than the gas pressure drop, so the compressibility of the gas is not considered.
Table 1 lists the physical properties of the two-phase fluids.
The two-phase fluid flow interface in the level set method is described by the following advection equation:
A constant interface thickness is used for numerical modeling and smooth parameter variation is allowed. A phase reinjection scheme is added to the two-phase fluid flow interface equation to enhance the stability of the solution process. The transfer equation of the fluid interface is expressed as
where ε and
γ are model parameters.
ε is related to the mesh size and controls the thickness of the interface and
γ is related to the velocity field and controls the re-initialization and stabilization amount. In this study,
γ is 1 m/s.
2.2. Flow Channel
As this chapter aims to investigate the effect of the inhomogeneity of the fracture sur-face wettability on the flow behavior of the gas–liquid two-phase flow, a flow channel consisting of elements with different tortuosity is selected for this chapter [
26].
Figure 1 shows a schematic diagram of the flow channel. By dividing the two-phase fluid flow channel into ten regions from the middle and setting the wettability parameters in each of them, the two-phase flow behavior in a non-homogeneous surface wettability rock fracture can be simulated. Meanwhile, to consider the effect of fracture roughness on the slug structure, the distance between the center of the U-shaped channel and the outermost edge (
D) is changed to form different flow channel tortuosity to reflect the change in fracture surface roughness.
2.3. Boundary Condition
In the whole simulation process, the gas and liquid phases are injected from the left side and flow out from the right side of the horizontal fracture. Velocity control is used at the inlet side, where the liquid phase is fixed at 0.5 m/s and the gas phase is set according to the water–gas fluid flow ratio (
n) of 0.5, 1, 1.5, and 2, respectively, to meet the needs of different gas–liquid flow ratio simulations. The two-phase fluid water–gas flow ratio (
n) can be expressed as
The outlet side is uniformly set to a zero-pressure outlet boundary (
P = 0) by using pressure control. To analyze the effect caused by a different surface wettability on the two-phase fluid transport, the wall contact angles were set according to a hydrophilic fracture surface (30°), a no significant wettability fracture surface (90°), a hydrophobic fracture surface (120°), and a superhydrophobic fracture surface (150°). The contact angle formed at each stage of the simulation conforms to Young’s equation. The Young–Laplace equation is expressed as
where
Pc is capillary pressure,
σ is the interfacial tension, and
κ0 is the interface curvature. The slip boundary conditions in numerical modeling are directly related to the wettability of the fracture surface through the two-phase fluid flow. In this study, wall slip is set to the Navier slip condition and the minimum length factor of slip length is set to one.
For the non-homogeneous boundary condition setting of the fracture surface, we assign the surface contact angle parameters to the two fracture wall sections by generating the same mean value but with deviation values (d) of 0°, 5°, 10°, 15°, and 20° to make each section exhibit significant non-homogeneity.
For different fracture tortuosity boundary conditions, the distance (D) between the center of the U-shaped channel and the outermost edge was set to 0.1, 0.2, 0.3, and 0.4 mm to make the flow channels show different tortuosity. In this study, the non-homogeneity of the fracture pore size distribution is disregarded, and all fracture pore sizes are set to 0.1 mm.
2.4. Grid Division and Solver Settings
The generated rough fractures are modeled using Comsol Multiphysics 6.0. In the numerical simulation, the flow parameters of each phase of the two-phase fluid were calculated and the control equations were solved using a separation solver. A generalized minimal residual solver (GMRES) was used to calculate the model parameters with an error of 10−3. To optimize convergence, the GMRES was solved in combination with a preprocessing method. The parallel sparse direct linear solver (PARDISO) was used as a hyperfine mesh preprocessor to investigate the phase initialization and time dependence problems. Since the calculations are time-dependent, the simulation time range is 0–0.02 s with a time step of 0.001 s. To avoid the effect of different rough fractures due to different mesh sizes, a mesh size of 0.3 μm~5 μm is used for the whole simulation. The maximum unit growth rate was 1.13 and the curvature factor was 0.3.
2.5. Model Verification
The validity and scientific validity of using the level set method to study the structural evolution of two-phase fluid flow has been reported in our previous studies. Meanwhile, the structural evolution of two-phase fluid flow in a trickle bed reactor by the level set method was also reported by [
19]. To verify the effectiveness of the level set method in characterizing the slug structure of two-phase fluids, we used the simulation results to compare with different slug structures reported by [
26]. The obtained slug structures are in good agreement with the slug structures found in the experiments (
Figure 2). In this paper, the red color represents water, and the blue color represents gas.
3. Evolutionary Law of Slug Structure
3.1. Effect of Water–Gas Flow Ratio
Figure 3 shows the evolution of the slug structure with the increase in the water–gas flow ratio when the mean value of the contact angle (
θ) of the fracture surface is 90° and the deviation value (
d) is 0°. It can be seen that the two-phase fluid flow structure in the fracture is slug flow under several two-phase fluid flow ratio conditions, but the slug section structure in the fracture also shows variability as the two-phase flow ratio increases. The average length of the gas slug section in the fracture is shortened by nearly two times when the water–gas flow ratio rises from 0.5 to 2. The number of slug sections in the fracture gradually increases and the length of the single liquid slug section gradually decreases as the two-phase fluid flow ratio rises. In addition, when
n = 0.5, the longest gas slug section present in the fracture is about double the length of the shortest gas slug section, while when
n = 2, the lengths between individual gas slug sections in the fracture are basically the same. This indicates that the structural length of gas slug segments in the fracture is more stable when the water–gas flow rate is relatively high. Detailed data for the simulations in this paper are presented in
Appendix A.
For the slug flow structure, the length of a slug section is certainly a quantitative description index and likewise, the number of slug sections is also an important parameter to evaluate the slug flow structure in the fracture. It can be seen that the number of slug segments increases significantly when the water–gas flow ratio increases from 0.5 to 1.5, and the number of slug segments in the fracture begins to stabilize from region to region when the water–gas flow ratio increases from 1.5 to 2. This is because the change in the slug flow structure in the fracture when the water–gas flow ratio increases from 0.5 to 1.5 is mainly reflected in the decrease in the length of the gas slug segments and the increase in the number of segments, while the liquid slug segments are much less affected than the gas slug segments at this stage. At this stage, the slug flow structure shows a significant increase in the number of slug sections. But when the water–gas flow ratio is higher than 1.5 and continues to rise, the fracture is gradually dominated by water when the liquid slug structure is stable, and the main change trend in the slug flow structure in the fracture changes in the direction of increasing the length of the gas slug section and decreasing the number of segments.
3.2. Effect of Mean Contact Angle
Figure 4 shows the two-phase fluid flow structure distribution in the fracture under the effect of the mean surface contact angle when the water–gas flow (
n) ratio is one and the deviation value (
d) is 0°. It can be seen that the two-phase fluid flow structure in the fracture is also slug flow under several average surface contact angle conditions. However, unlike the slug flow structure under the action of the water-to-gas flow ratio, the structural instability of the slug section in the fracture increases significantly with the increase in the average wall contact angle, and the slug section in the fracture starts to become different lengths.
Taking the average contact angle of 90° as the benchmark, it can be found that when the contact angle decreases, the fracture presents as hydrophilic, and the transformation of the slug flow structure in the fracture is mainly reflected in the elevated length of the gas slug section and the transient stratified flow structure of the two-phase fluid near the inlet. Meanwhile, when the contact angle of the fracture is higher than 90°, many short slug segments and even the presence of some slug segments in the fracture begin to convert to a droplet structure.
The evolution of the flow structure for two-phase fluids is analyzed on the basis of the apparent flow velocities of different two-phase fluids. Normally, the two-phase fluid flow structure in the fracture develops gradually along the bubble flow–slug flow–annular flow sequence as the apparent gas flow rate increases. However, the results in
Figure 4 show that the flow structure in the fracture similarly develops gradually along the bubble flow–slug flow–annular flow sequence with the decrease in the average contact angle.
3.3. Effect of Non-Homogeneous Wettability
Table 2 lists the value of the surface contact angle at each stage of fracture.
Figure 5 shows a schematic diagram of the change in the slug structure due to surface wetting heterogeneity.
Figure 6 shows the two-phase fluid flow structure for different surface wettability heterogeneities. There is some variability in the slug flow structure in different regions, and the structural distribution characteristics exhibited are generally consistent with the evolution of the two-phase fluid flow structure under different mean contact angle effects found in the previous section.
The flow structure of the two-phase fluid in the fracture exhibited locally depends mainly on the two-phase fluid flow ratio and the surface contact angle at the local location, and when the fracture surface wettability cannot be reflected by a single contact angle size, the two-phase fluid flow structure in the fracture also shows obvious inhomogeneity. However, unlike the water–gas flow ratio and the average surface contact angle on the number of slug sections in the channel, there is a clear pattern of evolution as the contact angle variance in the fracture surface increases, causing the slug structure to change, but there is no clear trend in this change.
3.4. Effect of Fracture Tortuosity
Figure 7 shows the two-phase fluid flow structure distribution in the fracture under the influence of fracture tortuosity when the water–gas flow ratio is one and the mean value of the wall contact angle is 90° with zero deviation values. It can be seen that under several fracture tortuosity conditions, the two-phase fluid flow structure in the fracture is also slug flow; with an increase in fracture tortuosity, the structural instability of the slug section in the fracture increases significantly, and the slug section in the fracture starts to become different lengths. Meanwhile, when the fracture tortuosity decreases, the transformation of the slug flow structure in the fracture is mainly reflected in the decrease in the number of gas slug segments and the decrease in the gas content in the fracture.
4. Quantitative Characterization Method for Two-Phase Fluid Slug Structure
To investigate the evolution of the slug structure in the fracture, the ratio of the total length of the slug section in the fracture to the total length of the flow channel (
LD) is introduced in this section to describe the slug structure.
Figure 8 shows a schematic diagram of gas plunger segments within a flow channel.
LD in the figure is equal to the ratio of the total length of the two slugs to the total length of the flow channel.
4.1. Effect of Water–Gas Flow Ratio on Slug Structure
To investigate the evolution of the slug structure in the fracture under the effect of the water–gas flow ratio and to better study the flow characteristics of two-phase fluid flow in the fracture, the evolution of the slug structure in the fracture under the effect of different flow ratios is studied.
Figure 9 shows the relationship between the water–gas flow ratio and the total length of the slug section to the total length of the flow channel (
LD) in the fracture.
LD increases linearly as the water–gas flow ratio increases from 0.5 to 2.
Table 3 lists the values of the correlation coefficient (
R2).
These correlation coefficients are all higher than 0.9, indicating a significant linear relationship between the water–gas flow ratio and
LD. The above linear relationship can be expressed as the following equation:
where
a1 and
b1 are fitting coefficients.
4.2. Effect of Mean Contact Angle on Slug Structure
Figure 10 shows the relationship between the mean contact angle and the total length of the slug section to the total length of the flow channel (
LD) in the fracture.
LD decreases linearly as the mean contact angle increases from 30° to 150°.
Table 4 lists the values of the correlation coefficient (
R2). These correlation coefficients are all higher than 0.9, indicating a significant linear relationship between the mean contact angle and
LD. The above linear relationship can be expressed as the following equation:
where
a2 and
b2 are fitting coefficients.
4.3. Effect of Non-Homogeneous Wettability on Slug Structure
Figure 11 shows the evolution of the ratio of the total length of the slug section to the total length of the flow channel in the fracture under the action of different heterogeneities of surface wettability. It can be seen that
LD is of a high level when d is zero and decreases with an increase in
d. However, it is difficult to establish a uniform functional relationship between these two parameters.
Therefore, the effect of the non-homogeneity of surface wetting on the two-phase fluid slug flow structure can be summarized qualitatively. The non-homogeneity of the contact angle of the fracture wall can lead to the difficulty of stable formation of the slug structure and LD decreases; however, LD is difficult to be characterized quantitatively by constructing a function of the non-homogeneous wettability on the fracture surface.
4.4. Effect of Fracture Tortuosity on Slug Structure
Figure 12 shows the relationship between the fracture tortuosity and the total length of the slug section to the total length of the flow channel (
LD) in the fracture.
LD increases linearly as the fracture tortuosity increases from 0.1 to 0.4 mm.
Table 5 lists the values of the correlation coefficient (
R2). These correlation coefficients are all higher than 0.9, indicating a significant linear relationship between the mean contact angle and
LD. The above linear relationship can be expressed as the following equation:
where
a3 and
b3 are fitting coefficients.
4.5. Quantitative Characterization Method for the Evolution of Slug Structure
Equations (10)–(12) show that the relationship between the
LD and the water–gas flow ratio, the mean contact angle, and the flow channel tortuosity satisfies the linear functional form. The
LD of the slug structure can be expressed by Equation (13) [
26].
where
A1,
A2,
A3, and
B1 are the best-fitting coefficients. Based on all simulation results, the best-fit equation is
The correlation coefficient (
R2) is 0.96, showing that Equation (14) can describe the simulation data well. It is worth mentioning that in combining the value of
LD with the flow structure in the corresponding state, it can be found that
LD ranges from 0.4 to 0.6 when there is a slug structure in the fracture. When
LD is lower than 0.4, sparse bubbles or droplet structures appear in the fracture (
Figure 4a), while a large number of gas–liquid layered flow structures exist in the fracture when it is higher than 0.6 and less than one (
Figure 4d).
Based on the above studies, it is known that the water–gas flow ratio, surface wettability, and tortuosity change the slug structure in the fracture by altering the gas-to-water ratio, spatial distribution, and gas-to-water interfacial forces within the fracture.
5. Optimization of Natural Gas Extraction Decongestion Technology
The experiment results show that the formation of slug flow is closely related to the two-phase fluid flow ratio, fracture tortuosity, and fracture surface wettability. To ensure efficient gas production, the gas–water two-phase fluid in the fracture should be maintained in an annular flow structure as far as possible (0.6 <
LD < 1).
Figure 13 shows the relationship between the flow structure and
LD. There are several issues that should be considered for the efficient and stable extraction of natural gas.
The slug structure is more likely to be found in tortuous fractures than in straight fractures. Optimizing the pumping process to reduce the curvature of fractures in complex fracture networks can therefore help to prevent the formation of slug flow. In addition, optimization of the proppant composition is also effective in reducing the surface roughness of the fracture and thus preventing the formation of a slug structure.
The elevated surface wettability drives the formation of the slug structure. Therefore, the effect of surface wettability on the slug structure is reduced by selecting a more suitable proppant to ensure that the fracture opening is increased.