Lattice Boltzmann method for diffusion-limited partial dissolution of fluids

A lattice Boltzmann model for two partially miscible fluids is developed. By partially miscible we mean that, although there is a definite interfacial region separating the two fluids with a surface tension force acting at all points of the transition region, each fluid can nonetheless accept molecules from the other fluid up to a set solubility limit. We allow each fluid to diffuse into the other with the solubility and diffusivity in each fluid being input parameters. The approach is to define two regions within the fluid: one interfacial region having finite width, across which most of the concentration change occurs, and in which a surface tension force and color separation step are allowed for and one miscible fluid region where the concentration of the binary fluids follows an advection-diffusion equation and the mixture as a whole obeys the Navier-Stokes incompressible flow equations. Numerical examples are presented in which the algorithm produces results that are quantitatively compared to exact analytical results as well as qualitatively examined for their reasonableness. The model has the ability to simulate how bubbles of one fluid flow through another while dissolving their contents as well as to simulate a range of practical invasion problems such as injecting supercritical ${\mathrm{CO}}_2$ into a porous material saturated with water for sequestration purposes.

Figure 1

Schematic representation of how the separate algorithm regions are defined. The figure shows the concentration profile of a circular blue bubble surrounded by red fluid.

Figure 2

Interface position $s$, where $\phi ={\alpha }_2=0.02$, as a function of time $t$. The position and time are measured, respectively, in terms of the lattice constant $\mathrm{\Delta }x$ and the time discretization $\mathrm{\Delta }t$. Open symbols show the numerical results, while black lines represent the analytical solutions. We have used linear interpolation to estimate subgrid positions of the interface.

Figure 3

Concentration profiles for $D_0=0.1$ and ${\alpha }_2=0.02$ as a function of position. Open circles represent the numerical results. Solid black lines show the analytical solution. Vertical dashed lines indicate, from the left, the position of the fluid interface at times $t=[57.0\times {10}^3,226\times {10}^3,505\times {10}^3,900\times {10}^3]$, respectively. The horizontal dashed line marks the saturation concentration ${\alpha }_2$.

Figure 4

Concentration profiles of a blue fluid bubble at different times during evaporation. The concentrations $\phi$ are given as functions of position measured in terms of the lattice constant $\mathrm{\Delta }x$. In this simulation the solubility levels are given by ${\alpha }_1=0.95$ and ${\alpha }_2=0.02$. The bubble radius was here initiated to be $r_0=20$ lattice units. The concentration profiles are plotted at time steps $t^{\prime }=t-t_0=[0.0,1.0\times {10}^3,1.0\times {10}^4,5.0\times {10}^4,1.0\times {10}^5,5.0\times {10}^5]$, where $t_0={10}^4$ was the initiation time for the evaporation process. To better see the behavior of the diffusion process taking place, two separate regions of the concentration profiles have been cut out and magnified. The initial and final concentration profiles are shown in solid lines.

Figure 5

Concentration profiles of an evaporating blue fluid bubble at different times as functions of position. Here ${\alpha }_1=0.8,{\alpha }_2=0.2$, and the initial bubble radius $r_0=25$ lattice units. The concentration profiles are plotted for the same time steps as in Fig. 4. The initial and final concentration profiles are as in Fig. 4 represented by solid lines.

Figure 6

Pressure profiles of an evaporating bubble as functions of position. This figure shows the same system as Fig. 5 and the profiles plotted are at the same time steps as in Figs. 4 and 5. The pressure difference $\mathrm{\Delta }p$ is the difference between the measured pressure at position $x$ and the pressure measured at a point outside the bubble, far away from the interface region between the two phases. The pressure is given in lattice units.

Figure 7

Snapshot of an invading fluid being injected from the left side of a straight channel, while it partially dissolves in a defending fluid. The figure shows the invading fluid concentration at $t=7.5\times {10}^4$ with $\mathrm{Pe}=3\times {10}^{-1}$ and $\mathrm{Ca}={10}^{-2}$ ($J_x={10}^{-4},\sigma ={10}^{-3}$, and $D_0={10}^{-2}$). The color scale (grayscale) is a linear scale showing the invading fluid concentration $\phi$. The thin white line represents the interface between the two fluids.

Figure 8

Snapshot of the invading fluid concentration $\phi$ at $t=1.5\times {10}^5$ with $\mathrm{Pe}=3$ and $\mathrm{Ca}={10}^{-2}$ ($J_x={10}^{-4},\sigma ={10}^{-3}$, and $D_0={10}^{-3}$). See Fig. 7 for further description.

Figure 9

Snapshot of the invading fluid concentration $\phi$ at $t=3.0\times {10}^4$ with $\mathrm{Pe}=30$ and $\mathrm{Ca}={10}^{-1}$ ($J_x={10}^{-3},\sigma ={10}^{-3}$, and $D_0={10}^{-3}$). The capillary number was increased by increasing the rate of injection. See Fig. 7 for further description.

Figure 10

For a capillary number $\mathrm{Ca}={10}^{-2}$, two invasion processes, where partial fluid dissolution occurs, are compared at different times. In both processes, an invading fluid is, while being injected from the left side of a porous medium, dissolving into the defending fluid. The porous structure is shown as white areas outlined in black. Each column shows one experiment where the images are ordered, from the top, by increasing time. The two pictures in each row depict the same instants in time. The invasion structure, i.e., the pure invading fluid, is presented in blue (dark gray), extending from the left side of the system. The invading fluid dissolved in the defending one is represented by a linear color scale (grayscale) given below the last image in each column. The darkest shade represents $0.5\ge \phi \ge {\alpha }_2=0.05$, while the lightest shade represents $\phi =0$. The left column shows the structure evolution of an invasion process characterized by $\mathrm{Pe}=7\times {10}^{-1}$ ($J_x={10}^{-3},\nu ={10}^{-1},\sigma ={10}^{-2}$, and $D_0={10}^{-2}$). In the right column, a similar process with $\mathrm{Pe}=7\times {10}^{-2}$ ($J_x={10}^{-3},\nu ={10}^{-1},\sigma ={10}^{-2}$, and $D_0={10}^{-1}$) is presented.