Next Article in Journal
Data-Driven Modeling of Geometry-Adaptive Steady Heat Convection Based on Convolutional Neural Networks
Next Article in Special Issue
Dynamics of a Laser-Induced Bubble above the Flat Top of a Solid Cylinder—Mushroom-Shaped Bubbles and the Fast Jet
Previous Article in Journal
Evaluation of the Thermofluidic Performance of Climatic Chambers: Numerical and Experimental Studies
Previous Article in Special Issue
LES of Particle-Laden Flow in Sharp Pipe Bends with Data-Driven Predictions of Agglomerate Breakage by Wall Impacts
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Equation of State’s Crossover Enhancement of Pseudopotential Lattice Boltzmann Modeling of CO2 Flow in Homogeneous Porous Media

by
Assetbek Ashirbekov
1,
Bagdagul Kabdenova
1,
Ernesto Monaco
2 and
Luis R. Rojas-Solórzano
1,*
1
Department of Mechanical and Aerospace Engineering, School of Engineering and Digital Sciences, Nazarbayev University, Nur-Sultan 010000, Kazakhstan
2
Engineering Software Steyr (ESS), 4400 Steyr, Austria
*
Author to whom correspondence should be addressed.
Fluids 2021, 6(12), 434; https://doi.org/10.3390/fluids6120434
Submission received: 23 June 2021 / Revised: 17 November 2021 / Accepted: 19 November 2021 / Published: 1 December 2021
(This article belongs to the Collection Advances in Flow of Multiphase Fluids and Granular Materials)

Abstract

:
The original Shan-Chen’s pseudopotential Lattice Boltzmann Model (LBM) has continuously evolved during the past two decades. However, despite its capability to simulate multiphase flows, the model still faces challenges when applied to multicomponent-multiphase flows in complex geometries with a moderately high-density ratio. Furthermore, classical cubic equations of state usually incorporated into the model cannot accurately predict fluid thermodynamics in the near-critical region. This paper addresses these issues by incorporating a crossover Peng–Robinson equation of state into LBM and further improving the model to consider the density and the critical temperature differences between the CO2 and water during the injection of the CO2 in a water-saturated 2D homogeneous porous medium. The numerical model is first validated by analyzing the supercritical CO2 penetration into a single narrow channel initially filled with H2O, depicting the fundamental role of the driving pressure gradient to overcome the capillary resistance in near one and higher density ratios. Significant differences are observed by extending the model to the injection of CO2 into a 2D homogeneous porous medium when using a flat versus a curved inlet velocity profile.

1. Introduction

Carbon Dioxide (CO2) and other greenhouse gases (i.e., methane, nitrous oxide, among others) are key contributors to climate change. According to the Energy Information Administration (EIA) estimates, the global emission of CO2 alone can rise to 6.41 billion tonnes by 2030 [1]. Therefore, various mitigation measures have been proposed to tackle global climate change. Among those proposed solutions, we use CO2 in different applications such as Hot-Dry Rock (HDR) systems, heat engines, selective extraction processes, reclamation processes of metals from industrial effluents, and chemical reactions [2,3]. Nevertheless, CO2 sequestration is recognized as one of the key strategies for reducing CO2 emission rates. Due to its low critical temperature (31 °C), CO2 is in supercritical (for example, in CO2 sequestration) or in near-critical (as in air conditioning and cooling systems [4,5]) state.
The wide usage of CO2 around the critical state requires a deeper investigation of the fluid thermodynamics under and above critical conditions, affecting fluid transport properties [1]. In addition, it facilitates understanding the industrial and laboratory processes wherever the fluid is involved since the behavior of fluids in supercritical or near-critical conditions dramatically differs from those with constant properties.
The fluid flow under different conditions can be studied experimentally and numerically. However, as the experimental study may not always be feasible and is time-consuming, the numerical analysis can be used as a starting point. It can even replace the experiments to provide valuable insights into the complex flow patterns.
Computational Fluid Dynamics (CFD) techniques such as the finite volume, volume of fluid, and level-set methods solve the discrete Navier–Stokes equations. However, these methods become computationally expensive when applied to the interface tracking of transient multiphase flows [6,7]. The mesoscale Lattice Boltzmann Model (LBM), specifically the LBM employed in the current work, does not explicitly track the phase interface; instead, it includes inter-particle forces naturally preserving the interfaces [8]. Additionally, it accommodates complex boundaries such as porous media and narrow channels. Furthermore, the pseudopotential LBM can incorporate real Equations of State (EoS) such as Van der Waals (VdW), Carnahan–Starling (CS), and Peng–Robinson (PR) EoS [9], which, in turn, allows the model to match the fluid characteristics more closely.
Despite its advantages, the pseudopotential LBM still faces issues when high density and viscosity ratios are included in multicomponent-multiphase flows in complex geometries. For example, studies on displacement patterns in porous media do not consider density contrast and differences in the fluids’ critical properties [10,11,12,13,14]. Hence, most numerical studies do not consider that CO2 and H2O are in different thermodynamic states. For example, CO2 is usually in a supercritical state during sequestration, while H2O is in a liquid state, limiting the accuracy of the numerical results.
In addition, most studies apply classical EoS to analyze the thermodynamic properties of the fluid in different ranges of temperature [9,15,16,17]. Though classical cubic EoS predict well the fluid properties away from the critical region, they fail to accurately capture the fluid behavior in the near-critical region or at a supercritical state because of long-scale fluctuations in density that cover a wide range of temperatures [1,18,19,20,21,22]. In this region, analytical solutions provided by classical EoS are no longer accurate. As a remedy for this pitfall, a so-called “crossover” formulation of the EoS was proposed, incorporating the non-analytical scaling laws and predicting fluid properties in the critical region more accurately [1,18,21]. Any classical EoS can be re-formulated in a crossover fashion. In this work, the crossover formulation is incorporated into PR EoS. As a result, PR EoS better characterizes the specific fluid through an acentric factor. This approach’s advantage is a hybrid equation that accounts for both non-analytical and analytical fluid behavior, close and far from the critical point, respectively. It can also transform back to ideal gas EoS at a zero-density limit.
Our previous study [23] incorporated the crossover PR EoS into the standard pseudopotential multicomponent LBM. We showed that crossover PR EoS improved the model accuracy for fluids at near/supercritical conditions compared to the classical PR EoS. Another previous work [24] showed that a higher water wettability promotes a higher effective CO2 penetration through a porous medium. We then applied the improved model formulation to study the transport of a supercritical CO2 bubble within a simplified hydrophilic porous medium saturated with H2O, proving the capacity of the newer model crossover formulation.
This work extends our work by introducing the PR EoS multicomponent LBM to study the penetration mechanisms of a continuous front of supercritical CO2 through a pore-scale homogeneous domain saturated with H2O. Firstly, a 2D narrow channel is considered to model the penetration of a supercritical CO2 bubble through the channel with a significant density ratio. Secondly, the sequestration of a continuous flow of CO2 is modeled by injecting it through a water-saturated 2D homogeneous hydrophilic porous medium (contact angle of 70°), varying the pressure gradient around the pressure drop needed to overcome the capillary resistance. Finally, we study the CO2 injection under a fixed flow regime and a mildly hydrophobic 2D porous medium (contact angle 85°) to demonstrate the effect of the injection profile in the CO2 penetration into the water-saturated porous medium. The porous media was constructed by placing circle-obstacles homogeneously. Hence, this study represents the first investigation, to the very best of our knowledge, featuring a crossover PR EoS formulation of LBM to shed some light on the CO2-H2O-surface interaction and CO2 sequestration in a porous domain at near supercritical conditions.
This paper is organized as follows: Section 2 introduces the multicomponent LBM and crossover PR EoS formulation. Section 3 presents the model’s performance for the static H2O droplet immersed in CO2 and dynamic conditions of CO2 injection at supercritical conditions into H2O-saturated capillary domains. The final discussion and conclusions are presented in the last section of the paper.

2. Methodology

2.1. General Lattice Boltzmann Model

The LBM equation represents the evolution in time of distribution functions because of fluid particles streaming and colliding.
The most general formulation of LBM for a multicomponent system is presented in Equation (1):
f i j x + e i t , t + t f i j x ,   t = 1 τ j f i j x ,   t f i j , e q x , t
where f i j x ,   t is the density distribution function at point x and time t with the discrete velocity in i direction for component j. In Equation (1), the left-hand side represents the streaming step, while the right-hand side is the collision step, which can assume different forms: in this work, we opted for the simplest one, the Bhatnagar–Gross–Krook (BGK) collision operator [25], which each fluid species is characterized by a single relaxation time τ j which is related to the whole mixture kinematic viscosity through the relation (Equation (2)) [26]:
v = c s 2 j χ j τ j 0.5 ,
where c s is the lattice sound speed and χ is the local value of the mass fraction. The collision operator relaxes the local distribution function towards its equilibrium value, f i j , e q x , t which depends on macroscopic quantities, as prescribed by Kinetic Theory, and is given by Equation (3):
f i j , e q x , t = ω i ρ j 1 + e i u j c s 2 + e i u j 2 2 c s 4 u j 2 2 c s 2 .
In this work, we used a two-dimensional formulation with nine local speeds at each computational node. This set of speeds, referred to as a lattice, must be carefully chosen to recover the Navier–Stokes equation. If the lattice vectors possess different modules, weighting factors w i need to be defined to achieve this goal. The number of dimensions and speeds defines a generic lattice. We chose the standard D2Q9 one, depicted in Figure 1, characterized by weighting factors 4 9 i = 0 ,   1 9 i = 1 , 2 , 3 , 4 ,   1 36 i = 5 , 6 , 7 , 8 .
Macroscopic quantities, including density ρ j , can also be obtained as discrete moments of the particle distribution f i :
ρ j = i f i j ,
ρ j u j = i f i j e i ,
where u is the macroscopic velocity vector. The difference between continuous distribution function f and f i is the discrete nature of the last one; subscript i refers to the discrete set of velocities e i , and superscript j refers to a component.

2.2. Pseudopotential Lattice Boltzmann Model

In Shan-Chen’s or pseudopotential model, the interaction between the particles is modeled using pairwise interaction potentials. These interaction potentials are functions of density for each component. If we define the interaction matrix  G j k between two generic species j and k at a generic node x , ψ as a pseudopotential (“effective mass”), then the force acting on jth component can be written as [27]:
F i n t j x , t = ψ j x , t k G j k x i ω i ψ j x + e i t , t e i
The interaction force is given by Equation (6), and it is applied between the nearest neighbor particles and leads to the pressure P j in the EoS:
P j = ρ j c s 2 + k g j k 2 c s 2 ψ k 2 ,
where the first term represents the ideal EoS and the second term accounts for repulsive (positive g j k ) and attractive forces (negative g j k ) between the fluid particles. For a two-component system, usually the diagonal terms of interaction matrix, accounting interactions between particles of same species are equal to zero, while interparticle terms are chosen positive to enforce repulsion, as shown in Equation (8).
G x j k = g j k = g k j ,   j k   0 ,   j k .
The pseudopotential ψ in Equation (7) can be tailored to obtain a given equation of state. Pseudopotential can also be obtained by rearranging Equation (7):
ψ j ρ j = 2 P j ρ j c s 2 g c s 2
Yuan and Schaefer incorporated different EoS such as VdW, PR, and CS into the pseudopotential model through Equation (9). In this case, the term g is used to ensure the positive sign under the square root [9].
The fluid spreading on a solid surface depends on the strength of fluid–solid interaction. The surface wettability is modeled using Equation (10) proposed by Martys and Chen [27] in multicomponent form:
F w e t j x = g w a l l j ρ j x i ω i s x + e i e i
The g w a l l represents the fluid–wall interaction strength. For example, when g w a l l is positive, the wettability decreases, while if the g w a l l is negative, the wettability increases. s x + e i is the “switch” function equal to 1 in the presence of the solid surface at a lattice node x ; otherwise, it is 0. The following method is based on density instead of potential ψ x [28].
The total force acting on a fluid particle is described with the equation:
F t o t j = F w e t j + F i n t j + F b
The F b is the external body force, which can have a constant value if applied to drive the flow, for example, or zero, if no external body force is present. The present work adopts the velocity-shift force scheme, originally used by Shan and Chen [8]. In this scheme, forces are incorporated by updating the velocity u in the equilibrium distribution function, replacing it with u e q , giving:
u e q j = u j + τ j ρ j F t o t j Δ t ,
where the actual physical velocity u is taken from Equation (5). Other available schemes are the Exact Difference Model (EDM) [17,29] and Guo’s [30]. These schemes result in a further term at the right-hand side of (1), the source term S i x , t .

2.3. Crossover Peng–Robinson EoS

According to Landau’s theory, free energy is an analytic function of the order parameter. However, phase transition disregards the order of the parameters, and the partition function becomes the sum of all possible values. Consequently, an analytical solution in the region near the critical point cannot be formulated [18,31]. On the other hand, the non-dimensional expression of Helmholtz’s free energy is expressed in Equation (13) in the crossover formulation of EoS. Moreover, it is made up of two terms [18], where the first term represents the non-analytical part A ¯ , and the second term considers the analytical part μ ¯ 0 T that depends on temperature.
A ¯ T , V = A T , V R T = A ¯ T , V V V c P ¯ 0 T + μ ¯ 0 T ,
where P ¯ 0 T = P T , V c V c RT is the non-dimensional pressure at the critical isochore ( V = V c ), and both T = T r 1 and V = V / V c 1 are dimensionless deviations from classical critical parameters. The critical parameters are found by solving Equation (14):
P V T c = 0 ,   2 P V 2 T c = 0 ,   P c V c R T c = Z c ,
where Z c is the critical compressibility factor.
After imposing the condition expressed in Equation (15), the critical part of the Helmholtz free energy takes the form shown in Equation (16) with the introduction of dimensionless coefficients b and critical exponent α :
A ¯ T , V = 0 = 0 ,   and   A ¯ V T V = 0 = 0 ,
A ¯ T , V = ln V b 1 + 1 T c T 2.077873 α T ln V / b 2 + 1 V / b 3 + 1 + 1 + V b 1 T c T 0.45724 α T V Z c b 2   b 3
The analytical part is given as:
μ ¯ 0 T = P ¯ 0 T ln b 1 + T c T 2.077873 α T ln b 2 b 3 + A ¯ 0 T
The term A ¯ 0 T in Equation (17) is expressed as A ¯ 0 T = A 0 T / R T and defined as a non-dimensional function of temperature that plays a role in reproducing ideal gas behavior in the limit of zero density. Helmholtz free energy in crossover form uses the renormalized form of classical reduced parameters T and V as shown in Equations (18) and (19):
τ ¯ = τ Y α 2 1 + 1 + τ τ c r i t Y 2 2 α 3 1 ,
η = η Y γ 2 β 4 1 + 1 + η η c r i t Y 2 α 2 1 ,
where τ = T / T c r i t 1 is the temperature difference from the real critical temperature ( T c r i t ) and η = V / V c r i t 1 is the deviation of the volume from the real critical volume ( V c r i t ). Both Equations (19) and (20) use the real critical temperature and volume, and their second term accounts for the difference between classical and real critical parameters. Moreover, these shift parameters are expressed as τ c r i t = T c r i t / T c 1 and η c r i t = V c r i t / V c 1 . As in [22], we assume that T c = T c r i t and therefore, τ c r i t = 0 . Another important parameter in Equations (17) and (18) is a crossover function, Y , that accounts for volumetric changes in the EoS depending on how far the system is from the critical region. If Y q 0 , when the system is at the critical point, the crossover function should approach zero, Y q 0 , and Y q 1 when the system is away from the critical point. The crossover function proposed by Feyzi et al. [21] satisfies these requirements, and we use the same formulation in Equation (20).
Y q = q q + 1 2 1 ,   q = q   exp 12.2 12.2 T T c r i t ,
where q is the renormalized parameter that measures the distance between the system volume and the critical one; q is found solving the Equation (21), as in [20].
q 2 τ G i 1 p 2 4 b 2 1 τ q 2 G i = b 2 η 1 + v 1 exp 10 η + d 1 τ m 0 G i β 2 Y 1 2 β 1 ,
b 2 = p 2 = 1.359 and G i is the Ginzburg number, while d 1 , v 1 , and m 0 are system-dependent parameters. Equation (21) has been solved numerically using MATLAB, and as it has several solutions, we take the highest positive value among all possible solutions. After finding q and Y q , the renormalized parameters τ ¯ and η ¯ are substituted into Equation (16) to obtain the critical form of the Helmholtz free energy.
Finally, the crossover PR EoS is obtained after substituting all the parameters and taking the differentiation with respect to the volume, Equation (22):
P = A V T = R T V c V c V c r i t A ¯ η T + P ¯ c T + V c V c r i t
As mentioned earlier, the critical shift τ c r i t is set to 0 and η c r i t , G i , d 1 , v 1 , m 0 , a 20 , and a 21 are parameters that depend on the type of fluid, and their values are shown in Table 1 for CO2 and H2O. The crossover PR EoS is incorporated into the LBM through Equation (22). For CO2, all parameters are taken from [21], while all other fluid parameters are taken from [20]. Furthermore, the physical critical parameters are derived from the NIST database [32].

3. Results

Our previous study [23] conducted several 2D-droplet tests at different reduced temperatures, and their respective coexistence curves have been constructed. The detailed comparison of the thermodynamic behavior of CO2 and H2O at various temperatures predicted by PR and crossover PR EoS showed that PR EoS dramatically improves the model accuracy, which is more noticeable in the liquid phase. PR generally overestimates the density of the curve’s liquid branch, although it predicts the gas phase accurately. As temperature reduces, the difference in liquid density predicted by the PR EoS increases, following the experimental data’s behavior [23]. Similar results are also noted in other studies [33,34]. It is also observed that the crossover PR EoS matches experimental data much better than the classical PR EoS at different temperatures close to the critical point. The results show that crossover PR EoS is more accurate than classical PR EoS because it considers the non-linear behavior of fluid species in the critical regions [23].

3.1. Multicomponent System

Modeling the multicomponent system with actual fluid parameters is challenging mainly due to a lack of thermodynamic consistency and potential large spurious currents at interfaces. Moreover, the magnitude of these spurious currents increases with the density ratio. Therefore, studies that involve complex systems such as flow in porous media are mostly limited to fluids with similar densities [10,35,36].
In our previous study [23] and the current work, we implement the improvements to the MCMP model following Stiles and Xu’s [15] approach, which improves the model in terms of numerical stability and allows a density ratio higher than unity. The improved model is then applied to simulate the stationary 2D droplet immersed in supercritical CO2, accounting for density contrast and the fluid’s critical temperature difference. We use a so-called scalar multiplier as suggested by Stiles and Xue [15] with the assumption that t = x = 1 . The pressure term in the numerator of Equation (9) is replaced with P i * = k P i . We ensure stability of the model by using relaxation time of 1, which is a requirement of the formulation.
In addition, using crossover PR EoS helps better match the binodal curve. It considers fluid properties by applying the acentric factor and system-dependent parameters to improve the critical region’s curve adjustment. The acentric factor for H2O is 0.344, and the set of parameters used for CO2 is given in Table 1.
Stiles and Xu investigated the effect of scaling parameters on the interaction force model and the coexistence curve (P versus V). As expected, the curve’s amplitude decreases when k < 1 and the volume of the liquid phase remains constant. This parameter only affects the volume of the gas phase. This effect becomes more intense at lower temperatures, thus limiting the range of application to moderate and high temperatures [15,23]; in this study, we assume that k = 0.75.
To account for differences in real fluids, critical temperatures T c ,   H 2 O T c ,   C O 2 2 , we use the set of parameters presented in Table 2. Previous validation [23] shows supercritical CO2 co-existing with subcritical H2O at 0.071 (lattice units—LU—temperature), corresponding to T r of 2.07 and 0.97, respectively. Hence, we first set the critical lattice temperature of H2O to 0.07292 and then use the relation in Equation (20) to find the critical temperature of CO2.
T C O 2 ,   L U = T H 2 O ,   L U T C O 2 ,   p h y s T H 2 O ,   p h y s
Figure 2 presents a 2D static water droplet surrounded by CO2 in the center of a rectangular domain with the size of 201 × 201 LU2, where the system’s temperature is 0.071 in LU. Here, g i j = g j i = 0.001 , as it gives good stability and is recommended in [15], and the drop radius is 30 LU.
Figure 2 shows a slight trace of H2O in the CO2 and vice versa, but the CO2, being the primary phase, mainly dictates the interfacial dynamics. The presence of trace elements in small amounts helps to reduce the initial gradient at the phase interface and improves numerical stability.

3.2. Contact Angle Test

This section presents the contact angle tests and shows that the improved model is valid for simulating fluid flows involving fluid–solid interactions.
Figure 3 shows that the multicomponent multiphase model can reproduce a wide range of contact angles. The properties of both fluids are the same as described in Section 3.1. The interaction strength parameter between the solid and fluid g CO 2 , wall is 0.2, and for the g H 2 O , wall is −0.2, accordingly.
In this case, the H2O droplet behaves as a wettable fluid, and the CO2 is non-wettable. The values of the g CO 2 , wall and g H 2 O , wall used in this study to perform the contact angle tests are presented in Table 3. The results show that the pseudopotential model can model a wide range of contact angles by adjusting the solid–fluid interaction parameters, g CO 2 , wall and g H 2 O , wall .

3.3. Penetration Process in 2D Narrow Channel

In this section, we analyze the performance of the multicomponent LBM model using the crossover PR EoS formulation in the visualization of the penetration of CO2 into an H2O-saturated 2D funneled channel that helps approximate the contractions in a pore network. CO2 is a non-wetting phase, and H2O is set as a wetting phase. The geometry of the domain is depicted in Figure 4a–c filled with immiscible H2O and CO2 fluids subject to a streamwise pressure gradient (body force) and inlet/outlet periodic boundary condition. The domain size is 41 × 151 LU2, and the channel is narrowed by placing two-block obstacles with a height of 40 LU and a width of 10 LU. Figure 4a also presents the initial condition with two large CO2 bubbles immersed in H2O and placed initially at the entrance and exit of the channel. H2O and CO2 fractions are 0.61 and 0.39, respectively, and the system’s temperature is 0.071 in LU.
The density ratio H2O/CO2 is 4, and H2O is set as the wetting phase with a contact angle of 70°. The top and bottom sides of the domain are set as periodic and left, and the right walls are set with a no-slip boundary condition. The streamwise body force is then applied throughout the domain with magnitude G f = 4.8 × 10−5 just enough to overcome the capillary resistance. Figure 4b–d show the evolution of the flow with the bubble squeezing in a periodic pseudo-steady regime. The results show that the velocity of the CO2 bubble is almost twice faster in a narrower section of the channel, as observed in Figure 5. From further analysis, it is observed that bubbly flow tends to form when CO2 saturation is lower or equal to 0.4. At CO2 saturation between 0.4 and 0.7, the annular flow is formed where wetting phase water stays closer to the channel wall. These findings are similar to what was observed in another similar study where free energy LBM was used instead [37,38].
The model is further challenged by varying the pressure gradient from insufficient to a much larger value than needed to overcome the capillary resistance. Figure 6 presents the snapshots after a long transient simulation (297,500 ts) to reach a pseudo-steady solution. The boundary conditions and wettability are like those in the previous case. Figure 6a shows that under a pressure gradient of 3.5 × 10−5, the CO2 cannot enter into the narrow channel due to insufficient driving force to overcome the capillary resistance. The body force is further increased to 4 × 10−5 (Figure 6b), slightly above the needed capillary pressure. The CO2 squeezes through the narrow channel until an annular flow is developed with a stable symmetrical shape. The centerline velocity in the middle (narrow) domain section is 0.016 LU/ts, twice higher than in the broader section of the domain. When the body force increases to 6 × 10−5, the multicomponent flow becomes unstable, and the inner CO2 starts oscillating in the spanwise direction, depicting a wavy shape surrounded by a film of H2O (Figure 6c).

3.4. Penetration Process in the 2D Pore Network

The modeling CO2 penetration into a 2D porous media under different body forces was also studied using the crossover PR EoS formulation embedded in the multicomponent LBM. The porous media is constructed by placing several previously saturated cylinders with H2O. The schematic diagram of the porous region is shown in Figure 7a, which has a porosity of 0.78 and a pore throat of 10 LU. The domain size is 191 × 81 LU2, with the horizontal length of the pore network being 90 LU. Top and bottom walls are set as a no-slip boundary condition, while both left and right sides are set as periodic. The computational domain is subject to a streamwise pressure gradient to drive the CO2 into the porous medium saturated with H2O (wettable fluid) with a contact angle of 70°. The density ratio is set to 1, and the system’s temperature is 0.071 in LU. Figure 7b presents the case with body force 2.5 × 10−5 at 113,400 ts. In this case, the driving force is not enough to overcome the capillary resistance.
We further increase the driving force to 4.5 × 10−5, and it is observed from Figure 8 that CO2 can then easily intrude the porous media with a flow pattern that resembles a stable displacement process.
Next, we compare the effect of the front profile of the intruding fluid (CO2). According to Fakhari et al. [39], the distance from the inlet to the porous section of the micromodel has a significant effect on numerical simulations. The parabolic velocity profile better matches the initial conditions between the numerical simulations and experiments. To achieve a parabolic profile, we extended the length of the inlet region to 85 LU. To observe the effect of the front shape of intruding fluid on displacement patterns, we also placed the inlet region closer to the pore network at 7 LU from the first row. The initial conditions for both cases are shown in Figure 9a and Figure 10a. In both cases, the overall domain size is 401 × 201 LU2, with the pore network length of 229 LU in the streamwise direction. The inlet velocity is 0.02 LU/ts as proposed by Zou and He [40], applied on the left bound of the domain (inlet), and the outflow boundary condition is applied on the right end. Both top and bottom walls are set as a no-slip boundary condition, and the contact angle is set to 85°. The density ratio is again set to 1, and the system’s temperature is 0.071 in LU.
Figure 9a shows the initial condition for the case where the front of the intruding fluid (CO2) is initially placed 75 LU away from the pore network (i.e., 10 LU inside the domain). The parabolic profile is formed at 3300 ts (Figure 9b). Tiny H2O droplets are observed at 6100 ts (Figure 9c). Droplets evaporate after 2000 ts, and Figure 9d shows that the flow pattern resembles a fingering displacement at 9000 ts. The process occurs at near-constant pressure conditions around the evaporation region. Thus, the model is expected to be capable of capturing the phenomenon despite the inherent isothermal assumption. Figure 10a presents the case with the flat front of the intruding fluid placed 7 LU away from the porous section. Compared with the previous case, the front profile of the velocity remains flat. In contrast, with the parabolical front profile, the velocity is higher in the middle section of the porous media. These results suggest that the distance from the inlet to the porous section can affect the displacement pattern.
A grid-independence analysis is developed by introducing two more refined models of 601 × 301 LU2 and 801 × 401 LU2, with geometry, inlet speed, viscosity, and timestep span scaled accordingly. The average velocity magnitude along a line drawn over the last row of obstacles, as shown in Figure 11, was utilized as a comparison parameter for the three grid sizes. Results of the velocity probe are shown in Figure 12, taken from lattice time of 0 to 7100.
As shown in Figure 12, velocity magnitude is very similar in all grid sizes, deviating only slightly in later time-steps, demonstrating the suitability of any of the three grids in our analysis.
In another comparison between the grids, we used the CO2 flux through the probe line for an extended period to permit the CO2 to reach the probe line. This parameter was calculated as the fraction of CO2 multiplied by the time-average velocity magnitude integrated over the whole probe line. Results, shown in Table 4, highlight the small error between contiguous grid resolutions, justifying our initially selected grid for the analysis.

4. Discussion and Concluding Remarks

In this study, the crossover PR EoS was incorporated into the pseudopotential LBM, and an improved model was further applied to study the CO2 penetration process in capillary water-saturated domains. The improved model’s performance is analyzed through several tests. Firstly, a static weightless water droplet is suspended in equilibrium and immersed in CO2 to prove the stability of the numerical model. Secondly, a sessile water droplet in equilibrium with the surrounding CO2 is presented to demonstrate the stability and accuracy of the model in determining the static contact angle. Both static simulations showed the model’s ability to describe different wettability regimes with sufficient numerical stability.
Then, several cases are presented to analyze the injection of supercritical CO2 in a water-saturated capillary domain. In the first place, a validation of the model around critical conditions is assessed through a single pore or capillary channel, demonstrating the accurate prediction of the threshold pressure gradient needed to overcome constraints to allow a CO2 bubble to squeeze through the narrow channel. Hence, the injection process was analyzed at different pressure gradients and contact angles. After the necessary pressure gradient was met, the non-wettable CO2 traveled through the center of the channel surrounded by the wettable fluid (H2O), depicting a centerline velocity at the narrow section almost twice higher than through the main channel. These results agree with previous work [37]. Increasing the pressure gradient much farther beyond the minimum capillary pressure affected the hydrodynamic stability of the CO2 stream, demonstrating the existence of a secondary stability threshold of pressure gradient above which a transient spanwise oscillation sets in.
Next, the CO2 sequestration process was further analyzed around supercritical conditions where the non-wetting fluid is injected into a complex, but homogeneous 2D pore network saturated with wetting H2O. Thus, we investigated the effect of a supercritical CO2 front profile on its displacing through a 2D water-saturated pore network. The results show that the intruding fluid’s front profile plays an important role and influences the overall CO2-H2O flow pattern through the pore networks. That observation was also highlighted in [41], where it was observed that the parabolical profile is more likely to happen in actual CO2 sequestration, and this velocity profile returns results that are closer to experimental data than obtained for the flat inlet profile. Therefore, our results proved that with enough development length allocated in the inlet region of the computational domain, it is reasonable to expect a more realistic parabolical profile at the entrance of porous networks during the CO2 sequestration numerical modeling.
These preliminary results show that the improved model is numerically stable when applied to reproduce moderately complex geometries and physical flow patterns, as expected in realistic conditions. However, our study is limited to qualitative results and a viscosity ratio of one and isothermal conditions. In future studies, we plan to continue exploring the numerical stability of the new formulation for more complex 3D porous media and even more realistic fluid–fluid conditions.
For a detailed expansion of the crossover PR EoS, please refer to Appendix 4 of [18].

Author Contributions

A.A.: developed part of the new model set-up, most of the simulations, and the analysis as lead author; B.K.: developed the conceptualization and most of the model set-up, reviewed the model formulation, and collaborated in the analysis of the results; E.M.: reviewed the model formulation and the literature review and collaborated in the analysis of the results; L.R.R.-S.: reviewed the model formulation and hypothesis statement and collaborated in the analysis of the results and overall content as supervisor. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by Nazarbayev University Ph.D. scholarship of Bagdagul Kabdenova and M.Sc. scholarship of Assetbek Ashirbekov, funded by the Ministry of Education and Science of the Republic of Kazakhstan. This work is also funded within the project NU Faculty Development Competitive Research Grants 2018, “Simulation of CO2 flow in porous media using Lattice Boltzmann Model”, #090118FD5321.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in the displayed figures and tables within the article but can also be obtained on request from the corresponding author within a reasonable time-frame.

Acknowledgments

The results reported in this study were obtained by using a modified version of the open-source platform DL_MESO LBM package https://www.scd.stfc.ac.uk (accessed on 10 January 2021), which is used in this work as a base to build a new, improved model. The authors acknowledge M. Seaton for providing the original version of the code.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

χLocal value of the mass fraction
A ¯ Non-dimensional Helmholtz’s free energy per mole
V c Critical volume
R Gas constant
T T c
T r Dimensionless temperature
T c Critical temperature
Z c Critical compressibility factor
b 1 ,   b 2 , b 3 Dimensionless coefficients for Helmholtz free energy for the PT equation of state
1 Critical exponent
a 20 ,   a 21 Coefficients of Landau expansion
G i Ginzburg number
d 1 Coefficients of rectilinear diameter
α ,   β , γ Critical exponents
G f Body force

References

  1. Polikhronidi, N.; Batyrova, R.; Aliev, A.; Abdulagatov, I. Supercritical CO2: Properties and Technological Applications-A Review. J. Therm. Sci. 2019, 28, 394–430. [Google Scholar] [CrossRef]
  2. Brown, D. A Hot Dry Rock Geothermal Energy Concept Utilizing Supercritical CO2 Instead of Water. In Proceedings of the Twenty-Fifth Workshop on Geothermal Reservoir Engineering, SGP-TR-165, Stanford, CA, USA, 24–26 January 2000. [Google Scholar]
  3. Holdych, D.J.; Georgiadis, J.G.; Buckius, R.O. Hydrodynamic instabilities of near-critical CO2 flow in microchannels: Lattice Boltzmann simulation. Phys. Fluids 2004, 16, 1791–1802. [Google Scholar] [CrossRef]
  4. Huai, X.L.; Koyama, S.; Zhao, T.S. An experimental study of flow and heat transfer of supercritical carbon dioxide in multi-port mini channels under cooling conditions. Chem. Eng. Sci. 2005, 60, 3337–3345. [Google Scholar] [CrossRef]
  5. Pitla, S.S.; Robinson, D.M.; Groll, E.A.; Ramadhyani, S. Heat transfer from supercritical carbon dioxide in tube flow: A critical review. HVACR Res. 1998, 4, 281–301. [Google Scholar] [CrossRef]
  6. Katopodes, N. Free-Surface Flow. Computational Methods, Butterworth-Heineman, Oxford. 2019. Available online: https://doi.org/10.1016/b978-0-12-815485-4.00002-4 (accessed on 10 January 2021). [CrossRef]
  7. Osher, S.; Sethian, J.A. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys. 1988, 79, 12–49. [Google Scholar] [CrossRef] [Green Version]
  8. Shan, X.; Chen, H. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev. E 1993, 47, 1815. [Google Scholar] [CrossRef] [Green Version]
  9. Yuan, P.; Schaefer, L. Equations of state in a lattice Boltzmann model. Phys. Fluids 2006, 18, 042101. [Google Scholar] [CrossRef]
  10. Liu, H.; Zhang, Y.; Valocchi, A.J. Lattice Boltzmann simulation of immiscible fluid displacement in porous media: Homogeneous versus heterogeneous pore network. Phys. Fluids 2015, 27, 052103. [Google Scholar] [CrossRef] [Green Version]
  11. Liu, H.; Valocchi, A.J.; Werth, C.; Kang, Q.; Oostrom, M. Pore-scale simulation of liquid CO2 displacement of water using a two-phase lattice Boltzmann model. Adv. Water Resour. 2014, 73, 144–158. [Google Scholar] [CrossRef] [Green Version]
  12. Liu, H.; Kang, Q.; Leonardi, C.R.; Schmieschek, S.; Narváez, A.; Jones, B.D.; Williams, J.R.; Valocchi, A.J.; Harting, J. Multiphase lattice Boltzmann simulations for porous media applications: A review. Comput. Geosci. 2016, 20, 777–805. [Google Scholar] [CrossRef] [Green Version]
  13. Zhao, H.; Ning, Z.; Kang, Q.; Chen, L.; Zhao, T. Relative permeability of two immiscible fluids flowing through porous media determined by lattice Boltzmann method. Int. Commun. Heat Mass Transf. 2017, 85, 53–61. [Google Scholar] [CrossRef]
  14. Boek, E.S.; Venturoli, M. Lattice-Boltzmann studies of fluid flow in porous media with realistic rock geometries. Comput. Math. Appl. 2010, 59, 2305–2314. [Google Scholar] [CrossRef] [Green Version]
  15. Stiles, C.D.; Xue, Y. High density ratio lattice Boltzmann method simulations of multicomponent multiphase transport of H2O in air. Comput. Fluids 2016, 131, 81–90. [Google Scholar] [CrossRef]
  16. Ikeda, M.K.; Rao, P.R.; Schaefer, L.A. A thermal multicomponent lattice Boltzmann model. Comput. Fluids 2014, 101, 250–262. [Google Scholar] [CrossRef]
  17. Kupershtokh, A.L.; Medvedev, D.A.; Karpov, D.I. On equations of state in a lattice Boltzmann method. Comput. Math. Appl. 2009, 58, 965–974. [Google Scholar] [CrossRef] [Green Version]
  18. Kiselev, S.B. Cubic Crossover Equation of State. Fluid Phase Equilib. 1998, 147, 7–23. [Google Scholar] [CrossRef]
  19. Kiselev, S.B.; Ely, J.F. Generalized corresponding states model for bulk and interfacial properties in pure fluids and fluid mixtures. J. Chem. Phys. 2003, 119, 8645–8662. [Google Scholar] [CrossRef] [Green Version]
  20. Kiselev, S.B.; Ely, J.F. Generalized crossover description of the thermodynamic and transport properties in pure fluids. Fluid Phase Equilib. 2004, 222–223, 149–159. [Google Scholar] [CrossRef]
  21. Feyzi, F.; Seydi, M.; Alavi, F. Crossover Peng-Robinson equation of state with introduction of a new expression for the crossover function. Fluid Phase Equilib. 2010, 293, 251–260. [Google Scholar] [CrossRef]
  22. Feyzi, F.; Riazi, M.R.; Shaban, H.I.; Ghotbi, S. Improving cubic equations of state for heavy reservoir fluids and critical region. Chem. Eng. Commun. 1998, 167, 147–166. [Google Scholar] [CrossRef]
  23. Kabdenova, B.; Rojas-Solórzano, L.R.; Monaco, E. Lattice Boltzmann simulation of near/supercritical CO2 flow featuring a crossover formulation of the equation of state. Comput. Fluids 2020, 216, 104820. [Google Scholar] [CrossRef]
  24. Atykhan, M.; Kabdenova, B.; Monaco, E.; Rojas-Solórzano, L.R. Modeling Immiscible Fluid Displacement in a Porous Medium Using Lattice Boltzmann Method. Fluids 2021, 6, 89. [Google Scholar] [CrossRef]
  25. He, X.; Luo, L. Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation. Phys. Rev. E 1997, 56, 6811–6817. [Google Scholar] [CrossRef] [Green Version]
  26. Huang, H.; Sukop, M.; Lu, X.-Y. Multiphase Lattice Boltzman Methods Theory and Application; John Wiley & Sons, Ltd.: Sussex, UK, 2015. [Google Scholar]
  27. Martys, N.S.; Chen, H. Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top. 1996, 53, 743–750. [Google Scholar] [CrossRef]
  28. Shan, X. Analysis and reduction of the spurious current in a class of multiphase Lattice Boltzmann models. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 2006, 73, 6–9. [Google Scholar] [CrossRef]
  29. Monaco, E.; Brenner, G.; Luo, K.H. Numerical simulation of the collision of two microdroplets with a pseudopotential multiple-relaxation-time lattice Boltzmann model. Microfluid. Nanofluidics 2014, 16, 329–346. [Google Scholar] [CrossRef]
  30. Guo, Z.; Zheng, C.; Shi, B. Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top. 2002, 65, 6. [Google Scholar] [CrossRef] [PubMed]
  31. Landau, L.; Lifshitz, E.M. Statistical Physics; Pergamon Press Ltd.: Headington Hill Hall, Oxford, England, UK, 1980. [Google Scholar]
  32. Thermophysical Properties of Fluid Systems in NIST Chemistry Webbook, NIST Standard Reference Database No.69, National Institute of Standards and Technology, (n.d.). Available online: https://webbook.nist.gov (accessed on 10 October 2018).
  33. Li, Q.; Luo, K.H.; Li, X.J. Forcing scheme in pseudopotential lattice Boltzmann model for multiphase flows. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 2012, 86, 016709. [Google Scholar] [CrossRef] [Green Version]
  34. Huang, H.; Krafczyk, M.; Lu, X. Forcing term in single-phase and Shan-Chen-type multiphase lattice Boltzmann models. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 2011, 84, 046710. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Liu, Y.; Mutailipu, M.; Jiang, L.; Zhao, J.; Song, Y.; Chen, L. Interfacial Tension and Contact Angle Measurements for the Evaluation of CO2-Brine Two-Phase Flow Characteristics in Porous Media. Environ. Prog. Sustain. Energy 2015, 34, 1756–1762. [Google Scholar] [CrossRef]
  36. Fakhari, A.; Rahimian, M.H. Phase-field modeling by the method of lattice Boltzmann equations. Phys. Rev. E 2010, 81, 036707. [Google Scholar] [CrossRef]
  37. Suekane, T.; Soukawa, S.; Iwatani, S.; Tsushima, S.; Hirai, S. Behavior of supercritical CO2 injected into porous media containing water. Energy 2005, 30, 2370–2382. [Google Scholar] [CrossRef]
  38. Swift, M.; Osborn, W.; Yeomans, J. Lattice Boltzmann Simulation of Nonideal Fluids. Phys. Rev. Lett. 1995, 75, 3824–3827. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  39. Fakhari, A.; Li, Y.; Bolster, D.; Christensen, K.T. A phase-field lattice Boltzmann model for simulating multiphase flows in porous media: Application and comparison to experiments of CO2 sequestration at pore scale. Adv. Water Resour. 2018, 114, 119–134. [Google Scholar] [CrossRef]
  40. Zou, Q.; He, X. On pressure and velocity flow boundary conditions and bounceback for the lattice Boltzmann BGK model. Phys. Fluids 1997, 9, 1591–1598. [Google Scholar] [CrossRef] [Green Version]
  41. Fakhari, A.; Rahimian, M. Investigation of deformation and breakup of a moving droplet by the method of lattice Boltzmann equations. Int. J. Numer. Methods Fluids 2009, 64, 827–849. [Google Scholar] [CrossRef]
Figure 1. Discrete velocity model D2Q9.
Figure 1. Discrete velocity model D2Q9.
Fluids 06 00434 g001
Figure 2. 2-Dimensional static H2O droplet at equilibrium surrounded by CO2. Contours of the domain with a size of 201 × 201 LU. Density is shown in LU [23]. (a) Grid one: H2O drop in red surrounded by an empty region in blue. (b) Grid two: CO2 in red surrounding an empty region in blue. (c) Density profile at the cross-section of the domain taken at the horizontal centerline (y = 100 LU).
Figure 2. 2-Dimensional static H2O droplet at equilibrium surrounded by CO2. Contours of the domain with a size of 201 × 201 LU. Density is shown in LU [23]. (a) Grid one: H2O drop in red surrounded by an empty region in blue. (b) Grid two: CO2 in red surrounding an empty region in blue. (c) Density profile at the cross-section of the domain taken at the horizontal centerline (y = 100 LU).
Fluids 06 00434 g002
Figure 3. Water droplet contour for different contact angles obtained through adjusting the interaction parameter g   w a l l (a) θeq = 70. (b) θeq = 90°. (c) θeq = 120°. (d) θeq = 130°.
Figure 3. Water droplet contour for different contact angles obtained through adjusting the interaction parameter g   w a l l (a) θeq = 70. (b) θeq = 90°. (c) θeq = 120°. (d) θeq = 130°.
Fluids 06 00434 g003
Figure 4. CO2 bubble passing through a narrow channel with a density ratio of 4 at Gf = 4.8 × 10−5 at (a) initial stage; (b) 36,600 ts; (c) 107,800 ts.
Figure 4. CO2 bubble passing through a narrow channel with a density ratio of 4 at Gf = 4.8 × 10−5 at (a) initial stage; (b) 36,600 ts; (c) 107,800 ts.
Fluids 06 00434 g004
Figure 5. The bubble velocity profile at narrow and wide sections of the channel in a periodic pseudo-steady regime.
Figure 5. The bubble velocity profile at narrow and wide sections of the channel in a periodic pseudo-steady regime.
Fluids 06 00434 g005
Figure 6. CO2 penetration into the channel initially filled with water at 297,500 ts and contact angle 70°. (a) Gf = 3.5 × 10−5. (b) Gf = 4 × 10−5. (c) Gf = 6 × 10−5.
Figure 6. CO2 penetration into the channel initially filled with water at 297,500 ts and contact angle 70°. (a) Gf = 3.5 × 10−5. (b) Gf = 4 × 10−5. (c) Gf = 6 × 10−5.
Fluids 06 00434 g006
Figure 7. CO2 displacement process with insufficient body force of 2.5 × 10−5. The same color scheme is used, as in Figure 6. (a) Geometry of the porous media with porosity 0.78. (b) CO2 penetration at 113,400 ts.
Figure 7. CO2 displacement process with insufficient body force of 2.5 × 10−5. The same color scheme is used, as in Figure 6. (a) Geometry of the porous media with porosity 0.78. (b) CO2 penetration at 113,400 ts.
Fluids 06 00434 g007
Figure 8. CO2 displacement process at body force 4.5 × 10−5 with a stable displacement pattern. The same color scheme is used, as in Figure 6. (a) 63,700 ts. (b) 186,900 ts.
Figure 8. CO2 displacement process at body force 4.5 × 10−5 with a stable displacement pattern. The same color scheme is used, as in Figure 6. (a) 63,700 ts. (b) 186,900 ts.
Fluids 06 00434 g008
Figure 9. CO2 penetration pattern with a parabolical profile and front initially at 75 LU from pore network. Same color scheme is used, as in Figure 6. (a) 0 ts. (b) 3300 ts. (c) 6100 ts. (d) 9000 ts.
Figure 9. CO2 penetration pattern with a parabolical profile and front initially at 75 LU from pore network. Same color scheme is used, as in Figure 6. (a) 0 ts. (b) 3300 ts. (c) 6100 ts. (d) 9000 ts.
Fluids 06 00434 g009
Figure 10. CO2 penetration pattern with a flat front profile initially at 7 LU from pore network. The same color scheme is used, as in Figure 6. (a) 0 ts. (b) 1600 ts. (c) 4800 ts. (d) 7100 ts.
Figure 10. CO2 penetration pattern with a flat front profile initially at 7 LU from pore network. The same color scheme is used, as in Figure 6. (a) 0 ts. (b) 1600 ts. (c) 4800 ts. (d) 7100 ts.
Fluids 06 00434 g010
Figure 11. Line probe location (shown in black) in the grid-independence analysis of CO2 penetration LBM model.
Figure 11. Line probe location (shown in black) in the grid-independence analysis of CO2 penetration LBM model.
Fluids 06 00434 g011
Figure 12. Average velocity magnitude along the probe line (shown in Figure 11) in the CO2 penetration LBM model.
Figure 12. Average velocity magnitude along the probe line (shown in Figure 11) in the CO2 penetration LBM model.
Fluids 06 00434 g012
Table 1. System-dependent parameters for the crossover PR EoS [23].
Table 1. System-dependent parameters for the crossover PR EoS [23].
Physical Critical ParametersCO2H2O
T c K 304.1 647.1
P c M P a 7.37 22.06
Z c 0.3046 0.3031
Critical shift η c r i t 0.10077 0.24483
Crossover parameters
G i 0.23136 0.05897
d 1 9.1179 4.32597
v 1 0.02437 0.00299
m 0 1.0459 1.4392
a 20 15.8502   10.20517
a 21 0.5333 5.2527
Table 2. Critical temperature of CO2 and H2O in physical and LU units [23].
Table 2. Critical temperature of CO2 and H2O in physical and LU units [23].
Physical T c for CO2 Physical   T c   for   H 2 O T c   for   CO 2   in   LU T c   for   H 2 O   in   LU
304.1 K647.1 K0.034280.07292
Table 3. Values of g CO 2 , wall and g H 2 O , wall for different contact angles.
Table 3. Values of g CO 2 , wall and g H 2 O , wall for different contact angles.
θeq g CO 2 , wall g H 2 O , wall
70°0.2−0.2
90°00
120°−0.20.2
130°−0.30.3
Table 4. Time-space average CO2 flux over a probe line (Figure 11) integrated over 7100 timesteps.
Table 4. Time-space average CO2 flux over a probe line (Figure 11) integrated over 7100 timesteps.
Grid Size, LU2CO2 Flux, × 10−7 LU/tsRelative Error (%)
401 × 2013.590-
601 × 3013.6511.7%
801 × 4013.6790.76%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ashirbekov, A.; Kabdenova, B.; Monaco, E.; Rojas-Solórzano, L.R. Equation of State’s Crossover Enhancement of Pseudopotential Lattice Boltzmann Modeling of CO2 Flow in Homogeneous Porous Media. Fluids 2021, 6, 434. https://doi.org/10.3390/fluids6120434

AMA Style

Ashirbekov A, Kabdenova B, Monaco E, Rojas-Solórzano LR. Equation of State’s Crossover Enhancement of Pseudopotential Lattice Boltzmann Modeling of CO2 Flow in Homogeneous Porous Media. Fluids. 2021; 6(12):434. https://doi.org/10.3390/fluids6120434

Chicago/Turabian Style

Ashirbekov, Assetbek, Bagdagul Kabdenova, Ernesto Monaco, and Luis R. Rojas-Solórzano. 2021. "Equation of State’s Crossover Enhancement of Pseudopotential Lattice Boltzmann Modeling of CO2 Flow in Homogeneous Porous Media" Fluids 6, no. 12: 434. https://doi.org/10.3390/fluids6120434

APA Style

Ashirbekov, A., Kabdenova, B., Monaco, E., & Rojas-Solórzano, L. R. (2021). Equation of State’s Crossover Enhancement of Pseudopotential Lattice Boltzmann Modeling of CO2 Flow in Homogeneous Porous Media. Fluids, 6(12), 434. https://doi.org/10.3390/fluids6120434

Article Metrics

Back to TopTop