Brief Review of Lattice Boltzmann Model of Fluid Phase Separation

The LB method of modeling fluid dynamics is actually a family 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 corrected for various deficiencies. In the present work, we utilize two LB models of complex fluids. The first and primarily studied method in this paper, was proposed by Shan and Chen [2,5]. It is particularly simple in form and adaptable to complex flow conditions like the presence of solid-fluid and air-fluid boundaries. For comparison, a second approach is studied [7] that directly incorporates the fluid/fluid interaction into a body-force term. This approach removes second order contributions with respect to the forcing (see below). The physical basis of such an approach is further described in the Appendix A. A third approach, not studied in this paper, is due to Yeomans et al. [3, 16] and is strongly rooted in the Cahn-Hilliard model of binary mixtures [12,13, 14,15,16]. A major criticism of this approach is the lack of Galilean invariance, but recent work suggests that the errors involved can be controlled in the description of the phase separation of fluids in the absence of shear [16]. It remains to be seen whether these higher order effects cause problems under conditions of fluid flow. All of these models, however, are not strictly energy conserving, which should be important under conditions of deep quenches and high rates of flow. Recent work [10] has introduced a framework for overcoming this technical problem and allows for systematic improvements of the LB calculation of fluid properties. In the present paper, we focus on isothermal flows and will examine the effect of energy conservation on the two-phase region of LB fluid mixtures in a future work.

We now present a brief description of the LB methods used in this study. A considerable computational advantage in modeling the fluid is obtained by restricting the particle positions to the sites of a lattice. The LB method extends the standard "quasi-crystalline" fluid model [19,20,21,22] of classical statistical mechanics by specifying the particle velocity distribution at each particle position so that both equilibrium and dynamical properties of the fluid can be calculated. Macroscopic variables such as density and fluid velocity are obtained by taking suitable moments of the velocity 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 number density of particles at node x, time t with velocity e_a where (a = 1,...,b). In this study, the particle velocities are represented in terms of a basis set defined on a cubic lattice. The velocity vectors are directed towards points that are nearest and next-nearest neighbors about a central lattice site and are in units of the lattice spacing divided by the time step. In the literature, this lattice is called a D3Q19 lattice where the 19 corresponds to the basis set size, b = 19 [23]. Time is counted in discrete time steps, and the fluid particles can collide with each other as they move under applied forces (surface tension, applied shear, etc.). The directions of the particle velocities are discretized in such a way that the isotropy of the pressure tensor and other properties are preserved [1]. The discretization of the time and spatial coordinates greatly reduces the information required to specify the fluid dynamics, thus extending the scale of the system that can be modeled by a given computational resource. This is the computational virtue of a lattice formulation. The price is that care must be taken to avoid artifacts that arise from the lattice model.

Macroscopic quantities such as density, $\rho^i({\bf x},t)$ ⁱ (x,t), and fluid velocity, uⁱ, of each fluid component, i, are obtained by the following moment sums,

$$\rho^i({\bf x},t)=m^in^{i}({\bf x},t) = m^i\sum_{a} n^{i}_a({\bf x},t),$$ (1)

and

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

where mⁱ is the molecular mass of the ith component. While the distribution function is defined only over a discrete set of velocities, the actual macroscopic velocity field of the fluid is nearly 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),$$ (3)

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 [23,24],

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

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 [23]:

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

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

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

and the weights are t_a=1/36 for 1 a 12 , t_a=1/36 for 13 a $13 \leq a \leq 18$ 18, and t₁₉=1/3 . The parameter d ⁱ_o can be related by self-consistency to an effective temperature, T, by the following moment of the equilibrium distribution,

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

In order that both fluids components have the same temperature, d ⁱ_o may be defined by the relation $d^i_o=1-2\frac{T}{m^i}$, where we choose units here such that the Boltzmann's constant k_B equals 1.

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

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

where c_i is the concentration of each component and the lattice constant c² = 2 for this model.