# A Pseudopotential Lattice Boltzmann Method for Simulation of Two-Phase Flow Transport in Porous Medium at High-Density and High–Viscosity Ratios

Authors: Eslam Ezzatneshan; Reza Goharimehr
Journal: Geofluids
 (2021)

---

## Abstract

In this work, the capability of a multiphase lattice Boltzmann method (LBM) based on the pseudopotential Shan-Chen (S-C) model is investigated for simulation of two-phase flows through porous media at high-density and high–viscosity ratios. The accuracy and robustness of the S-C LBM are examined by the implementation of the single relaxation time (SRT) and multiple relaxation time (MRT) collision operators with integrating the forcing schemes of the shifted velocity method (SVM) and the exact difference method (EDM). Herein, two equations of state (EoS), namely, the standard Shan-Chen (SC) EoS and Carnahan-Starling (CS) EoS, are implemented to assay the performance of the numerical technique employed for simulation of two-phase flows at high-density ratios. An appropriate modification in the forcing schemes is also used to remove the thermodynamic inconsistency in the simulation of two-phase flow problems studied at low reduced temperatures. The comparative study of these improvements of the S-C LBM is performed by considering an equilibrium state of a droplet suspended in the vapor phase. The solver is validated against the analytical coexistence curve for the liquid-vapor system and the surface tension estimation from the Laplace Law. Then, according to the results obtained, a conclusion has been made to choose an efficient numerical algorithm, including an appropriate collision operator, a realistic EoS, and an accurate forcing scheme, for simulation of multiphase flow transport through a porous medium. The patterns of two-phase flow transport through the porous medium are predicted using the present numerical scheme in different flow conditions defined by the capillary number and the dynamic viscosity ratio. The results obtained for the nonwetting phase saturation, penetration structure of the invading fluid, and the displacement patterns of two-phase flow in the porous medium are comparable with those reported in the literature. The present study demonstrates that the S-C LBM with employing the MRT-EDM scheme, CS EoS, and the modified forcing scheme is efficient and accurate for estimation of the two-phase flow characteristics through the porous medium.


---

## Body

## 1. Introduction

The lattice Boltzmann method (LBM) is known as a powerful mesoscopic numerical scheme for simulation of the multiphase flows transport through the porous medium [1–4]. The kinetic nature of the LBM allows to model the phase separation, interparticle interactions, and interfacial phenomena in multiphase flows [5, 6]. Due to these favorable features, numerical solution methods based on the multiphase LBM do not need an additional algorithm for tracking or capturing the interfacial dynamics, in contrast with the Navier-Stokes-based flow solvers. The efficiency and remarkable versatility of multiphase LBM make this method widely used for simulation of fluid flows in the applications involving interfacial dynamics in porous geometry.Several multiphase models developed in the LBM framework [7–11]. Among the existing models, the pseudopotential Shan-Chen (S-C) model [8, 9] is still attractive due to its capability for simulation of multicomponent multiphase flows and simplicity in implementation with appropriate accuracy and performance [12]. The S-C model has the ability to compute the phase transition [8] and predict properties of multicomponent multiphase systems [9, 13, 14]. Despite the advantages of the S-C model, this scheme suffers from spurious currents in the interfacial region of two-phase flows at high-density ratios [15]. The unphysical dependency of the magnitude of the density ratio to the value of the viscosity ratio is another restriction of the standard S-C model for the simulation of multiphase flows [16]. To apply the S-C LBM for the simulation of multiphase flow transport through porous media, it needs to review the improvements made in this scheme to eliminate the mentioned restrictions.Several techniques have been proposed to improve the stability and accuracy of the LBM incorporated with the pseudopotential S-C model for studying multiphase flows at high-density and –viscosity ratios. Implementation of different realistic equations of state (EoS) [16–18] and using higher-order spatial and temporal discretization [19, 20] are the attempts made for the numerical solution of multiphase flow problems at high-density ratios. Another category includes the extension of multiple relaxation time (MRT) technique for S-C LBM [21, 22] to eliminate the numerical instability of this approach for simulation of flow problems with high-viscosity ratios. The unphysical dependency of the viscosity-ratio to the value of the density-ratio in the S-C LBM is also resolved by employing appropriate forcing schemes to compute the interparticle interaction forces [23, 24]. More details of these improvements are provided below.A review of the literature relevant to the S-C LBM shows that the implementation of EoSs, e.g., Peng-Robinson (P-R) and Carnahan-Starling (CS) EoSs instead of the standard Shan-Chen (SC) EoS, is an attractive and efficient technique to achieve higher density ratios by the elimination of spurious currents. Although, Ezzatneshan and Vaseghnia [18] have recently represented that considering wall boundary conditions with wettability effects impacts the performance of the S-C LBM incorporated with the EoSs. A thermodynamic inconsistency can also appear in the solutions based on the standard S-C model at high-density ratios which leads to a considerable deviation of numerical results from the theoretical data in the phase change diagram [25, 26]. To solve this problem, Gong and Cheng [25] proposed a free parameter that is added into the interaction force computed based on the S-C model to control the strength of interparticle interactions. This technique is employed by Stiles et al. [27] to compute the coexistence curve for a realistic high-density ratio two-phase system, including the air and water. Recently, the method proposed by Gong and Cheng is also implemented for simulation of complex bubble dynamics in a cavitating flow problem [28].In the numerical solution of single-phase fluid flows at high Reynolds numbers using the LBM, the MRT collision operator is introduced as an efficient alternative approach to the standard single relaxation time (SRT) technique [29, 30]. In the same way, the MRT scheme is also implemented in the S-C LBM to eliminate the numerical instability of this model for simulation of two-phase flows with high-viscosity ratios. Kuzmin et al. [21] have investigated the advantages of dealing with the MRT collision operator in the study of capillary effects and the formation of the droplets. Yu and Fan [31] have shown that using the MRT-LBM leads to a considerable reduction in the spurious currents in the interfacial region of two-phase flows. Recently, the MRT S-C LBM is also implemented for studying various flow problems, e.g., the simulation of two-phase flow through porous media [32, 33], investigation of cavitation bubble dynamics [18, 34], and studying the heat transfer in multiphase and multicomponent fluid flows [34, 35]. Evaluation of the performance of the S-C LBM with the implementation of the MRT collision operator for simulation of multiphase flows is an interesting subject which this article deals within different flow conditions.In addition to the aforementioned improvements made for the development of high-density and -viscosity ratio S-C model, some unphysical properties of this model are also resolved by incorporating appropriate forcing schemes. The calculation method of the interparticle interaction force in the pseudopotential LBM has a significant effect on the accuracy and performance of the numerical solution of multiphase flows. A main restriction of the standard S-C LBM in this category is the unphysical dependency of the magnitudes of the density ratio and viscosity ratio of a multiphase system which can be eliminated by the implementation of forcing schemes. There are three well-known forcing schemes that are widely used by the S-C model. The standard Shan-Chen forcing scheme is known as shifted velocity method (SVM) [8, 9], the Guo et al. method [23], and the exact difference method (EDM) [24]. Evaluation of these schemes for the simulation of multiphase flows by using the S-C LBM shows that the EDM is numerically stable in a wide range of density ratios in comparison with the SVM and Guo schemes [16]. It has also been shown that the EDM forcing scheme gives more accurate predictions for high-density ratio multiphase flow characteristics with employing the Carnahan-Starling EoS [36]. Recently, Peng et al. [37] have studied the performance of the SVM, EDM, and Guo schemes for the simulation of the phase separation phenomenon in a two-phase flow system. They found that all these schemes are consistent with the Laplace law; however, the maximum density ratio which can be achieved by the EDM scheme is greater than that of obtained by implementation of other forcing schemes.This wide range of improvements shows that it is crucial to understand the capability and performance of the modified techniques of the S-C LBM, including the EoSs, collision operators, and forcing schemes, to employ an efficient algorithm for simulation of multiphase flow transport through a porous medium. However, to the best knowledge of the authors of this article, there is not a thorough comparative study in the literature to simultaneously evaluate the SRT and MRT collision operators incorporated with appropriate forcing schemes and realistic EoSs for the numerical solution of two-phase flows using the pseudopotential LBM. In this paper, the accuracy and robustness of the S-C LBM are investigated by the implementation of the SRT and MRT collision operators with integrating the SVM and EDM forcing schemes. The capability of the standard SC EoS for simulation of two-phase flows at a high-density ratio is examined in comparison with the CS EoS. In the present study, the method proposed by Gong and Cheng is also used to remove the thermodynamic inconsistency in the simulation of two-phase flow problems studied. The comparative study of these improvements of the S-C LBM is performed by considering an equilibrium state of a droplet suspended in the vapor phase. The numerical technique is validated against the analytical coexistence curve for the liquid-vapor system and the surface tension estimation from the Laplace Law. Then, according to the results obtained, a conclusion has been made to choose an efficient numerical algorithm, including an appropriate collision operator, a realistic EoS, and an accurate forcing scheme, for simulation of multiphase flow transport through the porous medium using the S-C LBM employed. Finally, the characteristics of two-phase flow transport through a porous medium are computed based on the present numerical scheme in different flow conditions defined by the capillary number and the dynamic viscosity ratio. The performance of the implemented numerical method based on the S-C LBM is also investigated by predicting the parameters of a flow transport through the porous medium at high-density and -viscosity ratios.The present paper is organized as follows: In Section2, the LBM with SRT and MRT collision operators is briefly presented. The pseudopotential multiphase LBM coupled with the EoSs and the incorporating of forcing schemes for computing the interparticle interaction forces are introduced in this section. The numerical results obtained for different liquid-vapor flow problems and the relevant discussions are presented in Section 3. Finally, some conclusions are summarized in Section 4.


## 2. Multiphase Lattice Boltzmann Method

In this section, the governing equations of LBM with SRT and MRT collision operators are introduced. Then, the multiphase pseudopotential S-C model [8, 9] with the standard formulation and its improved form with the realistic CS EoS is presented. Two forcing schemes, including SVM and EDM methods, are outlined for the incorporation of interparticle interaction force in the numerical solutions. The modification proposed by Gong and Cheng [25] is also reviewed to eliminate the thermodynamic inconsistency in the present computations.


### 2.1. SRT- and MRT-LBM Equations

The governing equation of the standard single relaxation time (SRT) LBM is given as [38]:
(1)fax+eaDt,t+Dt−fax,t=−Δtτfax,t−faeqx,t,where fax,t is the particle distribution function which demonstrates the distribution of particles with the velocity ea in position x at time t. The lattice speed c=Δx/Δt is the ratio of the lattice space Δx to the time step Δt, and τ is the relaxation time which is related to the kinematic viscosity ν=cs2τ−0.5Δt [3]. The subscript a indicates the discrete velocity directions. According to the two-dimensional (2D) lattice structure D2Q9 employed in the present work, the set of discrete velocities ea is defined as
(2)ea=0,0a=0,ccosa−1π2,sina−1π2a=1−4,2ccos2a−1π4,sin2a−1π4a=5−8.The right-hand side of Eq. (1) represents the SRT collision operator which leads the system to the local Maxwellian equilibrium on the time scale τ. The faeqx,t is the equilibrium distribution function which is defined as
(3)faeqx,t=waρ1+ea.ucs2+ea.u22cs4−u22cs2.The speed of soundcs for the LBM in D2Q9 framework is equal to 1/3, and the weighting coefficient wa is given by
(4)wα=4/9a=0,1/9a=1−4,1/36a=5−8.The macroscopic characteristics of the fluid flow then can be computed as(5)ρ=∑a=1,9fa,u=1r∑a=1,9eafa.For the implementation of MRT collision operator in the LBM presented in Eq. (1), the left-hand side of the equation is rewritten as
(6)fax+eaΔt,t+Δt−fax,t=−M−1Smax,t−max,t−maeqx,t.Herein,ma=fa∗M and maeq=faeq∗M are the distribution function and the equilibrium distribution function in moment space, respectively. The transformation matrix M and the relaxation matrix S can be defined as follows [39]:
(7)M=111111111−4−1−1−1−122224−2−2−2−21111010−101−1−110−20201−1−110010−111−1−100−20211−1−101−11−10000000001−11−1,S=diags1,s2,s3,s4,s5,s6,s7,s8,s9,where the values of s1, s4, and s6 should be nonzero in order to maintain the influence of the trapezoidal integration in Eq. (6) [29]. In the present work, these parameters set to be s1=s4=s6=1, s2=s3=0.8, s5=s7=1.2, and s8=s8=1/τ.


### 2.2. Pseudopotential Shan-Chen Model

Simulation of the phase transition process and the interfacial dynamics in the multiphase flows need to consider the adhesive and cohesive forces between the fluid-fluid and solid-fluid particles. The S-C model is known as an efficient method for computing these interparticle interaction forces in the multiphase LBM framework. Herein, the standard S-C model proposed by Shan and Chen [8, 9] and its improved form by incorporating the CS EoS are presented. The numerical technique proposed by Gong and Cheng [25] is also applied to remove the thermodynamic inconsistency in the present solutions.In the standard S-C LBM, the interaction and adhesive forces are imposed by shifting the macroscopic velocity in the equilibrium distribution functionfaeqx,t as
(8)ueq=u+τΔtFtotalρ,where the total force Ftotal is a summation of the interaction force Fint, adhesive force Fw, and external body force FB:
(9)Ftotal=Fint+Fw+FB.This form of incorporation of the forcing term into the S-C LBM is known as the standard shifted velocity method (SVM) which is discussed in Section2.3.Yuan and Schaefer [17] proposed the following relation for the interaction force Fint:
(10)Fint≅−c0gψ∇ψ,where c0 is a constant parameter and equals to 6 for the D2Q9 lattice structure. Also, the interparticle interaction potential g indicates the attractiveness (repulsiveness) of the interparticle interaction force for g<0 (g>0). In the D2Q9 model, the ∇ψ can be discretized as below:
(11)∂ψi,j∂x=c1ψi+1,j−ψi−1,j+c2ψi+1,j+1+ψi+1,j−1−ψi−1,j+1−ψi−1,j−1,where c1 and c2 are equal to 1/3 and 1/12, respectively. The mean-field potential function ψ defines the relevance between the interaction strength and the mass of particles. In the S-C model, ψ depends on the local density ρ as [8].
(12)ψ=ρ01−exp−ρρ0.The physical velocity of the multiphase flow is given by(13)ureal=u+ΔtFtotal2ρ,and the macroscopic pressure is expressed as
(14)p=cs2ρ+g2cs2ψ2,which is known as the standard SC EoS in this paper. For this EoS, Tc=4.5ρ0 and ρ0=1. To incorporate a realistic EoS in the S-C model, Eq. (14) can be rewritten as
(15)ψ=2p−ρcs2gcs2.Then, the pressure obtained from the EoS in the form ofp=pρ can be used for computing the function ψ. For the simulation of multiphase flows at high-density ratios by using the S-C LBM, a variety of EoSs is proposed in the literature [18]. According to the recent studies in the literature [16, 18], it is shown that the CS EoS provides more realistic results, particularly when the two-phase flow problem includes a wall boundary condition with wettability effects. Also, it is reported that the CS EoS provides more numerical stability by reducing spurious currents in the interfacial region that leads to achieving higher density ratios. Herein, the noncubic CS EoS [40] is implemented to describe the relationship between the fluid properties by the following equation:
(16)p=ρRT1+bρ/4+bρ/42−bρ/431−bρ/43−aρ2.In the CS EoS,a=0.4963R2Tc2/pc and b=0.18727RTc/pc are the attraction and repulsion parameters. The critical temperature is equal to Tc=0.3773a/bR, where R is the gas constant. In the LBM framework, these parameters are set to be a=1, b=4, and R=1 [17]. Accordingly, the critical density ρc and critical temperature Tc will be 0.13045 and 0.0943, respectively.The S-C LBM coupled with any EoSs may lead to the thermodynamic inconsistency in the simulation of multiphase flows at high-density ratios [41]. In order to eliminate such a numerical problem, Gong and Cheng [25] have proposed a modification in computing the interparticle interaction force Fint based on Eq. (10). This equation can be rewritten as follows:
(17)Fint≅−c0g∇ψ2/2.In terms of the numerical implementation, theFint computed by Eq. (10) includes both the local and the nearest-neighbor lattice points, while the present form based on Eq. (17) only deals with the local lattice nodes. In the modification technique proposed by Gong and Cheng, Fint is calculated by a combination of Eqs. (10) and (17) as
(18)Fint≅−βc0gψ∇ψ−1−βc0g∇ψ2/2,where β is a weighting factor that controls the effectiveness of each relation in the calculation of the interaction force. The value of β depends on which EoS is employed in numerical solutions [25]. In this paper, the effect of parameter β is investigated on the accuracy and stability of the present numerical scheme based on the S-C LBM coupled with the CS EoS.Finally, the adhesive forceFw between the two-phase fluid particles and solid nodes is expressed as
(19)Fw=−gψi,jψwc1∑a=1∼4sx+ea+c2∑a=5∼8sx+ea,where ψw is calculated by considering a wall density ρw for solid surface. The desired wettability characteristic of the surface defines the value of ρw. A value of ρw close to the liquid density ρl corresponds to the hydrophilic wall boundary and consequently, the contact angle will be θw<90∘. In contrast, for a ρw value close to the vapor density ρv, the solid surface acts as a hydrophobic surface with θw>90∘. In Eq. (19), the indicator s is used to turn the adhesion force on for the solid nodes by s=1 and turn it off on the fluid nodes by s=0 in the numerical solutions performed.


### 2.3. Forcing Schemes

The total forceFtotal, including the interaction, adhesive, and body forces, is usually implemented via forcing schemes to the pseudopotential S-C LBM for simulation of multiphase flows. Therefore, employing a forcing scheme can affect the performance and numerical stability of the S-C model. In the present study, the standard SVM and the EDM approaches are considered to evaluate the effects of the forcing term incorporation into the S-C LBM.In the SVM forcing scheme with both of the SRT and MRT collision operators, theFtotal is imposed on the velocity field through Eq. (8). Then, the equilibrium distribution function faeq presented in Eq. (3) is computed by the ueq. However, in the EDM approach, the force term is directly added to the right-hand side of Eq. (1) as follows:
(20)fax+eaΔt,t+Δt=fax,t−Δtτfax,t−faeqρ,u+faeqρ,u+Ftotal.Δt/ρ−faeqρ,u.To impose the forcing scheme based on the EDM with MRT collision operator, the force term needs to define in the moment space. Then, Eq. (20) is rewritten for MRT-LBM with EDM scheme as
(21)fax+eaΔt,t+Δt=M−1max,t−Smax,t−maeqρ,u+maeqρ,u+Ftotal.Δt/ρ−maeqρ,u.In the present study, the effect of the implementation of the EDM scheme for the incorporation of the force term in the S-C LBM is investigated for simulation of multiphase flows at high-density and -viscosity ratios. For the sake of clarity, the flowchart of the present solution based on the governing equations is presented in Figure1.Figure 1
Flowchart of S-C LBM implemented in present study.


## 2.1. SRT- and MRT-LBM Equations

The governing equation of the standard single relaxation time (SRT) LBM is given as [38]:
(1)fax+eaDt,t+Dt−fax,t=−Δtτfax,t−faeqx,t,where fax,t is the particle distribution function which demonstrates the distribution of particles with the velocity ea in position x at time t. The lattice speed c=Δx/Δt is the ratio of the lattice space Δx to the time step Δt, and τ is the relaxation time which is related to the kinematic viscosity ν=cs2τ−0.5Δt [3]. The subscript a indicates the discrete velocity directions. According to the two-dimensional (2D) lattice structure D2Q9 employed in the present work, the set of discrete velocities ea is defined as
(2)ea=0,0a=0,ccosa−1π2,sina−1π2a=1−4,2ccos2a−1π4,sin2a−1π4a=5−8.The right-hand side of Eq. (1) represents the SRT collision operator which leads the system to the local Maxwellian equilibrium on the time scale τ. The faeqx,t is the equilibrium distribution function which is defined as
(3)faeqx,t=waρ1+ea.ucs2+ea.u22cs4−u22cs2.The speed of soundcs for the LBM in D2Q9 framework is equal to 1/3, and the weighting coefficient wa is given by
(4)wα=4/9a=0,1/9a=1−4,1/36a=5−8.The macroscopic characteristics of the fluid flow then can be computed as(5)ρ=∑a=1,9fa,u=1r∑a=1,9eafa.For the implementation of MRT collision operator in the LBM presented in Eq. (1), the left-hand side of the equation is rewritten as
(6)fax+eaΔt,t+Δt−fax,t=−M−1Smax,t−max,t−maeqx,t.Herein,ma=fa∗M and maeq=faeq∗M are the distribution function and the equilibrium distribution function in moment space, respectively. The transformation matrix M and the relaxation matrix S can be defined as follows [39]:
(7)M=111111111−4−1−1−1−122224−2−2−2−21111010−101−1−110−20201−1−110010−111−1−100−20211−1−101−11−10000000001−11−1,S=diags1,s2,s3,s4,s5,s6,s7,s8,s9,where the values of s1, s4, and s6 should be nonzero in order to maintain the influence of the trapezoidal integration in Eq. (6) [29]. In the present work, these parameters set to be s1=s4=s6=1, s2=s3=0.8, s5=s7=1.2, and s8=s8=1/τ.


## 2.2. Pseudopotential Shan-Chen Model

Simulation of the phase transition process and the interfacial dynamics in the multiphase flows need to consider the adhesive and cohesive forces between the fluid-fluid and solid-fluid particles. The S-C model is known as an efficient method for computing these interparticle interaction forces in the multiphase LBM framework. Herein, the standard S-C model proposed by Shan and Chen [8, 9] and its improved form by incorporating the CS EoS are presented. The numerical technique proposed by Gong and Cheng [25] is also applied to remove the thermodynamic inconsistency in the present solutions.In the standard S-C LBM, the interaction and adhesive forces are imposed by shifting the macroscopic velocity in the equilibrium distribution functionfaeqx,t as
(8)ueq=u+τΔtFtotalρ,where the total force Ftotal is a summation of the interaction force Fint, adhesive force Fw, and external body force FB:
(9)Ftotal=Fint+Fw+FB.This form of incorporation of the forcing term into the S-C LBM is known as the standard shifted velocity method (SVM) which is discussed in Section2.3.Yuan and Schaefer [17] proposed the following relation for the interaction force Fint:
(10)Fint≅−c0gψ∇ψ,where c0 is a constant parameter and equals to 6 for the D2Q9 lattice structure. Also, the interparticle interaction potential g indicates the attractiveness (repulsiveness) of the interparticle interaction force for g<0 (g>0). In the D2Q9 model, the ∇ψ can be discretized as below:
(11)∂ψi,j∂x=c1ψi+1,j−ψi−1,j+c2ψi+1,j+1+ψi+1,j−1−ψi−1,j+1−ψi−1,j−1,where c1 and c2 are equal to 1/3 and 1/12, respectively. The mean-field potential function ψ defines the relevance between the interaction strength and the mass of particles. In the S-C model, ψ depends on the local density ρ as [8].
(12)ψ=ρ01−exp−ρρ0.The physical velocity of the multiphase flow is given by(13)ureal=u+ΔtFtotal2ρ,and the macroscopic pressure is expressed as
(14)p=cs2ρ+g2cs2ψ2,which is known as the standard SC EoS in this paper. For this EoS, Tc=4.5ρ0 and ρ0=1. To incorporate a realistic EoS in the S-C model, Eq. (14) can be rewritten as
(15)ψ=2p−ρcs2gcs2.Then, the pressure obtained from the EoS in the form ofp=pρ can be used for computing the function ψ. For the simulation of multiphase flows at high-density ratios by using the S-C LBM, a variety of EoSs is proposed in the literature [18]. According to the recent studies in the literature [16, 18], it is shown that the CS EoS provides more realistic results, particularly when the two-phase flow problem includes a wall boundary condition with wettability effects. Also, it is reported that the CS EoS provides more numerical stability by reducing spurious currents in the interfacial region that leads to achieving higher density ratios. Herein, the noncubic CS EoS [40] is implemented to describe the relationship between the fluid properties by the following equation:
(16)p=ρRT1+bρ/4+bρ/42−bρ/431−bρ/43−aρ2.In the CS EoS,a=0.4963R2Tc2/pc and b=0.18727RTc/pc are the attraction and repulsion parameters. The critical temperature is equal to Tc=0.3773a/bR, where R is the gas constant. In the LBM framework, these parameters are set to be a=1, b=4, and R=1 [17]. Accordingly, the critical density ρc and critical temperature Tc will be 0.13045 and 0.0943, respectively.The S-C LBM coupled with any EoSs may lead to the thermodynamic inconsistency in the simulation of multiphase flows at high-density ratios [41]. In order to eliminate such a numerical problem, Gong and Cheng [25] have proposed a modification in computing the interparticle interaction force Fint based on Eq. (10). This equation can be rewritten as follows:
(17)Fint≅−c0g∇ψ2/2.In terms of the numerical implementation, theFint computed by Eq. (10) includes both the local and the nearest-neighbor lattice points, while the present form based on Eq. (17) only deals with the local lattice nodes. In the modification technique proposed by Gong and Cheng, Fint is calculated by a combination of Eqs. (10) and (17) as
(18)Fint≅−βc0gψ∇ψ−1−βc0g∇ψ2/2,where β is a weighting factor that controls the effectiveness of each relation in the calculation of the interaction force. The value of β depends on which EoS is employed in numerical solutions [25]. In this paper, the effect of parameter β is investigated on the accuracy and stability of the present numerical scheme based on the S-C LBM coupled with the CS EoS.Finally, the adhesive forceFw between the two-phase fluid particles and solid nodes is expressed as
(19)Fw=−gψi,jψwc1∑a=1∼4sx+ea+c2∑a=5∼8sx+ea,where ψw is calculated by considering a wall density ρw for solid surface. The desired wettability characteristic of the surface defines the value of ρw. A value of ρw close to the liquid density ρl corresponds to the hydrophilic wall boundary and consequently, the contact angle will be θw<90∘. In contrast, for a ρw value close to the vapor density ρv, the solid surface acts as a hydrophobic surface with θw>90∘. In Eq. (19), the indicator s is used to turn the adhesion force on for the solid nodes by s=1 and turn it off on the fluid nodes by s=0 in the numerical solutions performed.


## 2.3. Forcing Schemes

The total forceFtotal, including the interaction, adhesive, and body forces, is usually implemented via forcing schemes to the pseudopotential S-C LBM for simulation of multiphase flows. Therefore, employing a forcing scheme can affect the performance and numerical stability of the S-C model. In the present study, the standard SVM and the EDM approaches are considered to evaluate the effects of the forcing term incorporation into the S-C LBM.In the SVM forcing scheme with both of the SRT and MRT collision operators, theFtotal is imposed on the velocity field through Eq. (8). Then, the equilibrium distribution function faeq presented in Eq. (3) is computed by the ueq. However, in the EDM approach, the force term is directly added to the right-hand side of Eq. (1) as follows:
(20)fax+eaΔt,t+Δt=fax,t−Δtτfax,t−faeqρ,u+faeqρ,u+Ftotal.Δt/ρ−faeqρ,u.To impose the forcing scheme based on the EDM with MRT collision operator, the force term needs to define in the moment space. Then, Eq. (20) is rewritten for MRT-LBM with EDM scheme as
(21)fax+eaΔt,t+Δt=M−1max,t−Smax,t−maeqρ,u+maeqρ,u+Ftotal.Δt/ρ−maeqρ,u.In the present study, the effect of the implementation of the EDM scheme for the incorporation of the force term in the S-C LBM is investigated for simulation of multiphase flows at high-density and -viscosity ratios. For the sake of clarity, the flowchart of the present solution based on the governing equations is presented in Figure1.Figure 1
Flowchart of S-C LBM implemented in present study.


## 3. Results and Discussions

In this section, the numerical algorithms presented based on the standard and modified S-C LBM are evaluated for the simulation of multiphase flows. For this comparison study, a combination of numerical schemes with SRT and MRT collision operators, SC and CS EoSs, SVM and EDM forcing schemes, and also the effect of control parameterβ is considered. At first, the accuracy and performance of different algorithms applied based on the S-C LBM are investigated by the simulation of two benchmark two-phase flow problems, including the equilibrium state of a droplet suspended in the vapor phase and equilibrium state of a droplet on a solid wall with wettability effects. The results obtained based on the present solutions are compared with the analytical and available numerical results. According to this comparison study, a numerical scheme is proposed for simulation of two-phase flow problems at high-density and –viscosity ratios. Then, the efficiency of the proposed solution technique is validated by the simulation of two-phase flow transport through a duct. Finally, the numerical solution of the two-phase flow transport through a porous medium is performed by employing the proposed S-C LBM, and the flow characteristics are studied in various flow conditions.


### 3.1. Study of Accuracy and Performance

For a 2D liquid droplet suspended in the vapor phase, the surface tensionσ causes a pressure difference Δp between the outside and inside of the droplet. The Laplace law [42] defines a linear relation between σ and Δp for a certain droplet radius R as
(22)Δp=pin−pout=σR.In order to verify the present numerical results obtained based on the S-C LBM by the Laplace law, a stationary droplet placed in a fully periodic square domain is considered, and simulations are performed for various droplet radiuses. Herein, the droplet is initialized in the center of the flow fieldx0,y0, and the density distribution is defined by Eq. (23) to preserve the stability of numerical solutions at the beginning iterations [12].
(23)ρ=ρl+ρv2−ρl−ρv2tanh2x−x02+y−y02−r0w.Figure2 illustrates the results obtained for Δp versus 1/R by using the SRT-SVM, SRT-EDM, and MRT-EDM schemes coupled with the SC EoS. The linear relation between these two parameters in this figure shows that the present numerical results based on all the mentioned schemes satisfy the Laplace law. The slope of the line represents the value of the surface tension of the liquid-vapor system studied which is obtained σ=0.05923, 0.06442, and 0.06292 by the SRT-SVM, SRT-EDM, and MRT-EDM schemes, respectively. The satisfaction of the Laplace law in this study demonstrates the accuracy of the employed S-C LBM coupled with different forcing schemes for predicting the characteristics of a stationary two-phase flow system. The slight difference between the values obtained for the surface tension depends on the efficiency and accuracy of these schemes which is discussed below with more details.Figure 2
Verification of Laplace law for simulation of 2D droplet by employing different S-C LBM schemes coupled with SC EoS.The accuracy of the present numerical schemes based on the S-C LBM is evaluated by using SRT and MRT collision operators, SC and CS EoSs, and SVM and EDM forcing schemes with the simulation of the stationary droplet at high-density and –viscosity ratios. Figure3 demonstrates the density of the coexisting phases calculated in the equilibrium state of the droplet at different values of reduced temperature T/Tc with relaxation time τ=1.0. It should be noted that the parameter g defines the temperature effects in the SC EoS by g=−1/T, which gives T/Tc=−4.5ρ0/g. Herein, the analytical solution based on the theory of Maxwell equal-area construction [43] is used to verify the present coexistence curves. Also, to consider the effect of the β-modification to remove the thermodynamic inconsistency in computing the interaction force by Eq. (18), this parameter is set to be β=1 for the standard calculations and β≠1 for the modified S-C LBM.Figure 3
Coexistence curves obtained by present S-C LBM with SC EoS (a, b) and CS EoS (c, d) in comparison with analytical solutions atτ=1 for different forcing schemes and β parameters.
(a)(b)(c)(d)The numerical results presented in Figure3 indicate that all of the S-C LBM schemes applied based on both the SC and CS EOSs accurately predicted the density of the liquid phase in comparison with the analytical solution. Also, at higher reduced temperatures correspond to the low-density ratios, the numerical results for the vapor phase are in good agreement with the analytical solution. However, the density values computed for the vapor phase at higher density ratios (low reduced temperatures) have a significant deviation from the analytical solution by all the schemes coupled with the standard interparticle interaction force (β=1). The x−axis of Figure 3 is plotted in a log scale to have a better view of the difference between the results. It is obvious that below the reduced temperature T/Tc=0.85, the deviation between the numerical results and the analytical solution becomes substantial. Due to the higher value of the vapor density predicted by the SC EOS, the difference between the vapor density of this model with the analytical solution does not a significant impact on the value of the density ratio of the whole two-phase flow system. However, for the CS EOS, a slight deviation of the vapor density from the analytical solution leads to a considerable variation of the flow density ratio because of very small values computed for the density of the vapor phase by this EOS. Therefore, it is crucial to improve the consistency of these EOSs with the analytical solutions for accurate simulation of two-phase flow systems at high-density ratios. As mentioned before, the β-modification is considered in the present work to resolve this problem.As shown in Figure3, the results confirm that the improved interaction forcing scheme by the β parameter has a significant effect on the accuracy of the SC and CS EOSs. In this case, the appropriate value obtained for β by CS EOS is between 1.15 and 1.18 depending on the SRT and MRT collision operators employed. Although, there is no considerable difference between the numerical results obtained based on these two β values. For the SC EoS, the value of β=0.886 is employed. Using different values for β according to the type of EOSs is also reported by Ping and Cheng [25]. By employing the β-modification, the density values obtained by the present numerical simulations for the vapor phase is in good agreement with those obtained by the analytical solution in a wide range of reduced temperatures. Although, the results obtained by the SC EOS still demonstrate some deviation below the T/Tc≃0.55 which can be due to the strong spurious currents in the interfacial region at such density ratios [44]. It is worth to mention that the SRT and MRT collision operators have not significant effect on the accuracy of the results obtained for the density ratio in comparison with the analytical solution. However, the MRT collision operator influences the numerical stability by decreasing the spurious currents at low reduced temperatures.For more detail discussions regarding to the unphysical currents, the magnitude of maximum spurious velocityUsmax obtained by the standard and modified interaction forcing schemes based on the SC and CS EOSs is presented in Figure 4 for the equilibrium state of a droplet at various reduced temperatures. This figure indicates that with employing both the SC and CS EoSs, the SRT-SVM and SRT-EDM models lead to almost the same values for Usmax with the implementation of both the standard and improved interaction forces. However, with applying the MRT collision operator, the computed values for Usmax are lower than those of obtained by SRT collision operator for all the schemes at various reduced temperatures. Increment of the magnitude of the maximum spurious velocity with decreasing the reduced temperature is obvious in Figure 4. The results presented in this figure also show that by applying CS EoS for simulation of the droplet at high reduced temperatures (low-density ratios), the parameter β has no remarkable effect on the magnitude of spurious currents. In contrast, at high-density ratios, the Usmax is significantly decreased by using the improved interaction force, particularly with employing the MRT-EDM model. This result demonstrates that the CS EoS is more realistic and efficient than SC EoS for simulation of two-phase flows at high-density ratios due to producing less spurious currents at the interfacial region.Figure 4
Spurious currents obtained by present S-C LBM with SC EoS (a) and CS EoS (b) atτ=1 for different forcing schemes and β parameters.
(a)(b)As described in relationν=cs2τ−0.5, the viscosity value is determined by the value of relaxation time τ in numerical solutions based on the LBM. Then, τ⟶0.5 corresponds to lower viscosity magnitudes which can be used for simulation of two-phase flows at high Reynolds numbers. However, in flow simulations based on the standard LBM, serious numerical instabilities arise at the values of relaxation time close to 0.5 [45]. It is expected that this numerical limitation will be intensified when the two-phase flow simulations are performed at high-density ratios. To study the mentioned restriction and examine the capability of the present numerical method based on the S-C LBM with SRT and MRT collision operators, SC and CS EoSs and SVM and EDM forcing schemes for the simulation of two-phase flows at low viscosity magnitudes, the simulation of the stationary droplet is considered at τ=0.52. This value for τ is the minimum relaxation time that can be achieved a stable numerical solution of the stationary droplet at reduced temperature T/Tc=0.5 by the present MRT-EDM model and CS EoS.Figure5 presents the coexistence curve for the stationary droplet simulations at τ=0.52 by both the SC and CS EOSs and implementation of the numerical schemes employed. At the low viscosity magnitude considered (τ=0.52), it is obvious in Figure 5 that the MRT-EDM scheme is more efficient than SRT-SVM and SRT-EDM techniques to achieve lower reduced temperatures that lead to the simulation of the two-phase flow at high-density ratios. The obtained results show that the minimum achievable reduced temperature is T/Tc=0.6 by applying the SC EOS with the MRT-EDM scheme. However, by employing the CS EOS with MRT-EDM and modified interaction force, this magnitude is T/Tc=0.475, that consequents a density ratio of 1597.8. To study the numerical stability of the present methods, Figure 6 illustrates the spurious currents calculated in terms of the reduced temperature. It can be seen that the MRT-EDM with the improved forcing scheme produces minimum spurious currents among the other methods. By considering the results presented in Figures 5 and 6, the advantages of employing the MRT scheme with β− modification of the forcing scheme are obvious for the simulation of a two-phase flow system at low relaxation times (τ⟶0.5).Figure 5
Coexistence curves obtained by present S-C LBM with SC EoS (a) and CS EoS (b) in comparison with analytical solutions atτ=0.52 for different forcing schemes and β parameters.
(a)(b)Figure 6
Spurious currents obtained by present S-C LBM with SC EoS (a) and CS EoS (b) atτ=0.52 for different forcing schemes and β parameters.
(a)(b)It should be noted that using the EDM scheme in the present study eliminates the unphysical dependency of the magnitudes of the density ratio and viscosity ratio of the two-phase flow system considered. To demonstrate this feature, simulation of a droplet in the equilibrium state is performed by the S-C LBM with SRT-SVM, SRT-EDM, and MRT-EDM schemes atT/Tc=0.8. Figure 7 shows the variation of the density ratio of each model at different relaxation times. As can be seen, by using the SVM scheme, the density ratio is decreased by increasing the magnitude of the relaxation time at a defined reduced temperature, which is unphysical. However, by using the EDM scheme, the density ratio remains constant as the reduced temperature T/Tc changes. The notable difference between the MRT-EDM and SRT-EDM schemes is the minimum achievable value of the relaxation time for the simulation of the two-phase flow by using these techniques. As given in Table 1, the present study shows that for simulation of a droplet at the equilibrium state, the present numerical solution is still stable at τ=0.5001 (very low viscosity value) by using the MRT-EDM scheme, which is more less than other two methods considered.Figure 7
Comparison of dependency between density ratio and relaxation time by different forcing schemes of S-C LBM for simulation of a 2D droplet.Table 1
Minimum achievable relaxation time by using different forcing schemes for simulation of a 2D droplet.
ModelsτυminSRT-SVM0.60.034SRT-EDM0.5150.005MRT-EDM0.50013.34×10−6


### 3.2. Model Validation

Predicting the interfacial dynamics by considering wetting surfaces in a liquid-vapor system is essential for simulation of multiphase flow transport through porous media. In order to demonstrate the capability of the S-C LBM employed for the numerical solution of such a flow problem, the simulation results of the equilibrium state of a droplet on a solid wall and a two-phase flow through a channel are verified with considering wettability effects. According to the accuracy and performance study of the different techniques considered in the previous section, the MRT-EDM scheme is used in this validation study by applying the CS EoS and employing the modified interaction force with theβ parameter.The simulation of the equilibrium state of a droplet on a solid wall is performed in a square domain with the grid size300×300, in which the periodic condition is considered in the left and right sides, and the wall boundary condition is imposed on the top and bottom of the flowfield. The study is carried out for a droplet placed on the bottom wall with radius 45lu at different wetting conditions. To verify the numerical scheme employed, the contact angle θw between the droplet interface and the solid surface is computed for different normalized density Dρ which is defined by
(24)Dρ=ρw−ρυρl−ρυ.In Figure8, the results obtained based on the present solution algorithm are compared with those reported in Ref. [12]. As shown in this figure, the wall surface becomes hydrophilic by increasing the ρw, which leads to the lower contact angle θw, and vice versa. It should be noted that the linear dependency between the value of Dρ and the computed contact angle is also reported in the literature [20, 46, 47] which the present results show good agreement.Figure 8
Comparison of wetting angle of a droplet on surface obtained by employing SVM-MRT scheme.The second test case for model validation is the simulation of the two-phase fluid flow through a duct which can be considered as a simplified pore of a porous medium. The transmitting of a two-phase flow within a channel causes the wetting phase to flow near the hydrophilic wall which forms a film of wetting fluid on the surface. Consequently, the nonwetting phase is forced to flow through the central core of the channel. The schematic of the two-phase flow inside a 2D channel is demonstrated in Figure9. The parameters a and b in this figure indicate the height of the nonwetting phase and the height of the channel from the centerline, respectively. Herein, the periodic boundary condition is considered in the streamwise direction, and the bounce-back boundary condition is imposed on the solid walls.Figure 9
Geometry and boundary conditions for simulation of two-phase flow through a 2D channel.For this two-phase flow problem, the fully developed flow condition is obtained by applying an external body forceF. Therefore, the analytical solution of the fully developed flow velocity for the two-phase flow in the cross-section of the channel is defined as follows [46]:
(25)ux=F2υwρwb2−y2a<y<bux=F2υwρwb2−y2+F2υnwρnwa2−y20<y<a,where υw, ρw, υnw, and ρnw are the kinematic viscosity and density of the wetting and nonwetting phases. According to Eq. (25), the velocity profile of the fully developed two-phase flow in the channel is related to the dynamic viscosity ratio M=μnw/μw and the saturation of each phase [46], which can be defined as
(26)Sw=b−a/b,Snw=a/b.In the present study, the fully developed velocity profile of the two-phase flow through the channel obtained based on the numerical solution is compared with the analytical data in Figure10 at the flow conditions sw=0.33 and M≈0.1,0.32. The velocity profile of the nonwetting phase predicted by the present MRT-EDM scheme has about 9.5% and 6% deviation from the analytical profile at M=0.1 and M=0.32, respectively. Huang and Lu [48] have reported about 7% difference between their numerical solution and the analytical data for the velocity profile of the two-phase flow through the channel. This comparison shows that the present numerical technique is accurate enough for the simulation of the two-phase flow pattern in 2D pore scales.Figure 10
Distribution of streamwise velocity componentux in a cross-section in the middle of 2D channel with flow condition F=1.5×10−8 at Sw=0.333 and (a) M=0.1 and (b) M=0.32.
(a)(b)To assess the capability of the implemented S-C LBM for computing the two-phase flow characteristics in the pore scale, the relative permeability of the wetting and nonwetting phases through the channel is computed at different flow conditions. Herein, the definition proposed by Yiotis et al. [46] is used for calculating the relative permeability of each phase as
(27)kr,w=12Sw23−Sw,kr,nw=Snw32M+Snw21−32M,where kr,w and kr,nw denote relative permeability of wetting and nonwetting phases, respectively. These relations define that in a 2D channel, kr,w only depends on the saturation of the wetting phase, while kr,nw is related to the saturation of the nonwetting phase and the dynamic viscosity ratio. Figure 11 illustrates that the relative permeability obtained based on the present numerical scheme is in excellent agreement with those computed with the analytical relations at different wetting phase saturations. A similar study is performed at different relative temperatures T/Tc to examine the accuracy of the S-C LBM for predicting the relative permeability magnitude in various density ratios. Figure 12 shows the comparison of the relative permeability computed by the present numerical scheme with those obtained by the analytical solution as a function of the reduced temperature. This study shows a very good agreement between the numerical and analytical solutions with the implementation of the modified interaction forcing scheme with the parameter β=1.18, while using the standard interaction force (β=1) causes a serious deviation between the results at the flow conditions with T/Tc<0.7, which corresponds to the high-density ratios. Therefore, the choice of the interaction forcing scheme becomes important for such flow conditions. This validation study shows that the S-C LBM with employing the MRT-EDM scheme, CS EoS, and the modified forcing scheme with parameter β is efficient and accurate for estimation of the two-phase flow structures, and characteristics through the channel and can be used for simulation of such flow problems in porous media.Figure 11
Comparison of relative permeability of the wetting and nonwetting phases obtained based on present numerical scheme with analytical solution atM=0.1.Figure 12
Comparison of relative permeability of wetting and nonwetting phases at various reduced temperatures for both standard and modified interaction forces.


### 3.3. Two-Phase Flow Transport through a Porous Medium

The present numerical scheme developed based on the S-C LBM with employing the CS EoS, MRT-EDM, and modified forcing techniques is applied for studying the two-phase flow transport through a porous medium. The capillary numberCa and the dynamic viscosity ratio M are considered as two effective parameters that characterize the two-phase flow structure in this study [49]. Herein, the capillary number is defined by
(28)Ca=unwμnwσcosθ,where μnw and unw are the dynamic viscosity and velocity of the advancing nonwetting fluid, respectively, and θ is the fluid-fluid contact angle. Herein, a homogeneous porous medium is generated by a distribution of obstacles in a flow domain with 1000×300 grid points. The pore space area and the height of the throat between the obstacles include 20×20 lattice and 4∼7 grid points, respectively. These considered grid points are required for the spatial discretization of a pore to perform the solution accurately enough for capturing the interfacial dynamics in the pore scale. Hence, the existence of very small pores in the geometry of a porous medium needs a high-resolution grid size and increases the computational cost. The geometry of the porous medium and the boundary conditions of the flow domain are shown in Figure 13. The velocity inlet and pressure outlet boundary conditions are imposed in the left and right boundaries, respectively, and the top and bottom sides are considered periodic. The wetting fluid is defined by the higher density (red color) which invades the nonwetting fluid with lower density (blue color) in the flow domain. Such an invasion produces some fingering regimes that their origin can be the surface or geometrical characteristics of a porous medium in nature, e.g., the nonuniform surface wettability of a sandstone [50]. Herein, the fingering regime during the two-phase flow simulations through the considered porous medium is occurred due to the random generation of the shape of obstacles and the wetting condition of the solid surfaces.Figure 13
Geometry and boundary conditions of porous medium.Figure14 demonstrates the two-phase flow patterns obtained by the present numerical scheme based on the S-C LBM at the dynamic viscosity ratio M=1017 with different capillary numbers. To define the capillary number, the surface tension is calculated by using the Laplace law analysis, which is obtained σ=0.02316. Also, the fluid-fluid contact angle is measured about 175∘ which shows fully hydrophobicity of the invading fluid. This study is performed at the reduced temperature T/Tc=0.482, and the parameter β is set to be 1.241 to preserve the thermodynamic consistency of the solutions. From the results presented in Figure 14, it can be seen that the penetration of the invading fluid at a certain time t=50000 is decreased by reducing the capillary number. These patterns are discussed with more details along with the explanations of Figures 15 and 16 at the end of this section.Figure 14
Two-phase flow patterns in the porous medium obtained by present simulations for flow conditionlogM=3.07 with various capillary numbers at t=50000.Figure 15
Results obtained for displacement pattern of the invading fluid in porous medium by the present S-C LBM for flow conditionlogM=3.07 with various capillary numbers in comparison with the schematic of phase diagram presented by Zhang et al. [51] and Lenormand et al. [49].Figure 16
Instantaneous displacement pattern of two-phase flow through the porous medium obtained by present simulations for flow conditionlogM=3.07 and logCa=−1.85 at three different times.To compare the obtained results for the two-phase flow characteristics through the porous medium in a similar situation, a normalized timeT∗ is defined by division of the time step to the total simulation time of each case. The total simulation time is counted until the first finger of the invading fluid touches the outer side of the porous medium for each case, which is called breakthrough [49]. Figure 17 illustrates the variation of the pressure difference Δp between the inlet and outlet of the porous medium versus the normalized time T∗ at different capillary numbers. As seen in this figure, by increment of the capillary number in a constant dynamic viscosity ratio, the invading fluid has enough pressure difference to overcome the throat threshold pressure pc presented by Lenormand et al. [49], which is known as a resistance factor. The throat threshold pressure is defined by the following relation:
(29)pc=2σcosθr,where σ is the surface tension, θ is the fluid-fluid contact angle, and r is the throat diameter. This parameter also affects the nonwetting phase saturation. In the present study, by considering the minimum and maximum throat radius, the threshold pressure is altered between ‐0.006591 and ‐0.011535 for r=7 and r=4, respectively. The negative sign indicates the resistance pressure against the invading fluid.Figure 17
Predicted pressure differenceΔp between the two sides of the porous medium in streamwise direction versus normalized time T∗ for flow condition logM=3.07 with various capillary numbers.Figure18 shows the nonwetting phase saturation as a function of the normalized time T∗. The results obtained demonstrate that the saturation of the nonwetting fluid is enhanced by increasing the capillary number, which is consistent with the previous works [49–51]. The final values of the nonwetting phase saturation and the pressure difference at the breakthrough condition are illustrated in Figure 19 at different capillary numbers. It is obvious that the variation of the nonwetting phase saturation is almost linear which is in agreement with the results obtained by Zhang et al. [51] and Liu et al. [50]. The agreement of the present results in comparison with those reported in the literature confirms that the implemented S-C LBM with the MRT-EDM scheme besides the modified interaction forcing method and the CS EoS is capable for simulation of two-phase flows at high-density and -viscosity ratio through porous media at different flow conditions.Figure 18
Nonwetting phase saturation versus normalized timeT∗ for flow condition logM=3.07 with various capillary numbers.Figure 19
Comparison of nonwetting fluid saturation and pressure differenceΔp at the breakthrough time for flow condition logM=3.07 with various capillary numbers.The structure of the two-phase flow through the porous medium is also investigated at different flow conditions. The displacement pattern of the invading fluid is categorized as the capillary fingering, viscous fingering, and stable displacement, depending on the capillary number and the viscosity ratio between two phases [52]. A phase diagram is used to describe these flow structures by logCa and logM for a gas-water flow system (see Figure 15). In the present simulations based on the S-C LBM, the magnitude of the viscosity ratio is set as logM>0. Therefore, the capillary fingering and stable displacement patterns are expected in the porous medium considered. As can be seen in Figure 14, the two-phase flow pattern is altered from the stable displacement to the capillary fingering by a decrement of the capillary number. Figure 16 demonstrates the present results for instantaneous two-phase flow displacement at different times for logCa=−1.85, in which the pattern behavior is almost stable displacement. In this flow condition, the interface of the invading fluid penetrates into the porous medium along a straight line, despite some irregularities due to the pore-scale geometry. This behavior can be the result of the dominant effect of the viscous force and the higher pressure difference that cause to fill the porous medium uniformly and simultaneously. By decreasing the capillary number, it is observed that the invading fluid is extended in the shape of irregular fingers in the flow field due to the dominant effects of the capillary force in comparison with the viscous force. Such a capillary fingering of two-phase flow through the porous medium is shown in Figure 20 at logCa=−3.5. In this case, due to the uniform pressure distribution throughout the domain, the nonwetting fluid only penetrates the pore-scale area with a lower throat pressure threshold. As a result, the direction of the flow displacement in the porous medium is irregular. In Figure 20, it can be seen that the capillary fingers of the nonwetting fluid form several loops that trap some slugs of the wetting fluid among the porous medium. These two-phase flow structures obtained using the present numerical scheme are in consistent with those reported in the literature [49, 50].Figure 20
Instantaneous displacement pattern of two-phase flow through the porous medium obtained by present simulations for flow conditionlogM=3.07 and logCa=−3.5 at various times.The interfacial area between two phases is quantified for the two-phase flow through the porous medium at different capillary numbers. This parameter is an important factor that influences the mass and energy transfer among the phases [50, 51]. In the present 2D studies, the interfacial length is measured instead of the interfacial area. Herein, due to the extreme hydrophilicity of the wetting fluid, the length of the wetting films around the obstacles is also added to the length of the measured interface between the phases for computing the interfacial length. To have comparable results for this parameter, a specific interfacial length is calculated as the ratio of the measured interfacial length to the pore area of the porous medium. Figure 21 shows the variation of the specific interfacial length as a function of the nonwetting phase saturation at different capillary numbers. The present results demonstrate that there is a linear relation between two parameters in all the capillary numbers considered which is in good agreement with the previous experimental finding of Zhang et al. [51].Figure 21
Comparison of specific interfacial length as a function of nonwetting phase saturation at different capillary numbers for flow conditionlogM=3.07.


## 3.1. Study of Accuracy and Performance

For a 2D liquid droplet suspended in the vapor phase, the surface tensionσ causes a pressure difference Δp between the outside and inside of the droplet. The Laplace law [42] defines a linear relation between σ and Δp for a certain droplet radius R as
(22)Δp=pin−pout=σR.In order to verify the present numerical results obtained based on the S-C LBM by the Laplace law, a stationary droplet placed in a fully periodic square domain is considered, and simulations are performed for various droplet radiuses. Herein, the droplet is initialized in the center of the flow fieldx0,y0, and the density distribution is defined by Eq. (23) to preserve the stability of numerical solutions at the beginning iterations [12].
(23)ρ=ρl+ρv2−ρl−ρv2tanh2x−x02+y−y02−r0w.Figure2 illustrates the results obtained for Δp versus 1/R by using the SRT-SVM, SRT-EDM, and MRT-EDM schemes coupled with the SC EoS. The linear relation between these two parameters in this figure shows that the present numerical results based on all the mentioned schemes satisfy the Laplace law. The slope of the line represents the value of the surface tension of the liquid-vapor system studied which is obtained σ=0.05923, 0.06442, and 0.06292 by the SRT-SVM, SRT-EDM, and MRT-EDM schemes, respectively. The satisfaction of the Laplace law in this study demonstrates the accuracy of the employed S-C LBM coupled with different forcing schemes for predicting the characteristics of a stationary two-phase flow system. The slight difference between the values obtained for the surface tension depends on the efficiency and accuracy of these schemes which is discussed below with more details.Figure 2
Verification of Laplace law for simulation of 2D droplet by employing different S-C LBM schemes coupled with SC EoS.The accuracy of the present numerical schemes based on the S-C LBM is evaluated by using SRT and MRT collision operators, SC and CS EoSs, and SVM and EDM forcing schemes with the simulation of the stationary droplet at high-density and –viscosity ratios. Figure3 demonstrates the density of the coexisting phases calculated in the equilibrium state of the droplet at different values of reduced temperature T/Tc with relaxation time τ=1.0. It should be noted that the parameter g defines the temperature effects in the SC EoS by g=−1/T, which gives T/Tc=−4.5ρ0/g. Herein, the analytical solution based on the theory of Maxwell equal-area construction [43] is used to verify the present coexistence curves. Also, to consider the effect of the β-modification to remove the thermodynamic inconsistency in computing the interaction force by Eq. (18), this parameter is set to be β=1 for the standard calculations and β≠1 for the modified S-C LBM.Figure 3
Coexistence curves obtained by present S-C LBM with SC EoS (a, b) and CS EoS (c, d) in comparison with analytical solutions atτ=1 for different forcing schemes and β parameters.
(a)(b)(c)(d)The numerical results presented in Figure3 indicate that all of the S-C LBM schemes applied based on both the SC and CS EOSs accurately predicted the density of the liquid phase in comparison with the analytical solution. Also, at higher reduced temperatures correspond to the low-density ratios, the numerical results for the vapor phase are in good agreement with the analytical solution. However, the density values computed for the vapor phase at higher density ratios (low reduced temperatures) have a significant deviation from the analytical solution by all the schemes coupled with the standard interparticle interaction force (β=1). The x−axis of Figure 3 is plotted in a log scale to have a better view of the difference between the results. It is obvious that below the reduced temperature T/Tc=0.85, the deviation between the numerical results and the analytical solution becomes substantial. Due to the higher value of the vapor density predicted by the SC EOS, the difference between the vapor density of this model with the analytical solution does not a significant impact on the value of the density ratio of the whole two-phase flow system. However, for the CS EOS, a slight deviation of the vapor density from the analytical solution leads to a considerable variation of the flow density ratio because of very small values computed for the density of the vapor phase by this EOS. Therefore, it is crucial to improve the consistency of these EOSs with the analytical solutions for accurate simulation of two-phase flow systems at high-density ratios. As mentioned before, the β-modification is considered in the present work to resolve this problem.As shown in Figure3, the results confirm that the improved interaction forcing scheme by the β parameter has a significant effect on the accuracy of the SC and CS EOSs. In this case, the appropriate value obtained for β by CS EOS is between 1.15 and 1.18 depending on the SRT and MRT collision operators employed. Although, there is no considerable difference between the numerical results obtained based on these two β values. For the SC EoS, the value of β=0.886 is employed. Using different values for β according to the type of EOSs is also reported by Ping and Cheng [25]. By employing the β-modification, the density values obtained by the present numerical simulations for the vapor phase is in good agreement with those obtained by the analytical solution in a wide range of reduced temperatures. Although, the results obtained by the SC EOS still demonstrate some deviation below the T/Tc≃0.55 which can be due to the strong spurious currents in the interfacial region at such density ratios [44]. It is worth to mention that the SRT and MRT collision operators have not significant effect on the accuracy of the results obtained for the density ratio in comparison with the analytical solution. However, the MRT collision operator influences the numerical stability by decreasing the spurious currents at low reduced temperatures.For more detail discussions regarding to the unphysical currents, the magnitude of maximum spurious velocityUsmax obtained by the standard and modified interaction forcing schemes based on the SC and CS EOSs is presented in Figure 4 for the equilibrium state of a droplet at various reduced temperatures. This figure indicates that with employing both the SC and CS EoSs, the SRT-SVM and SRT-EDM models lead to almost the same values for Usmax with the implementation of both the standard and improved interaction forces. However, with applying the MRT collision operator, the computed values for Usmax are lower than those of obtained by SRT collision operator for all the schemes at various reduced temperatures. Increment of the magnitude of the maximum spurious velocity with decreasing the reduced temperature is obvious in Figure 4. The results presented in this figure also show that by applying CS EoS for simulation of the droplet at high reduced temperatures (low-density ratios), the parameter β has no remarkable effect on the magnitude of spurious currents. In contrast, at high-density ratios, the Usmax is significantly decreased by using the improved interaction force, particularly with employing the MRT-EDM model. This result demonstrates that the CS EoS is more realistic and efficient than SC EoS for simulation of two-phase flows at high-density ratios due to producing less spurious currents at the interfacial region.Figure 4
Spurious currents obtained by present S-C LBM with SC EoS (a) and CS EoS (b) atτ=1 for different forcing schemes and β parameters.
(a)(b)As described in relationν=cs2τ−0.5, the viscosity value is determined by the value of relaxation time τ in numerical solutions based on the LBM. Then, τ⟶0.5 corresponds to lower viscosity magnitudes which can be used for simulation of two-phase flows at high Reynolds numbers. However, in flow simulations based on the standard LBM, serious numerical instabilities arise at the values of relaxation time close to 0.5 [45]. It is expected that this numerical limitation will be intensified when the two-phase flow simulations are performed at high-density ratios. To study the mentioned restriction and examine the capability of the present numerical method based on the S-C LBM with SRT and MRT collision operators, SC and CS EoSs and SVM and EDM forcing schemes for the simulation of two-phase flows at low viscosity magnitudes, the simulation of the stationary droplet is considered at τ=0.52. This value for τ is the minimum relaxation time that can be achieved a stable numerical solution of the stationary droplet at reduced temperature T/Tc=0.5 by the present MRT-EDM model and CS EoS.Figure5 presents the coexistence curve for the stationary droplet simulations at τ=0.52 by both the SC and CS EOSs and implementation of the numerical schemes employed. At the low viscosity magnitude considered (τ=0.52), it is obvious in Figure 5 that the MRT-EDM scheme is more efficient than SRT-SVM and SRT-EDM techniques to achieve lower reduced temperatures that lead to the simulation of the two-phase flow at high-density ratios. The obtained results show that the minimum achievable reduced temperature is T/Tc=0.6 by applying the SC EOS with the MRT-EDM scheme. However, by employing the CS EOS with MRT-EDM and modified interaction force, this magnitude is T/Tc=0.475, that consequents a density ratio of 1597.8. To study the numerical stability of the present methods, Figure 6 illustrates the spurious currents calculated in terms of the reduced temperature. It can be seen that the MRT-EDM with the improved forcing scheme produces minimum spurious currents among the other methods. By considering the results presented in Figures 5 and 6, the advantages of employing the MRT scheme with β− modification of the forcing scheme are obvious for the simulation of a two-phase flow system at low relaxation times (τ⟶0.5).Figure 5
Coexistence curves obtained by present S-C LBM with SC EoS (a) and CS EoS (b) in comparison with analytical solutions atτ=0.52 for different forcing schemes and β parameters.
(a)(b)Figure 6
Spurious currents obtained by present S-C LBM with SC EoS (a) and CS EoS (b) atτ=0.52 for different forcing schemes and β parameters.
(a)(b)It should be noted that using the EDM scheme in the present study eliminates the unphysical dependency of the magnitudes of the density ratio and viscosity ratio of the two-phase flow system considered. To demonstrate this feature, simulation of a droplet in the equilibrium state is performed by the S-C LBM with SRT-SVM, SRT-EDM, and MRT-EDM schemes atT/Tc=0.8. Figure 7 shows the variation of the density ratio of each model at different relaxation times. As can be seen, by using the SVM scheme, the density ratio is decreased by increasing the magnitude of the relaxation time at a defined reduced temperature, which is unphysical. However, by using the EDM scheme, the density ratio remains constant as the reduced temperature T/Tc changes. The notable difference between the MRT-EDM and SRT-EDM schemes is the minimum achievable value of the relaxation time for the simulation of the two-phase flow by using these techniques. As given in Table 1, the present study shows that for simulation of a droplet at the equilibrium state, the present numerical solution is still stable at τ=0.5001 (very low viscosity value) by using the MRT-EDM scheme, which is more less than other two methods considered.Figure 7
Comparison of dependency between density ratio and relaxation time by different forcing schemes of S-C LBM for simulation of a 2D droplet.Table 1
Minimum achievable relaxation time by using different forcing schemes for simulation of a 2D droplet.
ModelsτυminSRT-SVM0.60.034SRT-EDM0.5150.005MRT-EDM0.50013.34×10−6


## 3.2. Model Validation

Predicting the interfacial dynamics by considering wetting surfaces in a liquid-vapor system is essential for simulation of multiphase flow transport through porous media. In order to demonstrate the capability of the S-C LBM employed for the numerical solution of such a flow problem, the simulation results of the equilibrium state of a droplet on a solid wall and a two-phase flow through a channel are verified with considering wettability effects. According to the accuracy and performance study of the different techniques considered in the previous section, the MRT-EDM scheme is used in this validation study by applying the CS EoS and employing the modified interaction force with theβ parameter.The simulation of the equilibrium state of a droplet on a solid wall is performed in a square domain with the grid size300×300, in which the periodic condition is considered in the left and right sides, and the wall boundary condition is imposed on the top and bottom of the flowfield. The study is carried out for a droplet placed on the bottom wall with radius 45lu at different wetting conditions. To verify the numerical scheme employed, the contact angle θw between the droplet interface and the solid surface is computed for different normalized density Dρ which is defined by
(24)Dρ=ρw−ρυρl−ρυ.In Figure8, the results obtained based on the present solution algorithm are compared with those reported in Ref. [12]. As shown in this figure, the wall surface becomes hydrophilic by increasing the ρw, which leads to the lower contact angle θw, and vice versa. It should be noted that the linear dependency between the value of Dρ and the computed contact angle is also reported in the literature [20, 46, 47] which the present results show good agreement.Figure 8
Comparison of wetting angle of a droplet on surface obtained by employing SVM-MRT scheme.The second test case for model validation is the simulation of the two-phase fluid flow through a duct which can be considered as a simplified pore of a porous medium. The transmitting of a two-phase flow within a channel causes the wetting phase to flow near the hydrophilic wall which forms a film of wetting fluid on the surface. Consequently, the nonwetting phase is forced to flow through the central core of the channel. The schematic of the two-phase flow inside a 2D channel is demonstrated in Figure9. The parameters a and b in this figure indicate the height of the nonwetting phase and the height of the channel from the centerline, respectively. Herein, the periodic boundary condition is considered in the streamwise direction, and the bounce-back boundary condition is imposed on the solid walls.Figure 9
Geometry and boundary conditions for simulation of two-phase flow through a 2D channel.For this two-phase flow problem, the fully developed flow condition is obtained by applying an external body forceF. Therefore, the analytical solution of the fully developed flow velocity for the two-phase flow in the cross-section of the channel is defined as follows [46]:
(25)ux=F2υwρwb2−y2a<y<bux=F2υwρwb2−y2+F2υnwρnwa2−y20<y<a,where υw, ρw, υnw, and ρnw are the kinematic viscosity and density of the wetting and nonwetting phases. According to Eq. (25), the velocity profile of the fully developed two-phase flow in the channel is related to the dynamic viscosity ratio M=μnw/μw and the saturation of each phase [46], which can be defined as
(26)Sw=b−a/b,Snw=a/b.In the present study, the fully developed velocity profile of the two-phase flow through the channel obtained based on the numerical solution is compared with the analytical data in Figure10 at the flow conditions sw=0.33 and M≈0.1,0.32. The velocity profile of the nonwetting phase predicted by the present MRT-EDM scheme has about 9.5% and 6% deviation from the analytical profile at M=0.1 and M=0.32, respectively. Huang and Lu [48] have reported about 7% difference between their numerical solution and the analytical data for the velocity profile of the two-phase flow through the channel. This comparison shows that the present numerical technique is accurate enough for the simulation of the two-phase flow pattern in 2D pore scales.Figure 10
Distribution of streamwise velocity componentux in a cross-section in the middle of 2D channel with flow condition F=1.5×10−8 at Sw=0.333 and (a) M=0.1 and (b) M=0.32.
(a)(b)To assess the capability of the implemented S-C LBM for computing the two-phase flow characteristics in the pore scale, the relative permeability of the wetting and nonwetting phases through the channel is computed at different flow conditions. Herein, the definition proposed by Yiotis et al. [46] is used for calculating the relative permeability of each phase as
(27)kr,w=12Sw23−Sw,kr,nw=Snw32M+Snw21−32M,where kr,w and kr,nw denote relative permeability of wetting and nonwetting phases, respectively. These relations define that in a 2D channel, kr,w only depends on the saturation of the wetting phase, while kr,nw is related to the saturation of the nonwetting phase and the dynamic viscosity ratio. Figure 11 illustrates that the relative permeability obtained based on the present numerical scheme is in excellent agreement with those computed with the analytical relations at different wetting phase saturations. A similar study is performed at different relative temperatures T/Tc to examine the accuracy of the S-C LBM for predicting the relative permeability magnitude in various density ratios. Figure 12 shows the comparison of the relative permeability computed by the present numerical scheme with those obtained by the analytical solution as a function of the reduced temperature. This study shows a very good agreement between the numerical and analytical solutions with the implementation of the modified interaction forcing scheme with the parameter β=1.18, while using the standard interaction force (β=1) causes a serious deviation between the results at the flow conditions with T/Tc<0.7, which corresponds to the high-density ratios. Therefore, the choice of the interaction forcing scheme becomes important for such flow conditions. This validation study shows that the S-C LBM with employing the MRT-EDM scheme, CS EoS, and the modified forcing scheme with parameter β is efficient and accurate for estimation of the two-phase flow structures, and characteristics through the channel and can be used for simulation of such flow problems in porous media.Figure 11
Comparison of relative permeability of the wetting and nonwetting phases obtained based on present numerical scheme with analytical solution atM=0.1.Figure 12
Comparison of relative permeability of wetting and nonwetting phases at various reduced temperatures for both standard and modified interaction forces.


## 3.3. Two-Phase Flow Transport through a Porous Medium

The present numerical scheme developed based on the S-C LBM with employing the CS EoS, MRT-EDM, and modified forcing techniques is applied for studying the two-phase flow transport through a porous medium. The capillary numberCa and the dynamic viscosity ratio M are considered as two effective parameters that characterize the two-phase flow structure in this study [49]. Herein, the capillary number is defined by
(28)Ca=unwμnwσcosθ,where μnw and unw are the dynamic viscosity and velocity of the advancing nonwetting fluid, respectively, and θ is the fluid-fluid contact angle. Herein, a homogeneous porous medium is generated by a distribution of obstacles in a flow domain with 1000×300 grid points. The pore space area and the height of the throat between the obstacles include 20×20 lattice and 4∼7 grid points, respectively. These considered grid points are required for the spatial discretization of a pore to perform the solution accurately enough for capturing the interfacial dynamics in the pore scale. Hence, the existence of very small pores in the geometry of a porous medium needs a high-resolution grid size and increases the computational cost. The geometry of the porous medium and the boundary conditions of the flow domain are shown in Figure 13. The velocity inlet and pressure outlet boundary conditions are imposed in the left and right boundaries, respectively, and the top and bottom sides are considered periodic. The wetting fluid is defined by the higher density (red color) which invades the nonwetting fluid with lower density (blue color) in the flow domain. Such an invasion produces some fingering regimes that their origin can be the surface or geometrical characteristics of a porous medium in nature, e.g., the nonuniform surface wettability of a sandstone [50]. Herein, the fingering regime during the two-phase flow simulations through the considered porous medium is occurred due to the random generation of the shape of obstacles and the wetting condition of the solid surfaces.Figure 13
Geometry and boundary conditions of porous medium.Figure14 demonstrates the two-phase flow patterns obtained by the present numerical scheme based on the S-C LBM at the dynamic viscosity ratio M=1017 with different capillary numbers. To define the capillary number, the surface tension is calculated by using the Laplace law analysis, which is obtained σ=0.02316. Also, the fluid-fluid contact angle is measured about 175∘ which shows fully hydrophobicity of the invading fluid. This study is performed at the reduced temperature T/Tc=0.482, and the parameter β is set to be 1.241 to preserve the thermodynamic consistency of the solutions. From the results presented in Figure 14, it can be seen that the penetration of the invading fluid at a certain time t=50000 is decreased by reducing the capillary number. These patterns are discussed with more details along with the explanations of Figures 15 and 16 at the end of this section.Figure 14
Two-phase flow patterns in the porous medium obtained by present simulations for flow conditionlogM=3.07 with various capillary numbers at t=50000.Figure 15
Results obtained for displacement pattern of the invading fluid in porous medium by the present S-C LBM for flow conditionlogM=3.07 with various capillary numbers in comparison with the schematic of phase diagram presented by Zhang et al. [51] and Lenormand et al. [49].Figure 16
Instantaneous displacement pattern of two-phase flow through the porous medium obtained by present simulations for flow conditionlogM=3.07 and logCa=−1.85 at three different times.To compare the obtained results for the two-phase flow characteristics through the porous medium in a similar situation, a normalized timeT∗ is defined by division of the time step to the total simulation time of each case. The total simulation time is counted until the first finger of the invading fluid touches the outer side of the porous medium for each case, which is called breakthrough [49]. Figure 17 illustrates the variation of the pressure difference Δp between the inlet and outlet of the porous medium versus the normalized time T∗ at different capillary numbers. As seen in this figure, by increment of the capillary number in a constant dynamic viscosity ratio, the invading fluid has enough pressure difference to overcome the throat threshold pressure pc presented by Lenormand et al. [49], which is known as a resistance factor. The throat threshold pressure is defined by the following relation:
(29)pc=2σcosθr,where σ is the surface tension, θ is the fluid-fluid contact angle, and r is the throat diameter. This parameter also affects the nonwetting phase saturation. In the present study, by considering the minimum and maximum throat radius, the threshold pressure is altered between ‐0.006591 and ‐0.011535 for r=7 and r=4, respectively. The negative sign indicates the resistance pressure against the invading fluid.Figure 17
Predicted pressure differenceΔp between the two sides of the porous medium in streamwise direction versus normalized time T∗ for flow condition logM=3.07 with various capillary numbers.Figure18 shows the nonwetting phase saturation as a function of the normalized time T∗. The results obtained demonstrate that the saturation of the nonwetting fluid is enhanced by increasing the capillary number, which is consistent with the previous works [49–51]. The final values of the nonwetting phase saturation and the pressure difference at the breakthrough condition are illustrated in Figure 19 at different capillary numbers. It is obvious that the variation of the nonwetting phase saturation is almost linear which is in agreement with the results obtained by Zhang et al. [51] and Liu et al. [50]. The agreement of the present results in comparison with those reported in the literature confirms that the implemented S-C LBM with the MRT-EDM scheme besides the modified interaction forcing method and the CS EoS is capable for simulation of two-phase flows at high-density and -viscosity ratio through porous media at different flow conditions.Figure 18
Nonwetting phase saturation versus normalized timeT∗ for flow condition logM=3.07 with various capillary numbers.Figure 19
Comparison of nonwetting fluid saturation and pressure differenceΔp at the breakthrough time for flow condition logM=3.07 with various capillary numbers.The structure of the two-phase flow through the porous medium is also investigated at different flow conditions. The displacement pattern of the invading fluid is categorized as the capillary fingering, viscous fingering, and stable displacement, depending on the capillary number and the viscosity ratio between two phases [52]. A phase diagram is used to describe these flow structures by logCa and logM for a gas-water flow system (see Figure 15). In the present simulations based on the S-C LBM, the magnitude of the viscosity ratio is set as logM>0. Therefore, the capillary fingering and stable displacement patterns are expected in the porous medium considered. As can be seen in Figure 14, the two-phase flow pattern is altered from the stable displacement to the capillary fingering by a decrement of the capillary number. Figure 16 demonstrates the present results for instantaneous two-phase flow displacement at different times for logCa=−1.85, in which the pattern behavior is almost stable displacement. In this flow condition, the interface of the invading fluid penetrates into the porous medium along a straight line, despite some irregularities due to the pore-scale geometry. This behavior can be the result of the dominant effect of the viscous force and the higher pressure difference that cause to fill the porous medium uniformly and simultaneously. By decreasing the capillary number, it is observed that the invading fluid is extended in the shape of irregular fingers in the flow field due to the dominant effects of the capillary force in comparison with the viscous force. Such a capillary fingering of two-phase flow through the porous medium is shown in Figure 20 at logCa=−3.5. In this case, due to the uniform pressure distribution throughout the domain, the nonwetting fluid only penetrates the pore-scale area with a lower throat pressure threshold. As a result, the direction of the flow displacement in the porous medium is irregular. In Figure 20, it can be seen that the capillary fingers of the nonwetting fluid form several loops that trap some slugs of the wetting fluid among the porous medium. These two-phase flow structures obtained using the present numerical scheme are in consistent with those reported in the literature [49, 50].Figure 20
Instantaneous displacement pattern of two-phase flow through the porous medium obtained by present simulations for flow conditionlogM=3.07 and logCa=−3.5 at various times.The interfacial area between two phases is quantified for the two-phase flow through the porous medium at different capillary numbers. This parameter is an important factor that influences the mass and energy transfer among the phases [50, 51]. In the present 2D studies, the interfacial length is measured instead of the interfacial area. Herein, due to the extreme hydrophilicity of the wetting fluid, the length of the wetting films around the obstacles is also added to the length of the measured interface between the phases for computing the interfacial length. To have comparable results for this parameter, a specific interfacial length is calculated as the ratio of the measured interfacial length to the pore area of the porous medium. Figure 21 shows the variation of the specific interfacial length as a function of the nonwetting phase saturation at different capillary numbers. The present results demonstrate that there is a linear relation between two parameters in all the capillary numbers considered which is in good agreement with the previous experimental finding of Zhang et al. [51].Figure 21
Comparison of specific interfacial length as a function of nonwetting phase saturation at different capillary numbers for flow conditionlogM=3.07.


## 4. Conclusion

In this work, several different combinations of collision operators, EoSs, and forcing schemes of the pseudopotential multiphase LBM are employed to investigate their performance and capability for computing some two-phase flow characteristics at high-density and –viscosity ratios with the aim of applying for porous media applications. Herein, the SRT and MRT collision models, Shan-Chen (SC) and Carnahan-Starling (CS) EoSs, and SVM and EDM forcing schemes are considered to study their accuracy and numerical stability. An appropriate modification in the forcing schemes is also used to remove the thermodynamic inconsistency in the simulation of two-phase flow problems studied at low reduced temperatures. The simulation results obtained for the equilibrium state of a droplet, the two-phase flow characteristics in a channel, and through a porous medium are presented and discussed. The conclusions and remarks regarding the present study are as follows:(1)
The Laplace law is satisfied for simulation of the stationary droplet by using the S-C LBM coupled with the SRT-SVM, SRT-EDM, and MRT-EDM schemes which demonstrates the accuracy of the numerical schemes employed for predicting the characteristics of a stationary two-phase flow system(2)
The results obtained for the stationary droplet at low-density ratios (high reduced temperatures) by the SRT-SVM, SRT-EDM, and MRT-EDM schemes with employing both the SC and CS EoSs are in good agreement with the analytical solution on the coexisting curve. However, at high-density ratios (low reduced temperatures), the results have a significant deviation from the analytical solutions with the implementation of the forcing schemes in the standard form. By employing theβ-modification in the forcing schemes, the density values obtained by the present numerical simulations are in good agreement with those obtained by the analytical solution in a wide range of reduced temperatures(3)
The results presented in this study show that with applying the CS EoS for simulation of considered two-phase flow problems at high-density ratios, the spurious velocity is significantly decreased in the flowfield. This result demonstrates that the CS EoS is more realistic and efficient than SC EoS for simulation of two-phase flows at high-density ratios due to producing less spurious currents at the interfacial region(4)
It is shown that the MRT-EDM with the improved forcing scheme produces minimum spurious currents among the other methods. By considering the results presented in this study, the advantages of employing the MRT scheme withβ− modification of forcing schemes are obvious for simulation of a two-phase flow system at low relaxation times (τ⟶0.5)(5)
According to the previous remarks, the present study demonstrates that the S-C LBM with employing the MRT-EDM scheme, CS EoS, and the modified forcing scheme is efficient and accurate for estimation of the two-phase flow characteristics. Therefore, this combination of the implemented numerical techniques is recommended and applied for studying two-phase flows through the channel and porous medium(6)
The validation study of two-phase flow through the channel shows that the S-C LBM with employing the MRT-EDM scheme, CS EoS, and the modified forcing scheme with parameterβ is efficient and accurate for estimation of the two-phase flow structures and characteristics in pore scales(7)
The agreement of the present results obtained for the simulation of two-phase flow through the porous medium at different conditions in comparison with those reported in the literature confirms that the suggested S-C LBM is capable for predicting the displacement patterns, the nonwetting phase saturation, and the penetration structure of the invading fluid in porous media applications. The extension of the present method in the three-dimensional framework is straightforward and will be presented in future works



---

*Source: 5668743-2021-08-06.xml*