A Lattice Boltzmann Model of Multicomponent Fluids

The LB method of modeling fluid dynamics is actually a family [3] of models with varying degrees of faithfulness to the properties of real liquids. These methods are currently in a state of evolution as the models become better understood and are corrected for various deficiencies. The approach of LB is to consider a typical volume element of fluid to be composed of a collection of particles that are represented in terms of a particle velocity distribution function at each point in space. The particle velocity distribution, $n^{i}_a({\bf x},t)$, is the number density of particles at node x, time t, with velocity, e_a, where (a=1,...,b) indicates the velocity direction and superscript i labels the fluid component. The time is counted in discrete time steps, and the fluid particles can collide with each other as they move under applied forces.

For this paper we use the D3Q19 (3 Dimensional lattice with b=19) lattice [4, 5]. The microscopic velocity, e_a, equals all permutations of (±1, ±1, 0) for 1 a 12, (± 1, 0, 0) for 13 $1 \leq a \leq 12$, a $(\pm1, 0, 0)$ 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 the density, n_i(x,t), and the fluid velocity, uⁱ, of each fluid component, i, are obtained by taking suitable moment sums of $n^{i}_a({\bf x},t)$. Note that while the velocity 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)-g^{i}_a,$$ (1)

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

$$\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],$$ (2)

where $n^{i (eq)}_a({\bf x},t)$ is the equilibrium distribution at $({\bf x}, t)$ and $\tau_{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, 8]:

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

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

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},$$ (5)

S is the number of fluid component, m,ⁱ is the molecular mass of the ith component, t_a = 1/36 for 1 a 12, t_a = 1/18 for 13 a 18 and t₁₉ = 1/3. The free parameter d_o can be related to an effective temperature, T, for the system by the following moment of the equilibrium distribution:

$$T({\bf x},t) = \frac{\sum_{a} n^{i(eq)}_{a}({\bf x},t) ({\bf e}_{a}-{\bf v})^2} {3n^{i}({\bf x},t)},$$ (6)

which results in T=(1-d_o)/2 (we take units such that the Boltzmann constant k_b)=1).

The above formalism leads to a velocity field that is a solution of the Navier-Stokes [7] equation with the kinematic viscosity, $\nu=\frac{c^2}{6}(\sum_{i}^{S}c_{i} \tau_i -\frac{1}{2})$ where c_i is the concentration of each component [8].

Phase Separation of Fluids

There are a variety of approaches to modeling the phase separation of fluids [3, 9]. In the Shan-Chen model, a force, $\frac{d{\bf p}^i}{dt}({\bf x})$, between the two fluids is introduced that effectively perturbs the equilibrium velocity [4, 5] for each fluid so that they have a tendency to phase separate:

$$n^i({\bf x}){\bf v}^{'}({\bf x})=n^i {\bf v}({\bf x})+ \tau_i \frac{d{\bf p}^i}{dt}({\bf x})$$ (7)

where ${\bf v}'$ is the new velocity used in Eqs. [3] and [4]. A simple forcing that depends on the density of each fluid, is as follows [4, 5]:

$$\frac{d{\bf p}^i}{dt}({\bf x})=-n^i({\bf x})\sum_{i'}^ {S}\sum_{a}G_{ii'}^{a} n^{i'}({\bf x}+{\bf e}_a) {\bf e}_a$$ (8)

with $G_{ii'}^{a}= 2G$ for $\vert{\bf e}^a\vert = 1$; $G_{ii'}^{a}= G$ for $\vert{\bf e}^a\vert = \sqrt{2}$; and $G_{ii'}^{a}= 0$ for $i=i'$. $G$ is a constant that controls the strength of the interaction. Clearly, the forcing term is related to the density gradient of the fluid. It has been shown that the above forcing term can drive the phase separation process and naturally produce an interfacial surface tension effect consistent with the Laplace law boundary condition [5].

Phase separation of fluid can also be modeled by directly incorporating the force, $\frac{d{\bf p}^i}{dt}({\bf x})$, into the body force term. First note that in the continuum Boltzmann equation, the body force term is written ${\bf a} \cdot \nabla_e n( {\bf x},{\bf e})$, where ${\bf a}$ is an acceleration field due to a body force. A representation [10] of this body force term to second order in Hermite polynomials, in the discrete velocity space of the D3Q19 lattice is

$$g_a=-3t_an({\bf x})\left[({\bf e_a-v})\cdot { {\bf a}} +3({\bf e_a \cdot v)( e_a \cdot a)}\right]$$ (9)

To first order, the body force term is written as $g_a=-3t_in({\bf x}){\bf e_a}\cdot { \bf a}$.

In both models, phase separation takes place when the mutual diffusivity of the binary mixture becomes negative. An analytical expression for the mutual diffusivity has been determined in a previous work [8]. For the case of a critical composition, the condition for the system studied to undergo phase separation is $G \geq \frac{T}{12(n^1+n^2)}$.