Abstract
To describe accurately the flow characteristic of fracture scale displacements of immiscible fluids, an incompressible two-phase (crude oil and water) flow model incorporating interfacial forces and nonzero contact angles is developed. The roughness of the two-dimensional synthetic rough-walled fractures is controlled with different fractal dimension parameters. Described by the Navier–Stokes equations, the moving interface between crude oil and water is tracked using level set method. The method accounts for differences in densities and viscosities of crude oil and water and includes the effect of interfacial force. The wettability of the rough fracture wall is taken into account by defining the contact angle and slip length. The curve of the invasion pressure-water volume fraction is generated by modeling two-phase flow during a sudden drainage. The volume fraction of water restricted in the rough-walled fracture is calculated by integrating the water volume and dividing by the total cavity volume of the fracture while the two-phase flow is quasistatic. The effect of invasion pressure of crude oil, roughness of fracture wall, and wettability of the wall on two-phase flow in rough-walled fracture is evaluated.
1. Introduction
Two-phase flow through fractured rock is encountered in many industrial activities such as oil and gas recovery, carbon dioxide storage, and underground oil storage. It is also an important phenomenon related to non-aqueous-phase-liquid (NAPL) contamination of groundwater. A three-dimensional numerical simulator simulating two-phase flow was developed for fractured reservoirs [1] and was later used to simulate water-oil displacement in a fractured network [2].
The research of two-phase flow in fracture networks contained in fractured rock is based on that for single fracture. The numerical models that describe the two-phase flow in a single rough-walled fracture use different size scales and theories. The one-dimensional variable-aperture fracture model was initially derived, and the model was the conceptual basis used to evaluate the effect of aperture variation on two-phase flow [3, 4]. This conceptual model employed a local cubic law (LCL). An active capillary pressure between the wetting and nonwetting phases was developed to simulate two-phase flow in rough-walled rock fractures, from which the relative permeability for the lognormal aperture distribution was calculated [5]. Based on a finite volume implementation of the cubic law and the conservation of mass, a conceptual model of two phases in a single fracture was developed, and the results of simulations were compared with a one-dimensional analytic solution [6]. A composition simulator based on Darcy’s law and Brooks–Corey functional relationship was used to simulate the migration of dense NAPLs through a single fracture in a clay aquitard and to evaluate the effect of fracture aperture [7, 8]. Based on a modified TOUGH2 simulator, a continuum-based approach was used to simulate two-phase flow in a single heterogeneous fracture, and the results were compared with the invasion-percolation model [9, 10]. The multicomponent multiphase lattice Boltzmann method was used to investigate the effect of wettability on the NAPL flow in a microscale single fracture [11, 12].
In recent years, the Navier–Stokes equations and Eulerian techniques have been used to simulate two-phase flow aided in computations from such methods as the level set method, volume of fluid method, and phase field method.
The level set method was adopted in tracking the movement of this interface under this complicated force [13–15]. The level set method was presented to simulate the motion of air bubbles in water and falling water drops in air; this method included large density and viscosity ratios as well as surface tension [16, 17]. The level set method was applied to track the motion of interfaces of pore-scale displacements of two-phase in real porous media [18]. The progressive quasistatic algorithm based on the level set method was applied in an analysis of the matrix-fracture transfer and the shape of the two-phase interface during drainage and imbibition in low permeability rock [19, 20]. The level set method was also employed to simulate quasistatic drainage and imbibition for different contact angles in porous media [21]. The level set method and phase field method were used to simulate two-phase flow with viscosity contrast through complex porous media using COMSOL software [22]. In addition, within the COMSOL platform the level set method was used to model two-phase flow in a 3 mm minichannel [23].
The foregoing methods including the ideal conceptual model, cubic law model, and continuum model have been used in simulating two-phase flow in a single fracture with reasonable success, but there were some shortcomings: () the ideal conceptual model based on the parallel plate theory cannot describe the rough surfaces of real rock fractures; () the classical LCL cannot calculate the fluid inertia and in general overestimates flow through real fractures without considering the wettability of the fracture wall; () the continuum model is unable to capture the characteristic of two-phase flow in an unfilled single fracture; () the lattice Boltzmann method is limited to model two-phase flow in a large-scale single fracture because this method needs very long calculation time; () the effect of the roughness on two-phase flow in an unfilled single fracture is not clearly evaluated.
To prevent crude oil from damaging the underground environment during the future operation period of underground oil storage caverns, the property of two-phase flow around oil storage caverns is a key research object. The two-phase flow in rough-walled fracture is not easily accessible experimentally, so an aperture-scale simulation of the two-phase displacement in a rough-walled fracture can provide valuable information. The main aim of the study was to simulate the two-phase (crude oil and water) flow in a rough-walled fracture using the level set method. A two-dimensional synthetic rough-walled fracture was created, based on the computational fluid dynamics adopted in COMSOL; simulations were performed. The moving interface between crude oil and water was tracked in a large-scale unfilled single fracture by overcoming the shortcomings of the previous methods, and the process of crude oil displacing water was presented under different conditions. The effect of the invasion pressure of crude oil, roughness of the fracture wall, and wettability of the rough wall on the flow characteristics during drainage was evaluated at the same time. The methods and results are described and discussed below.
2. Synthetic Rough-Walled Fractures Based on Fractal Dimension
For the purpose of modeling two-phase flow in rough-walled fractures, synthetic rough fractures were used. These fractures were generated using the software SynFrac, which enables a numerical configuring of fractures of different roughness and average apertures, and it provides much more realistic numerically synthesized fractures than previous methods [25–28]. SynFrac was also used to generate a series of similar fractures with increasing roughness, to allow the effect of roughness on transport in fractures to be evaluated [29].
A series of similar rough-walled fractures with increasing roughness was created by controlling its associated fractal dimension (FD). A 100 mm two-dimensional profile was selected from the data of a three-dimensional fracture, while ensuring that the two-dimensional fracture had no contact point. With FD ranging from 2.0 to 2.4, the mean aperture of the synthetic rough fractures is 1.5 mm. The input parameters for the SynFrac software are listed in Table 1. The resolution (1024 × 1024), standard deviation (1 mm), and anisotropy factor (1.0) are kept constant. The five different fractures are drawn in Figure 1.
3. Numerical Model of Two-Phase Flow with Level Set Method
When the two-phase immiscible fluids are in contact with the wall in an unfilled fracture, the pressure difference at the interface of the two fluids is defined as the capillary pressure written as [3]where is capillary pressure (Pa), the pressure in the NAPL (Pa), the pressure in aqueous phase liquid (Pa), the interfacial force (N/m), the contact angle of the two-phase fluid interface with the fracture wall (rad), and the average radius of curvature in the single rough-walled fracture (m).
To describe the process in which the NAPL displaces the aqueous phase liquid, a specific capillary pressure-saturation curve (Figure 2) at the local position of the fracture was proposed [10, 24]. At this local position for a single fracture, the NAPL entry pressure can be calculated:where is the local aperture of the single fracture (m).
To analyze similar process for a single unfilled rough-walled fracture around an oil storage cavern, it is necessary to simulate the two-phase flow considering the effect of the NAPL invasion pressure and the capillary pressure at the interface of the two-phase fluids. Interface tracking is the key element of the modeling of this two-phase flow.
Level set tracking of the interface uses an auxiliary function on the fixed finite element mesh. The differences in densities, viscosities, and interfacial forces of the two immiscible fluids are considered in level set method [30]. For the horizontal single fracture of this study, the effect of gravity on the two fluids is ignored. The incompressible formation of the Navier–Stokes equations is used to calculate the two-phase flow in an unfilled fracture [30, 31]:where is the average density of the two immiscible fluids (), the velocity vector, the position vector, the time (s), the pressure in fluid (), the average dynamic viscosity of two immiscible fluids (Pa·S), the interfacial force (N/m), and the Laplace operator.
The auxiliary function added to track the interface is defined by [23, 30, 31]where the auxiliary level set function representing the interface is a smooth continuous function; where , the fracture cavity regions are filled with NAPL; where , the regions are part of the interface; where , the regions are filled with aqueous phase liquid. The parameter controls the intensity of reinitialization (m/s) and controls the interface thickness (m).
The density and dynamic viscosity vary smoothly over the interface, and the average density and dynamic viscosity can be defined using the level set function [30, 31]:where and are constants representing the densities of NAPL and aqueous phase liquid; and are constants representing dynamic viscosities of NAPL and aqueous phase liquid.
The interfacial force of the right-hand side in (3) can be calculated using [30]where is the interfacial force coefficient (N/m), the unit normal to the interface, Dirac’s delta-function specifying the interface (1/m), and the second derivative of the level set function.
To avoid poor accuracy in calculating the interfacial force derived from , the formulation of this force is developed further [30, 32]:The initial condition imposed in this model of two-phase flow in a rough-walled fracture iswhere is the position in the fracture and the volume fraction of water in the fracture cavity.
The boundary conditions include those for the inlet and outlet boundaries and the wetted wall boundary of the fracture. The inlet boundary is a constant-pressure boundary, and a pressure equals the invasion pressure of the NAPL. The outlet boundary is also a constant-pressure boundary but with the pressure equal to zero.
The properties of the wetted wall boundary include adopting the slip wall condition and setting the contact angle between the interface of two-phase fluid and the fracture wall. In the level set method, the slip wall condition of the interface along the fracture wall adds a frictional force [30]:where is the unit normal to the wall, the fractional force on the wall (N/m), and the slip length (m) (see Figure 3). In simulations, the slip length is set equal to the size of the mesh element.
The contact angle between the interface and fracture wall (Figure 4) is used to reflect the wettability of the rough-walled fracture [12]. In the level set method, given the contact angle , a weak boundary term is added in the wetted wall boundary condition [30]:The modeling of this two-phase flow using the level set method employs the computational fluid dynamics module of COMSOL Multiphysics and was run on a 2.5 GHz PC with 8 G memory. The solver for the two-phase flow interface was automatically selected for this purpose in COMSOL, and the linear system solver with direct methods was chosen in this paper.
4. Effect of Important Parameters on Two-Phase Flow in Rough-Walled Fracture
Water-sealed underground oil storage caverns are generally artificially excavated in fresh rock below a certain groundwater level [33]. After the excavation of oil storage caverns some fractures in the rock will open or be generated; it is necessary to study drainage process of crude oil invasion of a saturated single rough-walled fracture around oil storage caverns. The effects of invasion pressure, fracture roughness, and wettability of the fracture wall were evaluated by calculating volume-fraction curves of water in the fracture cavity.
4.1. Main Parameters of Simulation
In simulating crude oil displacing water in single unfilled fractures, various roughness was set by allowing the FD to range from 2.0 to 2.4. Three different contact angles were chosen: 30°, 45°, and 60°.
The single fractures were subdivided with triangular meshes using the finite element method. The size of the maximum/minimum element is 0.2 mm/0.01 mm. Hence, for example, if the FD is 2.0, the number of elements is 22716 after triangular subdivision.
The two immiscible fluids, crude oil and water, and the values of their physical properties are listed in Table 2 [34].
4.2. Effect of Invasion Pressure of Crude Oil
Fracture 1 with FD of 2.0 was chosen for the analysis of the effect of the invasion pressure. During the simulation, the contact angle was maintained constant at 45°. The single rough-walled fracture was saturated initially with water. The crude oil invades the fracture from the left-hand side as the invasion pressure increases from 14.45 Pa to 115.6 Pa. The curves of the volume fraction of water in the fracture decrease with increasing invasion pressure (Figure 5). For an average aperture 1.5 mm, (2) gives an entry pressure of 28.9 Pa. While the invasion pressure of the crude oil is less than or equal to 28.9 Pa, the crude oil is unable to break through the fracture. With apertures less than the average aperture, crude oil is trapped in local positions after displacing some water from the fracture. Figure 6 presents the distribution of crude oil trapped in the fracture and the aperture distribution along the fracture for invasion pressures of 14.45 Pa and 28.9 Pa. At the low pressure, crude oil nearly does not displace water, whereas at the higher pressure crude oil is trapped at the position 11 mm where the local aperture is 1.07 mm, and the volume fraction of water is 0.896. As evident in Figure 5, the invasion pressure is greater than 28.9 Pa, the crude oil can break through the fracture. Increasing the invasion pressure from 57.8 Pa to 115.6 Pa leads to faster displacements with times falling from 10.5 s to 3.25 s. The plots of invasion pressure of crude oil and volume fraction of water in the fracture at different times are presented in Figure 7. As the invasion time increases, the shape of the curve approaches that for capillary pressure-saturation (Figure 2). With the level set method, the drainage of water from the saturated fracture has been simulated accurately, and the residual water saturation can be calculated directly.
4.3. Effect of Fracture Roughness
To analyze the effect of fracture roughness on the displacement of water by crude oil, five rough-walled fractures were generated using different FD values (Figure 1) from 2.0 to 2.4 in steps of 0.1; the contact angle of the wall is set to 45°. During drainage, the volume fraction of water decreases with different fracture roughness and FD from 2.1 to 2.4 (Figure 8). With invasion pressure less than or equal to 28.9 Pa crude oil becomes trapped in fracture 2 (FD = 2.1) and fracture 3 (FD = 2.2), and the volume fraction of the water in fracture 2 (FD = 2.1) is less than that in fracture 3 (FD = 2.2). Clearly, the crude oil is more easily trapped in fractures with high roughness. As the invasion pressure rises to 57.8 Pa, the crude oil begins to break through fracture 2 (FD = 2.1) and fracture 3 (FD = 2.2) but is trapped in fracture 4 (FD = 2.3) and fracture 5 (FD = 2.4). With invasion pressures exceeding 115.6 Pa, crude oil can break through all four fractures.
(a)
(b)
(c)
(d)
As crude oil is breaking through the fractures some water still exists in the fracture at that instant (Figure 8). Higher invasion pressures lead to rapid crude oil breaking through the fracture with more water existing in the fracture. As time elapses, the volume fraction of water in the fracture decreases slowly. While the crude oil is breaking through fracture 4 and fracture 5 (at times 1.6 s and 2.85 s), the distributions of the crude oil in the fractures are obtained (Figure 9). In the left half of the fracture 4 (FD = 2.3) almost all water is displaced by crude oil, whereas in the right half of the fracture some water is restricted on the rough wall where local apertures are comparatively large. Because the left half of the fracture 4 is less rough than the right half, the average aperture of the left hand is larger than that of the right hand. The water in the left half was more easily displaced by crude oil which needs low NAPL entry pressure. Unlike fracture 4 some of water is restricted on the rough wall along the fracture 5 (FD = 2.4). The results obtained by modeling the two-phase flow indicates that a large roughness leads to more water restricted between crude oil and the wall of the fracture.
(a)
(b)
To analyze the effect of roughness on the two-phase flow further, the relationship curves of invasion pressure of crude oil and the volume fraction of water in the fracture for different roughness at time 10 s are calculated (Figure 10). The results from modeling show that, with large fracture roughness, water is not easily displaced by crude oil, and increasing the roughness changes the rate of drainage noticeably.
4.4. Effect of Wettability of the Fracture Wall
To evaluate the effect of wettability of the fracture wall on the two-phase flow, three fractures (fractures 1, 3, and 5) with contact angles 30° and 60° were used in the numerical calculations.
The volume fraction of water in fracture 1 (FD = 2.0) shows a decline with contact angles 30° and 60° (Figure 11). A comparison of simulation results presented in Figures 11(a) and 11(b) indicates that as crude oil is trapped in the fracture, the volume fraction of water decreases with increasing contact angle from 30° to 60°. While the crude oil is breaking through the fracture, an increase in contact angle leads to a decrease of the breakthrough time. The simulation results are in accord with (2), from which entry pressure decreases as contact angle increases at each local position. The relation between invasion pressure of crude oil and the volume fraction of water in the fracture for different contact angles at time 15 s is presented in Figure 12. The change in contact angle does not make any significant difference on the shape of curves during drainage.
(a)
(b)
The decreasing curves for the volume fraction of water for fracture 3 (FD = 2.2) with contact angles 30° and 60° are shown in Figure 13. The relationship between the invasion pressure of crude oil and the volume fraction of water in the fracture for different contact angles at time 15 s is given in Figure 14. The increase in contact angle leads to a decrease in the breakthrough time but does not make a large difference on the curves of invasion pressure and volume fraction of water in the fracture. With invasion pressures of 28.9 Pa and 57.8 Pa, the volume fraction of water in the fracture decreases slightly for the contact angles of 30° and 60°.
(a)
(b)
Figure 15 shows a decreasing trend for the volume fraction of water in fracture 5 (FD = 2.4) with contact angle set at 60°. With contact angle of 30°, the numerical calculation of the two-phase flow was not stable. With the level set method, the small contact angle requires a large local entry pressure causing an oscillation in simulation results after crude oil being trapped in the very rough wall of the fracture. The relationship between the invasion pressure of crude oil and volume fraction of water in the fracture for different contact angles at time 10 s (Figure 16) indicates that after crude oil breaks through the fracture, the volume fraction of water decreases noticeably because of the addition 15° of the contact angle.
5. Conclusions
In this study, the level set method was used to simulate the process of crude oil invasion into an unfilled two-dimensional rough-walled fracture saturated by water and to investigate the influences of the fracture roughness and wettability of the fracture wall on the relationship between the invasion pressure of crude oil and the volume fraction of water.
Increasing the invasion pressure decreases the breakthrough time, whereas the water volume fraction in the fracture is inversely proportional to the invasion pressure of crude oil as it breaks through the fracture. Larger roughness leads to more water being restricted between crude oil and the wall of fracture. The breakthrough times of crude oil decrease for large contact angles (high wettability), along with the volume fraction of residual water. Compared with the wettability of the fracture wall, wall roughness significantly affects the relationship between the invasion pressure of crude oil and volume fraction of water.
These simulation results are useful to understanding of the characteristics of two-phase (crude oil and water) flow in fractures around underground oil storage caverns.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The study is financially supported by National Key R&D Program of China (Grant no. 2016YFC0402800), project supported by National Natural Science Foundation of China (Grant no. 51709186), Special Funds for Basic Research and Business Expenses of Central Level Public Welfare Research Institutes (Grant no. Y516030), project supported by Basic Research Program (Natural Science Foundation) of Jiangsu Province (Grant no. BK20170154), and China Postdoctoral Science Foundation funded project (Grant no. 2017M611863). The authors also thank Richard Haase, Ph.D., from Liwen Bianji, Edanz Group China (https://www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.