Modeling Phase Separating Binary Mixtures in Porous Media

The modeling of phase separating binary mixtures in porous media is a great research challenge. Recently, new cellular automata methods called lattice gas and lattice Boltzmann (LB) [4] have demonstrated the capability to model multiphase fluid flow while including interfacial surface tension effects between the two fluids and between the fluid and solid interface [2]. Lattice Boltzmann and lattice gas algorithms are also easily adaptable for parallel computers since they, in general, only need nearest neighbor information on the lattice on which they are defined. Such information is readily available in a digital image of a porous medium.

In this section we present a brief description of the LB method used in this study. A detailed description of the LB method is found in the paper by Martys and Chen [2]. The approach of LB is to consider a typical volume element of fluid to be composed of a collection of particles which are represented in terms of a particle velocity distribution at each point in space. Macroscopic variables such as density and fluid velocity can be obtained by taking suitable moments of the distribution function. The velocity distribution function, $n^{i}_a({\bf x},t)$, where superscript i labels the fluid component and the subscript a indicates the velocity direction, is the amount of particles at node x, time t with velocity e_a where (a=1,...,b). The time is counted in discrete time steps, and the fluid particles can colide with each other as they move under applied forces (surface tension, applied shear, etc.). The directions of the particle velocities are discretized reflecting the lattice and physical symmetries.

For this study we use the D3Q19 (3 Dimensional lattice with b=19) [5] lattice where the discrete particle velocities, e_a , equal all permutations of (±1, ±1, 0) for 1 a 12, (±1 , 0, 0) for 13 a 18 and (0,0,0) for a=19. The units of e_a are the lattice constant divided by the time step. Macroscopic quantities such as density, n ⁱ(x,t ), and fluid velocity, u ⁱ, of each fluid component, i, are obtained by the following moment sums,

$$n^{i}({\bf x},t) = \sum_{a} n^{i}_a({\bf x},t),$$ (3)

and

$${\bf u}^{i}({\bf x},t) = \frac{\sum_{a} n^{i}_{a}({\bf x},t){\bf e}_{a}} {n^{i}({\bf x},t)}.$$ (4)

Note that while the distribution function is defined only over a discrete set of velocities, the actual macroscopic velocity field of the fluid is continuous.

The time evolution of the particle velocity distribution function satisfies the following LB equation:

$$n^{i}_a({\bf x}+{\bf e}_a,t+1)-n^{i}_a({\bf x},t) = \Omega^{i}_a({\bf x},t),$$ (5)

where $\Omega^{i}_a$ is the collision operator representing the rate of change of the particle distribution due to collisions. The collision operator is greatly simplified by use of the single time relaxation approximation [6,5]

$$\Omega^{i}_a({\bf x},t) =-\frac{1}{\tau}_{i} \left[n^{i}_a({\bf x},t)-n^{i (eq)}_a({\bf x},t)\right],$$ (6)

where $n^{i (eq)}_a({\bf x},t)$ is the equilibrium distribution at (x, t ) and $\tau_{i}$ _i is the relaxation time that controls the rate of approach to equilibrium. The equilibrium distribution can be represented in the following form for particles of each type [5]:

$$n^{i(eq)}_a({\bf x})=t_an^{i}({\bf x})\left[1+ 3{\bf e}_a \c... ...ac{3}{2}({3\bf e}_a{\bf e}_a:{\bf v}{\bf v}-{\bf v}^2)\right]$$ (7)

$$n^{i(eq)}_{19}({\bf x})=t_{19}n^{i}({\bf x})\left[1- \frac{3}{2}{\bf v}^2\right],$$ (8)

where

$${\bf v}=\frac{\sum_{i}^{S}m^{i}\sum_{a} n^{i}_{a}{\bf e}_{a}/ \tau_i} {\sum_{i}^{S}m^{i} n^{i}({\bf x})/\tau_i},$$ (9)

mⁱ is molecular mass of the ith component, and t_a=1/36 for $1 \leq a \leq 12$ and t_a=1/36 for $13 \leq a \leq 18$ and t₁₉=1/3 .

It has been shown that the above formalism leads to a velocity field which is a solution of the Navier-Stokes [6] equation with the kinematic viscosity [2], $\nu$,

$$\nu=c^2\frac{\sum_{i}^{S}c_{i} \tau_i -\frac{1}{2}}{6}$$ (10)

where c_i is the concentration of each component.