Generalized three-dimensional lattice Boltzmann color-gradient method for immiscible two-phase pore-scale imbibition and drainage in porous media

This article presents a three-dimensional numerical framework for the simulation of fluid-fluid immiscible compounds in complex geometries, based on the multiple-relaxation-time lattice Boltzmann method to model the fluid dynamics and the color-gradient approach to model multicomponent flow interaction. New lattice weights for the lattices D3Q15, D3Q19, and D3Q27 that improve the Galilean invariance of the color-gradient model as well as for modeling the interfacial tension are derived and provided in the Appendix. The presented method proposes in particular an approach to model the interaction between the fluid compound and the solid, and to maintain a precise contact angle between the two-component interface and the wall. Contrarily to previous approaches proposed in the literature, this method yields accurate solutions even in complex geometries and does not suffer from numerical artifacts like nonphysical mass transfer along the solid wall, which is crucial for modeling imbibition-type problems. The article also proposes an approach to model inflow and outflow boundaries with the color-gradient method by generalizing the regularized boundary conditions. The numerical framework is first validated for three-dimensional (3D) stationary state (Jurin's law) and time-dependent (Washburn's law and capillary waves) problems. Then, the usefulness of the method for practical problems of pore-scale flow imbibition and drainage in porous media is demonstrated. Through the simulation of nonwetting displacement in two-dimensional random porous media networks, we show that the model properly reproduces three main invasion regimes (stable displacement, capillary fingering, and viscous fingering) as well as the saturating zone transition between these regimes. Finally, the ability to simulate immiscible two-component flow imbibition and drainage is validated, with excellent results, by numerical simulations in a Berea sandstone, a frequently used benchmark case used in this field, using a complex geometry that originates from a 3D scan of a porous sandstone. The methods presented in this article were implemented in the open-source PALABOS library, a general C++ matrix-based library well adapted for massive fluid flow parallel computation.

Figure 1

A boundary node example for the D2Q9 lattice. The dashed arrows represent the unknown populations.

Figure 2

Dimensionless quantity $\frac{h}{R}\mathrm{Bo}$ as function of the viscosity ratio ${\nu }_B/{\nu }_T$, the Laplace number $\mathrm{La}$, and the contact angle ${\theta }_c$ compared to the analytical two-dimensional Jurin law. Note that the lattice resolution $s=2$.

Figure 3

Dimensionless quantity $\frac{h}{R}\mathrm{Bo}$ as function of the contact angle ${\theta }_c$ and the lattice resolution $s$ compared to the analytical two-dimensional Jurin law. Note that the viscosity ratio ${\nu }_B/{\nu }_T=20$ and the Laplace number $\mathrm{La}=10$.

Figure 4

Dimensionless quantity $\frac{h}{R}\mathrm{Bo}$ as function of the viscosity ratio ${\nu }_B/{\nu }_T$, the Laplace number $\mathrm{La}$, and the contact angle ${\theta }_c$ compared to the analytical three-dimensional Jurin law. Note that the lattice resolution $s=2$.

Figure 5

Graphical representation of the rotated tube geometry with the bottom fluid highlight. As predicted by Jurin's law, the capillary rise is clearly visible. Note that the lattice resolution $s=2$, the viscosity ratio ${\nu }_B/{\nu }_T=20$, the Laplace number $\mathrm{La}=10$, the contact angle ${\theta }_c={45}^{\circ }$, and the geometry rotations $({\theta }_x,{\theta }_y,{\theta }_z)=({30}^{\circ },{30}^{\circ },{30}^{\circ })$.

Figure 6

Dimensionless quantity $\frac{h}{R}\mathrm{Bo}$ as function of the contact angle ${\theta }_c$ and the geometry rotations $({\theta }_x,{\theta }_y,{\theta }_z)$ compared to the analytical three-dimensional Jurin law. Note that the lattice resolution $s=2$, the viscosity ratio ${\nu }_B/{\nu }_T=20$, and the Laplace number $\mathrm{La}=10$.

Figure 7

Capillary intrusion of water into a single cylindrical pore initially filled with hexane, i.e., ${\rho }_I/{\rho }_O=1.5094$ and ${\mu }_I/{\mu }_O=3.2277$. The water-solid contact angle is ${40}^{\circ }$ and $t^{*}\approx 48.3$.

Figure 8

Dimensionless position of the fluid interface $l^{*}$ into the single cylindrical pore as function of the dimensionless time $t^{*}$ and the contact angle ${\theta }_c$. The numerical results of the proposed lattice Boltzmann color-gradient method are compared to the analytical three-dimensional Washburn law. Water-hexane density and viscosity ratios are simulated, i.e., ${\rho }_I/{\rho }_O=1.5094$ and ${\mu }_I/{\mu }_O=3.2277$. These simulations correspond to a capillary intrusion of water into a single cylindrical pore filled initially with hexane.

Figure 9

Dimensionless position of the fluid interface $l^{*}$ into the single cylindrical pore as function of the dimensionless time $t^{*}$ and the contact angle ${\theta }_c$. The numerical results of the proposed lattice Boltzmann color-gradient method are compared to the analytical three-dimensional Washburn law. Water-trichloroethylene density and viscosity ratios are simulated, i.e., ${\rho }_I/{\rho }_O=0.6828$ and ${\mu }_I/{\mu }_O=1.7583$. These simulations correspond to a capillary intrusion of water into a single cylindrical pore filled initially with trichloroethylene.

Figure 10

Dimensionless position of the fluid interface $l^{*}$ into the single cylindrical pore as function of the dimensionless time $t^{*}$ and the contact angle ${\theta }_c$. The numerical results of the proposed CGM are compared to the analytical three-dimensional Washburn law. The density and viscosity ratios are respectively, i.e., ${\rho }_I/{\rho }_O=10$ and ${\mu }_I/{\mu }_O=30$.

Figure 11

Real part of the oscillatory Reynolds number $\Re \left(\mathrm{Re}\right)$ as function of the relative viscosity $a$, the relative density $b$, and the dimensionless number $\mathrm{\Gamma }$; the quantity $b=1/10$ in the left plot and $a=1/10$ in the right plot.

Figure 12

Lenormand's phase diagram and the CGM simulation parameters corresponding to the “cross” and “plus” symbols. The “circle,” “diamond,” and “triangle” symbols correspond to parameters where an animation of the simulation is available in the Supplemental Material [94]. Note that the three regime Lenormand plateaus (viscous fingering, stable displacement, and capillary fingering) are illustrated only for qualitative comparison as the exact position of these plateaus depends on many factors [55].

Figure 13

Displacement regimes in random porous media: (a) viscous fingering, (b) stable displacement, and (c) capillary fingering. An animation of these simulations is available in the Supplemental Material [94].

Figure 14

Saturation $S$ of random porous media as function of the dynamic viscosity ratio $M$. The jump transition from viscous to stable regime is qualitatively captured by the CGM.

Figure 15

Saturation $S$ of random porous media as function of the capillary number $C$. The jump transition from capillary to stable regime is qualitatively captured by the CGM.

Figure 16

Pore matrix of the Berea sandstone [85]. The full pore matrix is discretized in a lattice of size $400\times 400\times 400$. This corresponds to a resolution of 5.3 $\mu {\text{m}}^{\text{3}}$ per voxel.

Figure 17

Saturation of the Berea sandstone as function of the stone depth $x^{*}$ in the viscous fingering regime. Imbibition (${\theta }_c={45}^{\circ }$) and drainage (${\theta }_c={135}^{\circ }$) are illustrated as well as the corresponding final global saturation of the sample.

Figure 18

Final view of a drainage CGM simulation in the viscous fingering regime with half the solid Berea sandstone matrix in the back. An animation of this simulation for a reduced size matrix is available in the Supplemental Material [94].

Figure 19

Final view of an imbibition CGM simulation in the viscous fingering regime with half the solid Berea sandstone matrix in the back. An animation of this simulation for a reduced size matrix is available in the Supplemental Material [94].

Figure 20

Saturation of the Berea sandstone as function of the stone depth $x^{*}$ in the stable displacement regime. Imbibition (${\theta }_c={45}^{\circ }$) and drainage (${\theta }_c={135}^{\circ }$) are illustrated as well as the corresponding final global saturation of the sample.

Figure 21

Logarithm of the saturation of the Berea sandstone as function of the stone depth $x^{*}$ in the stable displacement regime and drainage (${\theta }_c={135}^{\circ }$).

Figure 22

Final view of a drainage CGM simulation in the stable displacement regime with half the solid Berea sandstone matrix in the back (lattice size of $440\times 400\times 400$). It is possible to see in the black circle the isolated droplet that has reached the sample outlet before the main advancing front.

Figure 23

Saturation of the Berea sandstone as function of the stone depth $x^{*}$ in the capillary fingering regime. Imbibition (${\theta }_c={45}^{\circ }$) and drainage (${\theta }_c={135}^{\circ }$) are illustrated as well as the corresponding final global saturation of the sample.

Figure 24

Saturation of the Berea sandstone as function of the stone depth $x^{*}$ and the dimensionless quantity $\mathcal{L}$. The critical exit saturation threshold for percolation is also illustrated. Only the cases with $\mathcal{L}=0$ and $\mathcal{L}=0.5$ percolate the core sample under the chosen steady state simulation criterion.

Figure 25

Global saturation of the Berea sandstone as function of the time (in lattice units) and the dimensionless quantity $\mathcal{L}$.

Figure 26

Saturation of the Berea sandstone as function of the stone depth $x^{*}$ and the type of wetting boundary conditions. The dimensionless quantity $\mathcal{L}=10$ which means that there is a strong inverse pressure gradient. The standard wetting boundary conditions do not behave as expected with such a strong inverse pressure gradient.