Abstract
Field evidence indicates that cavities often occur in fractured rocks, especially in a Karst region. Once the immiscible liquid flows into the cavity, the cavity has the immiscible liquid entrapped and results in a low recovery ratio. In this paper, the immiscible liquid transport in cavity-fractures was simulated by Lattice Boltzmann Method (LBM). The interfacial and surface tensions were incorporated by Multicomponent Shan-Chen (MCSC) model. Three various fracture positions were generated to investigate the influence on the irreducible nonwetting phase saturation and displacement time. The influences of fracture aperture and wettability on the immiscible liquid transport were discussed and analyzed. It was found that the cavity resulted in a long displacement time. Increasing the fracture aperture with the corresponding decrease in displacement pressure led to the long displacement time. This consequently decreased the irreducible nonwetting phase saturation. The fracture positions had a significant effect on the displacement time and irreducible saturation. The distribution of the irreducible nonwetting phase was strongly dependent on wettability and fracture position. Furthermore, this study demonstrated that the LBM was very effective in simulating the immiscible two-phase flow in the cavity-fracture.
1. Introduction
The problems of immiscible two-phase flows in fractures are frequently encountered in environmental protection and petroleum industries, such as enhanced oil recovery, geological sequestration of carbon dioxide, and polluted groundwater remediation. Research on the influence of fluid properties on immiscible two-phase flow has been ongoing for decades and has shown that fluid properties (e.g., the capillary number, bond number, viscosity ratio, saturation, and wettability) are very important for immiscible two-phase flows [1–3]. However, in practice, the immiscible two-phase flow in fracture is not only controlled by fluid properties but also influenced by fracture structures.
Past research dealing with immiscible two-phase flow has focused primarily on fluid properties. Joekar-Niasar and Hassanizadeh [4] presented a dynamic pore-network model to investigate the effect of fluids properties on nonequilibrium capillarity effects and found that the immiscible two-phase flow displacement mechanism is highly related to the viscosity ratio and capillary number. The research on the influence of capillary number showed that the displacement time and irreducible saturation are strongly dependent on capillary number [5]. Dou et al. [6] investigated the influence of wettability on interfacial area during the immiscible liquid invasion into a rough fracture and the results showed that the wettability also affects the displacement mechanism and leads to the variation of interfacial area. In fractures, once the displacement pressure (, where is the interfacial tension, is the contact angle, and is the local aperture in fractures) overcomes the threshold pressure, one fluid phase can be displaced by another fluid phase. The displacement pressure varies with the aperture in fractures [7, 8]. The aperture distributions in fractures and the influence on flow have been a topic of considerable interest for many years, because an aperture distribution is essential to the studies of the mechanical and transport properties of rock fractures. Grasmueck et al. [9] used Ground-Penetrating Radar (GPR) in Karst region and found that there are some cavities in vertical fractures. However, only a few researches have investigated the immiscible two-phase flow in cavity-fractures. Zhou et al. [10] used smoothed particle hydrodynamics (SPH) to study the oil displacement in cavity-fractures and found that the ultimate oil recovery is determined by the height of the fracture. A further study reported by Zhou et al. [11] investigated the influences of the height of initial two-phase interface on the oil displacement. The lid driven cavity flow is commonly used as a benchmark problem in the field of traditional computational fluid dynamics (CFD). Patel et al. [12] presented the two-sided cubic lid-driven cavity flow with parallel lids moving in opposite direction at a large Reynolds number (). This study showed that the dissipative nature of the flow inside the cavity increases with increase in the size of corner eddies. Cai et al. [13] used a partial transfixion single fracture to investigate the influence of cavity flow on miscible solute transport and they found that the porous media in cavity have a significant effect on flow regime.
The well-known fractured karst reservoir is composed of many cavity-fractures. Those cavities are much bigger than the fractures that connect with them. For the typical cavity-fracture, the cavity is settled between inlet and outlet fractures (see Figure 1). Varied apertures in fracture always associate with a rough fracture wall. An extreme aperture distribution can be found in cavity-fracture. Due to the displacement pressure depended on the local aperture, the outlet fracture may play a significant role on displacement rate and irreducible saturation. Cai et al. [14] proposed a generalized model of spontaneous imbibition and the results indicated that the variably shaped apertures have an important effect on the spontaneous imbibition.
Modeling immiscible two-phase flow in fractures is a frontier challenge for petroleum resource science. Work by various researchers has presented different numerical methods to couple fluid flow with other phenomenon. Kordilla et al. [15] used the SPH model to simulate surface tension dominated flow on smooth fracture surfaces under various flow conditions. Yang et al. [16] proposed an adaptive circle-fitting approach to calculate the in-plane interfacial curvature during the drainage in fractures. Raeesi and Piri [17] used a 3D mixed-wet random pore-scale network model to investigate the impacts of wettability and trapping on the capillary pressure-saturation-interfacial area relationship during two-phase drainage and imbibition processes. However, it is still difficult for those methods to consider the factors (e.g., interfacial tension, surface tension, wettability, and viscosity) that influence the fluid flow. In recent years, the lattice Boltzmann method (LBM) has attracted much attention because of its simplicity and accuracy in handling complex immiscible two-phase flow. LBM is capable of modeling the fluid/fluid and fluid/solid interactions. There are several available models in LBM for the immiscible two-phase flow simulation, such as the chromodynamics model, the pseudopotential model, free energy model, and He-Shan-Doolen model.
In this paper, the immiscible two-phase flow in cavity-fracture was simulated using lattice Boltzmann method. The interfacial and surface tensions were incorporated by Multicomponent Shan-Chen model. Three various positions of the inlet and outlet fractures (so-called “lower case,” “middle case,” and “upper case”) were generated to investigate the influence on irreducible non-wetting phase saturation and displacement time. The effect of the aperture of the inlet and outlet fractures on the immiscible liquid transport was discussed and analyzed in detail. In addition, by adjusting the interfacial and surface tensions, the wettability of the cavity-fracture wall varied from strongly water-wet to weakly water-wet and then the effect of the wettability was investigated.
2. Methodology
2.1. LBM Fundamentals and Multicomponent Shan-Chen (MCSC) Model
The MCSC model is based on the basic lattice Boltzmann equations. In this study, two distribution functions were used to represent wetting phase and nonwetting phase in cavity-fractures, respectively. The lattice Bhatnagar-Gross-Krook (LBGK) model was implemented for the collision term approximation in each evolution equation. The LBGK model is most widely used model due to its simplicity and accuracy in solving the collision term approximation [18, 19]. The evolution equation of each distribution function satisfies the basic lattice Boltzmann equation of a single relaxation time model,where is the fluid particle distribution function in the th velocity at and denotes either the wetting phase or nonwetting phase. is the dimensionless relaxation time which is related to the kinematic viscosity as . For the D2Q9 model, the nine discrete velocities are given byIn (1), is the local equilibrium distribution which is a specially discretized Maxwell-Boltzmann function with a Taylor expansion and defined aswhere is the lattice sound speed and with being the ratio of lattice spacing and time step . In (3), the weights in the D2Q9 model are given byThe macroscopic fluid density, macroscopic velocity, and the velocity common to the various components are computed fromIt should be noted that is the common averaged velocity of all the fluid components without any additional force. In (6), is a change in velocity field due to the additional force.
In the MCSC model, the effect of total force is incorporated through adding acceleration into the velocity field. In the typical MCSC model, the total force includes the fluid-fluid cohesive force and the fluid-solid adhesive force . In (6), . The initial values of particle distribution function of each component are usually assumed as calculated local equilibrium density distribution for the same component based on initial macroscopic fluid density and velocity.
The fluid-fluid cohesive force is calculated in (8) based on the presence of the other fluid in neighboring lattice sites,where the and denote the wetting phase and the nonwetting phase, respectively, and is a parameter that controls the strength of fluid-fluid interaction.
The fluid-solid adhesive force, , is calculated in (9) based on the presence of the solid phase at neighboring lattice sites,where is an indicator function that is 1 and 0 for a solid and a fluid lattice sites, respectively, and is a parameter that controls the strength between each fluid and a wall. Due to the incorporating the fluid-fluid and fluid-solid interaction force, the actual whole fluid velocity and the whole fluid pressure [20] are defined asThe fluid-fluid and fluid-solid interaction forces in MCSC model can be easily fitted to obtain the interfacial tension and surface tension. The interfacial tension between two immiscible fluids is determined according to Laplace’s Law (, where is the interfacial tension, is the mean interface curvature, and is a pressure difference between the inside and outside of bubbles or droplets). In MCSC model, different interfacial tensions between two immiscible phases can be obtained by adjusting the coefficient [21].
In this paper, the contact angle for a water droplet was changed from 0 to 90 degrees, implying the solid surface ranged from strongly water-wet to weakly water-wet. The solid surface was assumed to be hydrophilic. However, if the contact angle is greater than 90 degree, the solid surface becomes hydrophobic.
2.2. Boundary Condition
Many boundary conditions are available in the literatures [22–24] for the MCSC model, such as the nonslip boundary, the velocity and pressure boundary, and the convection-diffusion boundary. In the MCSC model, any lattice site in the computational domain represents either a solid site or a fluid site. For the solid sites, the half-way bounce-back algorithm collision step is implemented instead of the collision step to provide a nonslip wall boundary condition. In order to simulate the displacement process of the immiscible two-phase flow, the pressure and velocity conditions can be implemented. The pressure reservoir boundary condition is applied at inlet and outlet fractures. The immiscible two-phase flow under this pressure boundary condition can be assumed as a quasiequilibrium process, where the capillary force dominated the flow behavior.
In this work, the algorithms for MCSC model were developed as a C++ code. All the variables are in lattice units such that one lattice unit () is defined as 1 lu, and one time step () is defined as 1 ts. The densities of both the wetting and nonwetting phase were set to 1. The effect of the gravitational field is neglected.
2.3. Model Validation: Immiscible Two-Phase Flows in a 2D Smooth Fracture
For the immiscible two-phase flows in fractures, the wetting-phase covers and moves along the fracture wall, while the nonwetting phase flows in the center of fracture (see Figure 2).
To validate the MCSC model, the velocity distributions of the immiscible two-phase flow in a 2D smooth fracture was investigated. A computational domain of lu2 representing the 2D smooth fracture was used for model validation of the MSCS model implementation. The periodical boundary condition was applied to the outlet and inlet of the fracture. A constant body force was added to each phase. The viscosity ratio has an effect on the velocity distribution of the nonwetting phase and wetting phase [21]. Thus, the kinematic viscosities of nonwetting phase and wetting phase were 1.667 lu2/ts and 0.1667 lu2/ts, respectively. The contact angle between the nonwetting phase and wetting phase was set to .
Assuming the wetting phase flows along the fracture in the region , the wetting phase and nonwetting phase saturations are defined as and , respectively. The analytical solutions for velocity distribution can be derived by solving the Navier-Stokes equations:where , , , and are the kinematic viscosities and densities of the wetting phase and nonwetting phase, respectively. Figure 3 shows that the results of LB simulation and analytical solution. The MCSC model simulation result is very consistent with the analytical solution.
(a) lu2 grid resolution
(b) lu2 grid resolution
(c) lu2 grid resolution
Moreover, to test the grid independence for MCSC model, three different lattice grid resolutions were generated. For the three different lattice grid resolutions, the body force, kinematic viscosity, contact angle, and boundary conditions were constant. The initial region of the wetting phase was set as lu in the different lattice grid resolutions. Figure 3 shows the velocity distributions for the lu2, lu2, and lu2 lattice grid resolutions. The MCSC model simulation results are very consistent with the analytical solution under the different lattice grid resolutions. It means that the MCSC model solution is not gird sensitive. In this study, the sets of grid were chosen depending on the capacity of computer. In Figure 3, it is found that although the body force is constant, the maxium velocity varies with the lattice grid resolution. This is because the fact that the body force is applied in each lattice grid. As the lattice grid resolution increases, the resultant body force increases.
3. Results and Discussion
3.1. Effect of Inlet and Outlet Fractures
To investigate the effect of inlet and outlet fractures on the irreducible nonwetting phase saturation and displacement time, the immiscible liquid transport in cavity-fracture with various fracture positions and apertures was simulated by the MCSC model. In this paper, the cavity was assumed to be a simple lu2 square that previously was occupied by the nonwetting phase. The inlet and outlet fractures had the smooth and paralleled fracture walls and symmetrically distributed on the right and left sides of the cavity, respectively. The cavity aperture was set to 98 lu. For all the simulations, the maximum wetting phase velocity was lu/ts by adjusting the pressure boundary condition and the kinematic viscosities of both phases were 0.1667 lu2/ts. Thus, the dimensionless Re number was less than 23, indicating that there was a laminar flow in the cavity-fracture [25].
Firstly, the processes of the wetting phase displacing the nonwetting phase (imbibition process) are simulated with three various positions of the symmetrical inlet and outlet fractures. The so-called “lower case,” “middle case,” and “upper case” for the three various fracture positions can be found in Figures 4 and 5, respectively. For the three various fracture positions, the apertures are equal to 5 lu. Furthermore, in order to investigate the effect of the fracture aperture on immiscible liquid transport, three apertures for inlet and outlet fractures are set as lu, lu, and lu. The contact angle in this section is equal to 12° by adjusting the interfacial and surface tensions, implying that the solid surface is weakly water-wet. If the nonwetting saturation change over 5000 ts is less than 0.1%, it is assumed that the process of the wetting phase invasion has finished and that the nonwetting phase saturation is equal to the irreducible saturation. This computation takes about five hours to run on the Intel CPU (3.2 GHz) computer for a computational domain of lu2.
(a) ts
(b) ts
(c) ts
(a) Middle case, ts
(b) Upper case, ts
As shown in Figure 4, the cavity is initially saturated by the wetting phase. Due to the pressure boundary condition of the inlet and outlet fracture, the flow direction in cavity-fracture is from left to right. The capillary pressure increases with time. Once the capillary pressure overcomes the entry pressure for the inlet and outlet fracture, the wetting phase invades the cavity. However, the displacement pressure is strongly dependent on the contact angle and aperture. The effect of contact angle and aperture on the displacement will be investigated in the following section.
In Figure 4, the wetting phase invades into the cavity from the inlet fracture and there exists an apparent interface between the nonwetting phase and wetting phase due to the interfacial tension. In this case, as the capillary pressure increases with time, the interface between nonwetting phase and wetting phase moves to the outlet fracture. This process is similar to the imbibition process in porous media. Thus, the process of the wetting phase displacing the nonwetting phase is strongly dependent on the interfacial, surface tensions, and fracture structure. As more wetting phase was injected into the inlet fracture, the nonwetting phase in the cavity decreased. It is found that the nonwetting phase still remains and becomes residual, when the displacement reaches the final steady-state. For this case, the irreducible nonwetting phase saturation is equal to . From Figure 5, it can be seen that the position of the inlet and outlet fracture has an effect on the irreducible nonwetting phase saturation. The irreducible nonwetting phase saturations for the middle and upper case are and , respectively. Thus, the position of the inlet and outlet fracture becomes close to the middle of the cavity, and the irreducible saturations decreases.
The effect of aperture on the irreducible nonwetting phase saturation is investigated. Other two apertures (7 lu and 9 lu) for the fractures are employed, respectively. Figure 6 shows the distributions of the irreducible nonwetting phase for the various positions of the inlet and outlet fractures. It can be seen that the aperture has little effect on the distributions of the irreducible nonwetting phase for the same position of the inlet and outlet fractures. However, the aperture has an important effect on the displacement time. Many studies [10, 12] indicate that the cavity leads to the long displacement time. However, in Figure 6, it is found that both of the aperture and fracture position have an important effect on the displacement time. For each case, as the aperture increases, the displacement time becomes shorter. This is because the fact that the capillary pressure increases with time for the pressure boundary condition on the inlet and outlet fracture. Depending on the definition of the displacement pressure, a larger aperture results in a smaller displacement pressure. This consequently decreases the displacement time. On the other hand, although the domain of the cavity for all the simulations is constant, the ratio of the aperture to cavity decreases due to the decreased aperture. Thus, the process of the wetting phase displacing the nonwetting phase takes less time because of the relative large volume capacity of the cavity. For the upper case, the position of the inlet and outlet fracture is far away from the nonslip boundary of the cavity. As a result, the fracture position moves far away from the middle of the cavity, and the displacement time increases. Furthermore, in Figure 6, due to the smaller displacement pressure, it can be seen that the small nonwetting phase droplet remains at the lower right corner of the cavity. This also increases the irreducible nonwetting phase saturation. However, as the fracture position moves to the bottom of the cavity, the small nonwetting phase droplet disappears. Comparing Figures 5 and 6, it is found that the aperture effects the occurrence of the small nonwetting phase droplet. In Figure 5, the nonwetting phase is continuous for three various fracture positions. As the aperture increases more than 5 lu, the nonwetting phases in the middle and upper cases are discontinuous. The small nonwetting phase droplets are captured in the lower right corner of the cavity. This is mainly due to the increased aperture and varied fracture position, and the displacement pressure decreases.
(a) ,
(b) ,
(c) ,
(d) ,
(e) ,
(f) ,
Figure 7 shows that both of aperture and fracture position have an important effect on the irreducible non-wetting phase saturation. From Figure 7, it can be seen that as the fracture position becomes close to the middle position of the cavity, the irreducible nonwetting phase saturation decreases. This result is consistent with the experimental study reported by Zhou et al. [10].
3.2. Effect of Wettability
The wettability plays an important role on the immiscible two-phase flow and is investigated in the following section. By adjusting the surface tension and keeping the interfacial tension constant, the solid surface becomes strongly water-wet (the contact angle is equal to 85°). The domain of the cavity is the same as that in Section 3.1.
Figure 8 shows the distributions of the irreducible nonwetting phase in various cavity-fractures with the strongly water-wet walls. For the three various case, there are different apertures and fracture positions. Comparing to the cases with the weakly water-wet fracture walls (Figures 4 and 5), it is found that the irreducible nonwetting phase droplets tend to be circles. This is due to the increased surface tension. Thus, the nonwetting phase tends to move away from the strongly water-wet fractures wall and the wetting phase is more easier to wet the fracture wall. For the three various cases in Figure 8, the displacement time decreases as the aperture increases and the fracture position becomes close to the middle of the cavity, which is similar to the result with the weakly water-wet fracture walls. In Figure 8, it is found that when the fracture wall is strongly water wet, the irreducible nonwetting phase saturation decreases as the fracture position becomes close to the middle of the cavity. This is also similar to the result with the weakly water wet fracture walls. For the lower case in Figure 8, the irreducible nonwetting phase saturation increases sharply when the fracture aperture is more than 5 lu. This result is different from the lower case in Figures 5 and 6. The increased contact angle leads to not only the increased displacement pressure but also the weak cohesive force between the nonwetting phase and fracture wall. The weak cohesive force however results in the increment of the irreducible nonwetting phase. In fact, the weak cohesive force is in favor of the movement of the nonwetting phase. Thus, the irreducible nonwetting phase saturations with the strongly water-wet fracture walls are generally lower than that with the weakly water-wet fracture walls. For the lower case, the fracture position is close to the bottom of the cavity, where the flow velocity is relatively high. Furthermore, due to the increased contact angle and the variation of aperture from cavity to fracture, the increased capillary pressure changed the displacement mechanism. The outlet fracture is important for the displacement process. The increased capillary pressure is actually caused by the decreased nonwetting phase pressure and the increased wetting phase pressure. Once the wetting phase occupied the outlet fracture, the nonwetting phase in the cavity will remain as the irreducible nonwetting phase. However, this process is strongly dependent on the wettability and fracture aperture.
(a) , ,
(b) , ,
(c) , ,
(d) , ,
(e) , ,
(f) , ,
(g) , ,
(h) , ,
(i) , ,
4. Conclusions
In this paper, a numerical model based on the lattice Boltzmann method has been developed to study the immiscible two-phase flow in the cavity-fracture. The effects of fracture apertures and positions on the irreducible nonwetting phase saturation and displacement time were discussed in detail. By adjusting the interfacial and surface tensions, the immiscible two-phase flow in the cavity-fracture with the strongly water-wet and the weakly water-wet fracture walls was simulated. This study demonstrated that the LBM was very effective in simulating the immiscible two-phase flow in the cavity-fracture.
Although the computational domain of cavity was constant, the cavity resulted in a long displacement time. For the pressure boundary condition on the inlet and outlet fracture, the capillary pressure in cavity-fractures increased with time. The displacement pressure decreased as the fracture aperture increased. This consequently decreased the irreducible nonwetting phase saturation. The fracture positions had a significant effect on the displacement time and irreducible saturation. The distribution of the irreducible nonwetting phase was strongly dependent on wettability and fracture position.
Nomenclature
| : | Lattice sound speed |
| : | th velocity |
| : | Local aperture in fractures |
| : | Fluid particle distribution function |
| : | Local equilibrium distribution function |
| : | Fluid-fluid cohesive force |
| : | Fluid-solid adhesive force |
| : | Parameter that controls the strength of fluid-fluid interaction |
| : | Parameter that controls the strength between each fluid and a wall |
| : | Body force |
| : | Contact angle |
| : | Pressure difference between the inside and outside of bubbles or droplets |
| : | Displacement pressure |
| : | Whole fluid pressure |
| : | Macroscopic fluid density |
| : | Wetting phase density |
| : | Nonwetting phase density |
| : | Interfacial tension |
| : | Indicator function |
| : | Wetting phase saturation |
| : | Nonwetting phase saturation |
| : | Time step |
| : | Dimensionless relaxation time |
| : | Common averaged velocity of all the fluid components |
| : | Change in velocity field |
| : | Macroscopic velocity |
| : | Actual whole fluid velocity |
| : | Nonwetting phase kinematic viscosity |
| : | Wetting phase kinematic viscosity |
| : | Kinematic viscosity |
| : | Weights |
| : | Lattice spacing. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The study is financially supported by the National Natural Science Foundation of China (Grant nos. 41402217, 41202177, and 51109139) and the Natural Science Foundation of Jiangsu (Grant nos. BK2011110 and BK2012814). The authors thank the two anonymous reviewers for their suggestions in improving the paper.