Comparative investigation of a lattice Boltzmann boundary treatment of multiphase mass transport with heterogeneous chemical reactions

Multiphase reactive transport in porous media is an important component of many natural and engineering processes. In the present study, boundary schemes for the continuum species transport–lattice Boltzmann (CST-LB) mass transport model and the multicomponent pseudopotential model are proposed to simulate heterogeneous chemical reactions in a multiphase system. For the CST-LB model, a lattice-interface-tracking scheme for the heterogeneous chemical reaction boundary is provided. Meanwhile, a local-average virtual density boundary scheme for the multicomponent pseudopotential model is formulated based on the work of Li et al. [Li, Yu, and Luo, Phys. Rev. E 100, 053313 (2019)]. With these boundary treatments, a numerical implementation is put forward that couples the multiphase fluid flow, interfacial species transport, heterogeneous chemical reactions, and porous matrix structural evolution. A series of comparison benchmark cases are investigated to evaluate the numerical performance for different pseudopotential wetting boundary treatments, and an application case of multiphase dissolution in porous media is conducted to validate the present models’ ability to solve complex problems. By applying the present LB models with reasonable boundary treatments, multiphase reactive transport in various natural or engineering scenarios can be simulated accurately.

Figure 1

A schematic diagram of the three interfacial mass transport treatment schemes: (a) the lattice-interface-tracking scheme, (b) the physical-interface-tracking scheme, and (c) the phase-fraction-indicating scheme. BC denotes the position where the boundary condition is enforced, and the numbers in (c) represent the phase fraction of each node. The blue lines denote the numerical interface shape under the corresponding interfacial scheme.

Figure 2

A schematic illustration of the heterogeneous chemical reaction boundary treatment [45]. $g_2^{*}$ can be calculated using Eq. (26) and is implemented in the streaming step of the LB algorithm.

Figure 3

A schematic diagram of the overall numerical implementation.

Figure 4

The initial phase distribution of the semicircular droplet covering the solid cylinder.

Figure 5

The phase fraction $x_{\mathrm{A}}$ distribution simulated with the multicomponent pseudopotential model using the local-average virtual density boundary scheme with (a) $\mathrm{\Delta }{\rho }_{\mathrm{A}}=1.0,{\lambda }_{\mathrm{B}}=1.5$; (b) ${\lambda }_{\mathrm{A}}=1.0,{\lambda }_{\mathrm{B}}=1.0$; and (c) ${\lambda }_{\mathrm{A}}=2.5,{\lambda }_{\mathrm{B}}=1.0$; corresponding to contact angles (a) $\theta \approx {124}^{\circ }$, (b) $\theta \approx {85}^{\circ }$, and (c) $\theta \approx {30}^{\circ }$.

Figure 6

The phase fraction $x_{\mathrm{A}}$ distribution simulated by the multicomponent pseudopotential model using the solid-fluid interaction boundary scheme with (a) $G_{\mathrm{As}}=-3.0,G_{\mathrm{Bs}}=0.0$; (b) $G_{\mathrm{As}}=-4.5,G_{\mathrm{Bs}}=0.5$; and (c) $G_{\mathrm{As}}=-6.0,G_{\mathrm{Bs}}=1.0$; corresponding to contact angles (a) $\theta \approx {127}^{\circ }$, (b) $\theta \approx {81}^{\circ }$, and (c) $\theta \approx {27}^{\circ }$.

Figure 7

The phase fraction $x_{\mathrm{A}}$ distribution simulated by the multicomponent pseudopotential model using the virtual density boundary scheme with (a) ${\rho }_{\mathrm{sA}}=1.0,{\rho }_{\mathrm{sB}}=0.5$; (b) ${\rho }_{\mathrm{sA}}=3.0,{\rho }_{\mathrm{sB}}=0.3$; and (c) ${\rho }_{\mathrm{sA}}=5.0,{\rho }_{\mathrm{sB}}=0.1$; corresponding to the contact angles (a) $\theta \approx {134}^{\circ }$, (b) $\theta \approx {86}^{\circ }$, and (c) $\theta \approx {26}^{\circ }$.

Figure 8

The phase fraction $x_{\mathrm{A}}$ distribution along the vertical central axis of the computational domain with (a) the local-average virtual density scheme, (b) the solid-fluid interaction scheme, and (c) the global virtual density scheme.

Figure 9

The schematic diagram of the physical problem in benchmark 1, one-dimensional multiphase diffusion with heterogeneous chemical reactions.

Figure 10

A comparison of the analytical solutions and numerical results for the concentration distribution of $P_1$ at $t=0.5\mathrm{s}$ for benchmark 1 ($H=0.5$). The numerical results were generated for (a) the local-average virtual density scheme, (b) the solid-fluid interaction scheme, and (c) the global virtual density scheme.

Figure 11

A comparison of the analytical solutions and numerical results for the concentration distribution of $P_1$ at $t=0.5\mathrm{s}$ for benchmark 1 ($H=0.2$). The numerical results were generated for (a) the local-average virtual density scheme, (b) the solid-fluid interaction scheme, and (c) the global virtual density scheme.

Figure 12

Relative errors under different Henry coefficients using the three wetting boundaries in the multiphase flow model with the large contact angle.

Figure 13

The evolution of the $P_1$ concentration and solid structure for benchmark 2 with contact angles (a) $\theta \approx {124}^{\circ }$, (b) $\theta \approx {85}^{\circ }$, (c) $\theta \approx {30}^{\circ }$ using the local-average virtual density scheme. The red line denotes the interface between fluid A and fluid B. The gray color represents the solid structure, and the black dotted line denotes the initial solid profile.

Figure 14

The evolution of $P_1$ concentration and solid structure for benchmark 2 with contact angle $\theta \approx {81}^{\circ }$ using the solid-fluid interaction scheme.

Figure 15

The evolution of $P_1$ concentration and solid structure for benchmark 2 with contact angle $\theta \approx {86}^{\circ }$ using the global virtual density scheme.

Figure 16

The absolute value of the spurious current produced by (a) the local-average virtual density scheme, (b) the solid-fluid interaction scheme, and (c) the global virtual density scheme. The red arrows denote the velocity vectors.

Figure 17

The profile of concentration along the red dotted line in the contour map at $t=0.01\mathrm{s}$ for the single-phase and multiphase problems in benchmark 2.

Figure 18

The computational domain for the reactive transport problem during drainage in a capillary tube.

Figure 19

The concentration evolution of the reactive transport during drainage in a capillary tube with $\mathrm{Da}=3.0$, $\mathrm{Pe}=6.0$, and $\mathrm{Ca}=1.75\times {10}^{-3}$. The numerical results for the CST-LB and VOF-CST methods are compared.

Figure 20

The concentration distribution along the horizontal central axis in Fig. 19.

Figure 21

The concentration evolution of the reactive transport during drainage in a capillary tube with $\mathrm{Da}=3.0$, $\mathrm{Pe}=6.0$, and $\mathrm{Ca}=1.75\times {10}^{-2}$.

Figure 22

The concentration distribution along the horizontal central axis in Fig. 21.

Figure 23

The initial phase distribution for the application case with fluid A saturations $S_{\mathrm{A}}=0.28$ and $0.55$. The black color represents the unreactive rock matrix, and the gray color the soluble solid. The red color represents fluid A and the blue color fluid B.

Figure 24

The evolution of the $P_1$ concentration for $S_{\mathrm{A}}=0.55$, $\mathrm{Pe}=0.07$, and $\mathrm{Da}=0.64$, computed with the local-average virtual density scheme. The red line represents the phase interfaces.

Figure 25

The evolution of the $P_1$ concentration for $S_{\mathrm{A}}=0.28$, $\mathrm{Pe}=2.49$, and $\mathrm{Da}=0.64$, computed with the local-average virtual density scheme. The red line represents the phase interfaces.

Figure 26

A comparison of the dissolution ratios and average concentration curves with and without the mass transport limitation for two $S_{\mathrm{A}}$ cases using the local-average virtual density scheme.

Figure 27

The concentration distribution of the reactive transport during drainage at $t=0.05\mathrm{s}$ for (a) the local-average virtual density scheme, (b) the solid-fluid interaction scheme, and (c) the global virtual density scheme. (d) The concentration profile at the vertical section at $x=40\mu \mathrm{m}$ for different numerical treatments.

Figure 28

The evolution of the $P_1$ concentration for $S_{\mathrm{A}}=0.55$, $\mathrm{Pe}=0.07$, and $\mathrm{Da}=0.64$, computed with the global virtual density scheme. The red line represents the phase interfaces. The unphysical dissolution in fluid B is marked by yellow circles.

Figure 29

A comparison of average concentration curves with the three wetting boundary treatments.