1. Introduction
The increase in greenhouse gas (GHG) concentration in the atmosphere over the past decades has generated strong interest in developing technologies that help in reducing CO
2 emissions, especially those generated by fossil fuel combustion in heat and power generation [
1,
2,
3]. Carbon dioxide sequestration, a fundamental part of a larger Carbon Capture and Storage (CCS) initiative [
4], is an emerging field of research that is viewed as one of the potential solutions to decrease CO
2 concentration in the atmosphere. CCS is based on two main stages [
5,
6]. Firstly, it requires the capture of carbon dioxide from industrial and energy sources, and, secondly, it needs to transport and store it for long-term sequestration in underground aquifers, which are essentially porous media saturated with saline water isolated from the atmosphere. The penetration of a displacing fluid into a porous medium depends on the capillary and viscous forces at work, as well as the physical properties of the porous medium itself [
7]. The displacement of the non-wetting fluid by the wetting phase is essential and affects CO
2 sequestration in terms of storage efficiency, capacity, and security. However, it is a rather complex phenomenon, difficult to study in laboratory or through purely analytical means. Numerical modeling is a complementary and powerful approach that can greatly help in understanding this phenomenon, provided that the discretized model is capable of capturing the rich physics involved in such a process.
Despite the continuous increase in computational power, modelling multiphase flow in a porous medium is challenging for conventional Computational Fluid Dynamics (CFD) tools due to the complex geometries involved, including small length-scale pores, and, hence, the required use of a computationally prohibitive adapted mesh [
8]. For such meso-scale phenomena, the multiphase Lattice Boltzmann Method (LBM) has proved to be a promising alternative to standard CFD techniques [
8,
9]. Among these techniques, the volume-of-fluid (VOF) [
10], the level set (LS) [
11] and the phase field (PF) [
12] methods are most common. However, according to Yang and Boek [
7], these techniques have numerical instabilities at small-scale interface regions.
The LBM reproduces the macroscopic continuous dynamics of fluids, starting from the statistical approach embedded in the Kinetic Theory, in which the microscopic, discrete nature of a fluid can be described by particle distribution functions while avoiding the necessity to track individual molecules [
9,
10,
11]. Additionally, in comparison to traditional CFD tools, the LBM algorithm is intrinsically parallel and effective in dealing with complex boundary geometries [
13,
14]. Since the very early stages of LBM, several multiphase models have been proposed, namely the color-fluid model, the pseudo-potential model, the free energy model and the mean-field theory model [
15,
16,
17]. Despite the many advantages of LBM, all these multicomponent multiphase models are confined to use for small liquid–gas density ratios due to the presence of spurious currents at the interface, which ultimately promote numerical instabilities [
18,
19,
20,
21]. According to Fakhari and Rahiman [
22], LBM multicomponent multiphase models have a limited range of gas–liquid density ratios (up to 15) due to these numerical instabilities. Generally, if the magnitude of the unphysical currents corresponds to the order of the local flow velocity, the model generates meaningless results [
14].
There have been several attempts to improve the performance of multiphase LBM at high density ratios. For example, Inamuro et al. [
23] proposed using a lattice kinetic scheme, which applies a single equilibrium distribution function calibrating the speed of sound, while Lee et al. [
24] introduced a new method, which is based on a pressure equation and can be applied to a wide range of density ratios. Additionally, in the approach by Lee et al., spurious velocities at the interface are drastically reduced by increasing the surface thickness. Moreover, Yuan and Schaefer [
25] studied the stability of LBM at high density ratios under various equation of states (EOS). It was found that Peng–Robinson (P-R) and Carnahan–Starling (C-S) EOS can model high-density ratios while minimizing the effects of spurious velocities because these EOS allow LBM to correctly capture the fluid flow phenomena in terms of temperature, in addition to other relevant thermodynamic parameters. Furthermore, Fakhari et al. [
26] have found that the flow through the pore domain is extremely sensitive to the choice and accuracy of the boundary conditions, especially the entrance conditions.
Well-known numerical studies on porous media apply multiphase color-fluid LBM [
6,
17,
27,
28] with a density ratio of 1. However, the pseudo-potential model has shown better performance in terms of stability for fluid flows involving higher density contrasts, as, for example, in the case of CO
2 flow in porous media [
25]. According to Li et al. [
29], Kupershtokh et al. [
30], Chen et al. [
31] and Yuan-Schaefer [
25], the performance of the model depends on EOS and the forcing scheme. There are different forcing schemes, such as velocity-shift, exact difference method (EDM) and Guo’s forcing scheme, depending on how interaction force is introduced into LBM. The accuracy and stability of these forcing schemes have been analyzed in several studies [
16,
29,
30,
31]. According to Li et al. [
29], EDM is more stable when relaxation time is less than 1, i.e.,
, while a velocity-shift forcing scheme is more stable when
. Li et al. [
29] also analyzed the effects of different values of
on the coexistence curves. However, in a previous study, Kupershtokh et al. [
30] showed that at
, EDM and velocity-shift forcing schemes are identical. Due to the fact that LBM is still a newly emerging method, there are contradictions between some studies. Therefore, in this study we compare these two popular forcing schemes in terms of numerical stability and accuracy and make recommendations regarding multiphase flows in complex geometries, where a model’s high stability is crucial.
We use P-R EOS in this study, as it is one of the accurate EOS, which helps specify the fluid thermodynamic properties through their acentric factor. We study the displacement process and its efficiency by applying pseudo-potential LBM instead of the widely used color fluid model. Throughout our work, the gas is injected with a uniform inlet profile into a saturated aquifer at different viscosity ratios, capillary numbers and pore geometry. By changing these parameters, their effects on displacement efficiency are investigated. The liquid–gas density ratio is fixed at 5 to ensure negligible spurious velocities at the interface. This apparently low-density ratio is also consistent with carbon dioxide injection condition in deep saline aquifers, which typically occurs at critical or near-critical regimes.
The paper is arranged as follows: The multiphase LBM with the P-R EOS model is proposed in
Section 2.
Section 3 is dedicated to validating the model, and
Section 4 addresses extensive analysis using different fluid properties and pore configuration. The conclusions are presented in
Section 5.
3. Validation of the Model
The validation of our model was carried out in consideration of a capillary pressure test with a known analytical solution. The test considered the injection of a gaseous phase from a plenum into two parallel passages (see
Figure 2). The inlet and outlet boundary conditions were of uniform velocity and constant pressure, respectively, while the top and bottom walls were prescribed as on-grid bounce-back. The passages had different widths (
r1 and
r2) and, correspondingly, different capillary pressures (
Pc1 and
Pc2). Capillary pressure can be calculated by the Young–Laplace equation
where
is the surface tension. Then, the obtained values are compared to the pressure difference (Δ
p) at the inlet (
pin) and the outlet (
pout). Depending on the conditions present, the fluid can either enter or not enter the passages (note that passage 1 is narrower than passage 2). When Δ
p is less than
Pc2, the gas cannot enter any of the passages, while if Δ
p is in the range from
Pc1 to
Pc2, the gaseous phase can only penetrate through passage 2. Finally, when Δ
p is larger than
Pc1, the gas can fully penetrate through both passages.
The simulations were performed with contact angle θ = 70°, viscosity ratio D = 1 and a liquid–gas density ratio of 5, while r1 and r2 were chosen as 15 and 30, respectively. Three distinct pressure differences were simulated, namely Δp1 = 1.04 × 0−3, Δp2 = 2.13 × 10−3 and Δp3 = 4.3 × 10−3, where Pc1 and Pc2 are 2.51 × 10−3 and 1.25 × 10−3, correspondingly. All these quantities are in non-dimensional lattice units (lu).
The computational domain was discretized using 280 lattice units (lu) in the x-direction and 140 lu in the y-direction. The passages 1 and 2 are l00 lu in length by 15 lu-width and 30 lu-width, respectively.
Figure 3 shows the excellent agreement between the numerical results and the expected flow behavior.
4. Results and Discussion
The schematic setup of the porous domain is shown in
Figure 4a,b. Inlet and outlet boundary conditions are uniform velocity and constant pressure, respectively, while the top and bottom walls are prescribed as periodic boundary conditions. The computational domain has a size of 401 × 401 (lu
2) and is filled with staggered columns of circular pores with a radius of 10 lu or square pores at 20 × 20 (lu
2), separated by 10 lu. Each column of pores is spatially arranged uniformly with just a small perturbation (+1 and −1 lu) among contiguous pores to facilitate the break-up of flow symmetry during the simulation and resemble a more realistic flow behavior.
Additionally, one pore is removed from the 1st column to observe its effect on the fingering phenomenon. The uniform size and position of the pores in our study allows us a more precise control of the problem, which is different to what happens in experiments where pores often have different positions, pore sizes or variations in wettability. Initially, the porous medium is filled with the liquid phase, while the gaseous phase is introduced from the left to the right at a constant and uniform velocity. All simulations are conducted at a liquid–gas density ratio of 5, in order to have negligible spurious currents at gas-liquid interfaces [
25].
The penetration of the gaseous phase is expected to depend on multiple dimensionless parameters, namely the capillary number (
Ca), the viscosity ratio (
D) and surface wettability. The capillary number
Ca describes the relationship between viscous and interfacial forces:
where
and
are the average velocity and the dynamic viscosity of the gaseous phase, respectively. The viscosity ratio is defined as the ratio between the viscosity of the liquid and gaseous phases:
where
is the dynamic viscosity of the liquid phase.
Furthermore, a common dimensionless time
defined in Equation (16) is used in the comparison of gas penetration for cases with different injection velocity.
where
is the total time-steps and
H is the width of the domain.
As a first stage in the modeling process, we performed a forcing scheme analysis to ensure its proper selection. For this purpose, the model with surface contact angle
θ = 70°,
Ca = 0.038 and circular pore geometry was selected. In this test, two different forcing schemes (namely velocity-shift and EDM) were analyzed under the same conditions at
τ = 1. As can be seen from
Figure 5, all simulations have approximately the same results, as previously proven in
Section 2.2. Thus, further calculations at
τ = 1 will be presented using the original velocity-shift forcing scheme. Nevertheless, the authors would like to suggest that for low-temperature and low-viscosity ratio fluid flow conditions not within the scope of this paper, the advantages of EDM over velocity-shift should be scrutinized, since the former might possess better numerical stability.
Multiple simulations were performed to identify the effect of Ca on the penetration of the gaseous phase into the porous media. Initially, the viscosity ratio was equal to 1 and the surface contact angle was set to θ = 70°, indicating a wettable condition. Three different capillary numbers were examined in this study, namely Ca1 = 0.038, Ca2 = 0.076 and Ca3 = 0.115. Different Ca can be obtained by varying velocity at the inlet boundary.
Figure 6 and
Figure 7 illustrate the different penetration lengths at the common non-dimensional time
= 0.175 at which the same amount of gas has been injected despite the injection velocity.
Figure 6 and
Figure 7 show that a larger
Ca favors faster penetration of the gaseous phase into the porous media, as it was expected. Additionally, the fingering can be observed for each selected array of pore geometry near the entrance of the domain. Furthermore, at a lower
Ca (
Figure 6a and
Figure 7a), it can be observed that the fingering favored the span-wise direction (with a larger pore throat) due to lower entry pressure, while it moved towards the stream-wise (front) faster with increasing
Ca. Nevertheless, backward fingering, also known as capillary fingering, is hardly appreciated in our simulations, despite previous works observing that phenomenon being driven by the capillary force [
28]. Lenormand et al. [
36] also pointed out that capillary fingering may occur through a large pore throat in backward direction due to low entry pressure, but this backward-fingering is more noticeable for non-homogeneous porous media where the difference in pore throat between the cylinders is higher than that of the spacing used in the present investigation.
It is evident that, depending on pore geometry, the effective gas penetration, measured by length T (
Figure 8), changes slightly, as shown in
Table 3.
Figure 8 illustrates the domain length (L), the slip distance (S) and the effective penetration length (T). In our case the slip distance is the same for all simulations.
Table 3 shows that, for low
Ca, gas penetration is almost the same for both geometries, but as
Ca increases, the penetration through the square-shaped pores is augmented significantly more than for the circular pores. We noted that extremely low
Ca promotes the growth of spurious currents associated with numerical instabilities, consistent with findings by Raeini et al. [
37].
Additionally, our results overcame the limitations found by Lenormand et al. [
36], where it was stated that in a heterogeneous pore domain, a high
Ca will force the fluid to occupy a limited fraction of the domain. Instead, our model has periodic boundary conditions, uniform size and uniform positioning of pores, removing that pitfall.
Furthermore, the effect of the viscosity ratio
D is examined by varying gaseous phase viscosity while keeping the viscosity of the liquid phase constant. Additionally,
Ca and
θ are set as constant and equal to 0.038 and 70°, respectively.
Figure 9 and
Figure 10 show that increasing the viscosity ratio makes the gas injection more difficult and increases the time needed to displace the liquid phase out of the domain. The results corroborate previous studies, indicating that, at low
D, the fingering narrows and spreads sideways faster, while it locally thickens when
D increases [
17]. Moreover, our results show the expected [
27] easier front penetration at lower viscosity ratios, depicting how the intruding phase progresses faster in a finger-like manner as shown in
Figure 9a and
Figure 10a. Our results also corroborated that increasing the viscosity ratio favors a uniform frontal displacement (
Figure 9c and
Figure 10c). It should be noted that pore geometry affects the effective gas penetration length. According to
Table 4, in the circular-shape pore domain, the gas penetration length is shorter than in the square-shape pore domain.
The effect of contact angle on gas penetration is investigated by changing the adhesion parameter between wall and fluid. As was previously discussed, depending on capillary pressure (Equation (13)), the gross evolution of the gaseous phase through the pores might be predicted. Consequently, the surface contact angle effectively affects the displacement of the gaseous phase into the liquid-filled porous domain. Previous studies [
6,
17] on porous media were done on a constant contact angle. Therefore, multiple simulations with various surface contact angles, namely
= 83°,
= 70°,
= 50° and
= 33° at fixed
Ca, were carried out.
Figure 9 demonstrates the evolution of the gaseous phase into the porous domain at
= 0.225. All simulations were performed for wetting conditions, i.e.,
< 90°.
Figure 9 shows how gas penetration length increases while the surface contact angle decreases, as was expected. However, the effective gas penetration in the circular-pore domain (
Figure 11) is higher than that in the square-pore domain (
Figure 12) due to the narrower inter-pore pas-sages (see
Figure 13).
5. Conclusions
A multiphase LBM model was employed to study the influence of capillary number, viscosity ratio and surface contact angle on the hydrodynamics of gas penetration into a homogeneous 2D porous medium saturated with liquid. Previous studies [
6,
17,
27,
28] were done based on the phase-field (PF) method in a framework of a color-fluid multiphase LBM model with a density ratio of 1. In the present work, the Peng–Robinson EOS was successfully integrated into LBM in order to increase the stability of the model. As a result, P-R EOS enables the modeling of the gas penetration phenomenon with a liquid–gas density ratio of 5, resulting in negligible instabilities (i.e., spurious currents) at the interface. Additionally, the validity of the prescribed forcing scheme was assessed as one the main strategies to increase the stability of the model. It was found that at
τ = 1 the velocity-shift and EDM forcing schemes performed similarly in terms of stability, and, therefore, the former was chosen to develop the whole set of simulations in porous media.
The numerical results demonstrated that increasing the capillary number (Ca) and surface wettability enhanced the effective gas penetration at constant viscosity ratio D, while an opposite effect is noticed when increasing D whilst maintaining constant Ca and wettability. Furthermore, the shape of the pores influences fingering length and effective gas penetration. In the present study, the effective gas penetration was higher in the square-pore domain than that in the circle-pore domain by approximately 10% and 4% for explored Ca,max and Dmin, respectively. Moreover, for the largest wettability explored in this study, the effective gas penetration in the square-pore domain was less than 3% than that in the circle-pore domain. Future work is planned to investigate gas–liquid penetration behavior in heterogeneous pore domains while increasing the viscosity ratio. What is more, in future study, we intend to use the EDM forcing scheme for better stability at low relaxation times.