PHYSICAL REVIEW E 88, 063305 (2013)
Multicomponent lattice Boltzmann equation method with a discontinuous hydrodynamic interface T. J. Spencer and I. Halliday* Materials and Engineering Research Institute, Sheffield Hallam University, Sheffield S1 1WB, United Kingdom (Received 4 October 2012; revised manuscript received 5 August 2013; published 6 December 2013) In the multicomponent lattice Boltzmann equation simulation method (MCLB), applied to the continuum regime of fluid flow, the finite width of the fluid-fluid interface introduces unphysical scales. We present a practical, robust, computationally efficient, and easy to implement solution to this problem which needs only low order interpolation to be stable and accurate and is applicable to any MCLB variant which uses a continuous phase field to distinguish between immiscible fluids with arrested coalescence. Our method extends the ideas of Kim and Pitsch, [Phys. Fluids 19, 108101 (2007)] and uses no external force distribution whatsoever to generate continuum interfacial physics, i.e., the Laplace law and no traction conditions on interfacial stresses. As such, it is amenable to the simplest form of Chapman-Enskog analysis used for lattice Boltzmann models. We assess our method and proceed to compare key results obtained with it against other equivalent data, obtained using the established continuum regime MCLB technique based upon the work of Lishchuk, Care, and Halliday, [Phys. Rev. E 67, 036701 (2003)] and Halliday, Hollis, and Care, [Phys. Rev. E 76, 026708 (2007)], quantifying performance in terms of the minimum feasible capillary available to simulation using that technique. DOI: 10.1103/PhysRevE.88.063305
PACS number(s): 47.11.−j
I. INTRODUCTION
Lattice Boltzmann method (LB) has advantages over conventional computational fluid dynamics in only a few areas, one being multicomponent flow. A facility with multicomponent flow is, alone, sufficient to justify the existence of the LB technique, however, for the flow of immiscible fluids with arrested coalescence, in the continuum hydrodynamic regime, it is of significance in a range of biological, microfluidic, and general engineering applications. The multicomponent lattice Boltzmann simulation method (MCLB) is a simple Eulerian numerical technique, in which a fluid-fluid interface is represented by a phase field, or order parameter, throughout denoted ρ N . Here we are concerned with the particular class of MCLB aimed at the continuum regime of hydrodynamics which uses a form of immersed boundary force [1], the action of which is directed by ρ N , to embed appropriate interfacial physics. The fact that ρ N cannot change discontinuously imposes extraneous length scales on MCLB simulations, characteristic of the distance over which the phase field changes between the values which represent the separated components. In mesoscopic applications, this is often acceptable but in large scale applications to continuum hydrodynamics it is inconvenient (since computational mesh is a scarce resource) and inaccurate (since the fluid-fluid interface is viewed as an infinitely thin surface, separating two fluids, upon which certain boundary conditions must be imposed [2]). In short, a diffuse hydrodynamic boundary obstructs detailed comparison with analytic calculations and calculations obtained using other numerical methods (such as boundary integral calculations), it introduces unphysical scales into simulations, and it is a drain on those computational resources which are associated with resolution. To address this problem we recount a means of imposing “discontinuous” interfaces between completely immiscible
* Corresponding author: [email protected]; FAX: (+44)114 225 3501.
1539-3755/2013/88(6)/063305(11)
LB fluids. For present purposes, we define a “discontinuous hydrodynamic interface” to be a closed, reducible, sublattice contour at which the separated fluids’ physical properties, such as pressure (or density) and viscosity but not, note, the phase field ρ N (which must remain differentiable) undergo a step change. The MCLB method we devise to achieve this is simple; it needs only low order interpolation to be stable and, judged in terms of the interfacial microcurrent (see below), the solution we advance introduces no penalty over the range of envisaged capillary numbers Ca to which MCLB may be accurately applied. Furthermore, our solution is computationally efficient and theoretically straightforward, since it requires no external force distribution, such as has been used previously [3,4] to generate appropriate continuum interfacial physics; the interface method advanced here more closely approximates a true boundary condition, which can be treated with the simplest form of Chapman-Enskog analysis [5]. In Sec. II selective background information is presented, upon which we shall build in Sec. III, where we introduce two closely related discontinuous hydrodynamic interface methods. In Sec. IV we will present results which we proceed to discuss and in Sec. V we present our conclusions and plans for further work. For the sake of simplicity, and to maintain parity with previous, related work [3,6], we consider the LB variant designated D2Q9 LBGK [7], which has lattice Q = 9 links denoted ci , i = 0–8, each with an associated weight denoted ti [5]. However, the work presented here will generalize very transparently to other multiple relaxation time LB variants, since no force distribution is used. We shall use common notation throughout. See Fig. 1 of Ref. [4] for the lattice link numbering convention in use. For simplicity, units of the lattice are used, so that x = t = 1 throughout. II. MCLB USING LISHCHUK’S METHOD
In this work, we rely upon certain aspects of our previous methods, i.e., that advanced by Lishchuk et al. in Ref. [3], combined with a distributed interface method after d’Ortona et al. [4,8]. We shall refer to this combination as “Lishchuk’s
063305-1
©2013 American Physical Society
T. J. SPENCER AND I. HALLIDAY
PHYSICAL REVIEW E 88, 063305 (2013)
method” [3,4] for brevity. In this section we review those of its features which are directly relevant and introduce related extensions required by the discontinuous hydrodynamic interface method. In the wake of the pioneering work of Gunstensen and Rothman [9], several two-component MCLB methods have arisen, each distinguished by the physical regime of its fluid-fluid interface. At the mesoscale, where the kinematics of phase separation must be considered, free energy methods [10,11] and the thermodynamically consistent extensions of Wagner et al. [12–14] (based, as they are, on the Cahn-Hilliard theory) are appropriate; for workers with a background in molecular simulation, the Shan-Chen method [15] is a natural choice. At the continuum scale, addressing completely immiscible fluids, the MCLB interface of Lishchuk et al. [3,4] provides a robust MCLB with assigned surface tension, appropriate dynamics, and, when used with an appropriate segregation [4], continuum interfacial kinematics. It is important to note that all the methods cited above contain distributed fluid-fluid interfaces, which extend over about >5 lattice sites (sites having −0.9 ρ N 0.9), with a variable body force applied on all of them, to produce a Laplace pressure step. Lishchuk’s method [3,4] imposes Laplace law and “no-traction” effects by impressing a curvature-dependent, immersed-boundary external force density, primarily in regions where fluid component phase parameter ρ N varies. The calculation of this force and the necessary correction to the model’s velocity [16] complicates the algorithm and is expensive. Let two fluid distributions which occupy lattice link i at position r and time t be described by distribution functions, Ri (r,t) and Bi (r,t), with the overall distribution function
3 (a)
ρ (lattice units)
2.5 2 1.5 1 0.5 0 -5
-4
-3
-2
-1 0 1 x (lattice units)
2
3
4
5
3 (b)
ρ (lattice units)
2.5 2 1.5 1 0.5 0 -5
-4
-3
-2
-1 0 1 x (lattice units)
2
3
4
5
3 (c)
ρ (lattice units)
2.5 2 1.5 1
fi (r,t) = Ri (r,t) + Bi (r,t) .
0.5
The nodal density of the red and blue fluid is R (r,t) = Ri (r,t) , B (r,t) = Bi (r,t) ,
0 -5
-4
-3
-2
-1 0 1 x (lattice units)
2
3
4
5
3
i
(d)
(2)
i
and the local phase field is [4]
2.5
ρ (lattice units)
(1)
ρ N (r,t) =
2 1.5
0.5 0 -4
-3
-2
-1 0 1 x (lattice units)
2
3
4
(3)
This definition is used without modification in the present work, to define an extended, smoothly varying phase field. Closed contour ρ N = 0 will define the discontinuous hydrodynamic interface, which, having no transverse extension, is always “sublattice.” Throughout the interfacial region, where |ρ N | < 1, a local interface normal is easily derived from derivative of ρ N :
1
-5
R (r,t) − B (r,t) . R (r,t) + B (r,t)
5
FIG. 1. The development of numerical instability as the interfacial pressure, or density step, , increases. Rest interface density profiles for red (triangles) and blue (squares) fluids obtained with various values of pressure step . The lines are included only to guide the eye. The interface is centered upon the sublattice position x0 = 0.5 lattice units, on a D2Q9 lattice, with density ρ0 = 2.0, interface segregation parameter β = 0.67, and LBGK collision parameter τ = 1. The fluids are at rest. Numerical differentiability degrades as increases. It is important to note that the stable system with = 0.4 depicted in [Fig. 1(b)] corresponds to a density step of 0%, which corresponds to a large Laplace law pressure step.
nˆ = −
∇ρ N . |∇ρ N |
(4)
The extended interfacial curvature field is obtained from the surface gradient of nˆ = (nˆ x ,nˆ y ), which, in two dimensions, is given by [3] ∂ nˆ y ∂ nˆ y ∂ nˆ x ∂ nˆ x − nˆ 2y + − nˆ 2x . (5) K(r) ≡ nˆ x nˆ y ∂x ∂y ∂x ∂y This field will be used to determine the curvature on the hydrodynamic interface ρ N = 0.
063305-2
MULTICOMPONENT LATTICE BOLTZMANN EQUATION . . .
PHYSICAL REVIEW E 88, 063305 (2013)
All the derivatives in Eqs. (4) and (5) may be computed to good accuracy by using the stencil ∂φ = 3 i ti φ(r + ci )ciα + o(x 3 ), (6) ∂xα r
(a)
cs2 ρ (lattice units)
in which φ(r) denotes ρ N , nˆ x , or nˆ y . The above result applies to the D2Q9 lattice. Shan has pointed out that the above stencil contains anisotropic lattice tensors in its error term [17]. Hence, normally directed forces calculated from it may contain tangential components which drive an interfacial microcurrent; we return to this issue in Sec. IV. The largest such anisotropic contribution in the error term of Eq. (6) may be eliminated, resulting in the following alternative stencil (again for the D2Q9 lattice): ∂φ = 3.2 i ti φ(r + ci )ciα ∂xα r (7)
u0 =
u(r )ρ N (r) − u(r)ρ N (r ) . ρ N (r) − ρ N (r )
ρ -1 -1.5
(c)
K0 (lattice units )
-1
|u| x10 (lattice units)
6
Figure 2(c) shows the value of curvature K0 interpolated using Eq. (10) on the contour ρ N = 0, for a circular drop of initial radius 35 lattice units. In Lishchuk’s method [3,4], continuum regime interfacial effects emerge by applying an external force density F = 1 ∇ρ N α0 K, where 12 ∇ρ N is the appropriate weight [3] and α0 2 is the interfacial tension. Impressing this force, one recovers a Laplace law pressure step [2] and the no-traction condition [3]. Inserting external forces, however, complicates LB’s underlying theory, its implementation, and reduces execution
0.03 0.025 0.02 0.015 0.01 0.005 0 -0.005 -0.01
(d)
(10)
(11)
0 -0.5
which, using Eq. (8) gives
A similar interpolation may be used to obtain the fluid velocity on the ρ N = 0 contour:
60
0.5
The local curvature K0 of the ρ N = 0 boundary, or fluid-fluid interface, at its intersection with ci may then be obtained from the interpolation: K(r ) − K(r) K0 = K(r) + d , (9) |ci | K(r )ρ N (r) − K(r)ρ N (r ) . ρ N (r) − ρ N (r )
40
1
We note that the phase field generates no hydrodynamic behavior whatsoever until some interaction between it and the LB distribution function fi is defined. It is for this essential reason that a continuous phase field may be used to produce a discontinuous hydrodynamic interface. In the present work, it will be necessary to evaluate the curvature at points on the contour ρ N = 0. This is achieved with sufficient accuracy, using a low order interpolation as follows. Identify all pairs of lattice sites r [with ρ N (r) > 0] and r [with ρ N (r ) < 0] using the criterion ρ N (r)ρ N (r ) < 0. Every such pair defines a lattice link ci which intersects the interface contour ρ N = 0 at distance d from the node at r, where |ci | d = N . (8) ρ N (r) ρ (r) − ρ N (r )
K0 =
20
140 120 100 80 60 40 20 80 100 120 140 0
(b)
N
− 0.1 i ti φ(r + 2ci )ciα + o(x ). 3
0.6668 0.66678 0.66676 0.66674 0.66672 0.6667 0.66668 0.66666 0.66664 0.66662 0.6666 0
3 2.5 2 1.5 1 0.5 0 -0.5
FIG. 2. For all these data the two-dimensional (2D) drop has initial radius of 35 lattice units. It was simulated on a lattice of 125 × 125 lattice units, without viscosity contrast, using LBGK collision parameter τ = 1.0, separation parameter β = 0.67, and an interfacial tension parameter α = 5.0 × 10−3 lattice units. All parts of this figure show nodal values, hence the smooth ρ N = 0 contour cannot be distinguished. The data represent the nodal value of steady state (a) discontinuous pressure field, (b) corresponding continuous phase field, ρ N (c) interpolated curvature value displayed at the lattice node closest to the ρ N = 0 contour, (d) interfacial microcurrent [note the scale factor in (d)]. The data of this figure are typical of both model variants described in Secs. III C1 and III C2.
speed, for to describe fluid subject to an external force density requires both the computation of that force and the insertion
063305-3
T. J. SPENCER AND I. HALLIDAY
PHYSICAL REVIEW E 88, 063305 (2013)
of a “source term” Fi into the evolution equation [16,18]: 1 † fi (r,t) = fi(0) (ρ,u) + 1 − fi(1) (∇ρ, . . . ∂x uy . . . ,F) τ (12) + Fi (τ,F,u) , 1 )( cic−u 2 2τ s
+ · F. Here we have conwith Fi = ti (1 − sidered the LBGK variant [7]; the dagger superscript indicates a postcollision, prepropagate quantity, fi(0) (ρ,u) denotes the equilibrium distribution function [5,7]:
uα uβ ciα ciβ − cs2 δαβ uα ciα (0) fi (ρ,u) = ti ρ 1 + 2 + , (13) cs 2cs4 (ci .u)(ci ) ) cs4
and density and velocity are defined as follows: 1 ρ≡ fi , ρu ≡ fi ci + F, 2 i i
τ together determine the viscous stress tensor, 1 i fi(1) (r)ciα ciβ , (r) = − 1 − σαβ 2τ
(20)
and that its weakly compressible hydrodynamics has a NavierStokes equation [20], D ∂ 2 ∂ ρuα = − cs ρ + σ . Dt ∂xα ∂xβ αβ
(21)
Note that lattice units have been used in Eqs. (19)–(21). III. A DISCONTINUOUS HYDRODYNAMIC INTERFACE FOR MCLB
(17)
Laplace law compliance and stress continuity boundary conditions require that two separated fluids’ hydrostatic pressure, p = cs2 ρ, evaluated on either side of a short section of interface, should differ by an amount α0 K0 , where K0 is the local interface curvature and α0 is the interfacial tension while tangential and normal viscous stress contractions should equate [2]. The kinematic condition of mutual impenetrability requires that the two fluids should move with the same velocity. Let us consider first the interface between two immiscible fluids, still designated red and blue, having identical shear viscosity, separated by an effectively discontinuous interface. Later we will consider viscosity contrast, i.e., a change in the value of collision parameter τ (which actually introduces no further complications). How are such dynamic and kinematic boundary conditions to be incorporated into a simple, efficient, and accurate MCLB simulation? It is appropriate to recall that there is, in fact, only a single “background” lattice fluid in all MCLB methods and that only physical observables like pressure and density need to change discontinuously.
in which parameter β (β 0.69, note) regulates phase field width. Furthermore, the red color density (say) of a locally flat interface with normal nˆ embedded between fluids in uniform motion u is described by [19] β ρ0 ˆ − ut) ; 1 − tanh R(x,t) = n.(r (18) 2 x
To comply with continuum physics, a LB distribution function which streams, or propagates, across a sublattice boundary ρ N = 0 must see only its equilibrium part, fi(0) (ρ,u), change by an amount which is consistent with the lattice fluid density change associated with the Laplace pressure step:
(14)
all other symbols having their usual meaning [4]. Streaming may be written †
fi (r + ci ,t + 1) = fi (r,t).
(15)
In the present work, we will obtain the same physics without recourse to an external force density F. Note that for F = 0, Fi = 0 and Eq. (12) simplifies considerably. Physical interfacial kinematics and isotropy emerge from carefully chosen postcollision color segregation rules. We have previously shown [4] that correct motion of the phase field boundary accrues from the use of d’Ortona’s segregation [8], R(r,t)B(r,t) R(r,t) † † ˆ Ri (r,t) = (16) f (r,t) + βti ci .n, ρ(r,t) i ρ(r,t) †
†
†
Bi (r,t) = fi (r,t) − Ri (r,t),
here x is lattice spacing. Recall that in this work x = 1. This result shows that the time development of a flat interface of the type used in this work corresponds to advection at the same speed as the underlying fluid, which is consistent with the kinematic condition of mutual impenetrability, and demonstrates the Galilean invariance of the color algorithm. Note that while the variation of R and B may be influenced by β, in Sec. III F the variation of ρ = R + B will be shown to be independent of β. It is relevant to recall that the LBGK method’s kinematic viscosity is determined by the collision parameter: 2τ − 1 η , (19) ν= = ρ 6 that the principal part of its distribution function, fi = Ri + Bi , is the equilibrium contribution fi(0) (ρ,u) which, note, is linear in the density ρ, as is clear in Eq. (13), that its nonequilibrium contribution to fi , fi(1) and collision parameter
A. Kinematic and dynamic compliance
δP = α0 K0 .
(22)
It will often be convenient to think of the corresponding density step: ≡ R0 − B0 ,
(23)
where R0 and B0 are the density of the red and blue fluids at rest, at a large distance from any interface. Clearly, =
α0 K0 , cs2
(24)
by the Laplace law. Any distribution function which streams across a sublattice boundary should see its equilibrium part change by an amount fi(0) (,u), where = α0 K0 /cs2 , while in the absence of viscosity and density contrasts, its nonequilibrium part should not change (so that viscous stresses and, hence, their normal and tangential contractions, will be continuous).
063305-4
MULTICOMPONENT LATTICE BOLTZMANN EQUATION . . . B. Background
cs2 ρ
The aim of adjusting the value of in accord with the Laplace law while maintaining f (1) continuous across our MCLB interface is facilitated by previous work, which we now summarize. To simulate pressure driven internal flows, Kim and Pitsch devised a periodic boundary closure which produces a spatially and temporally fixed pressure step of constant magnitude cs2 at, say, x = x0 , with the velocity and fluid stresses continuous across that boundary [21]. This closure may be thought of as the exchange of a constant “quantum” of density ti between appropriately chosen LB distribution function component pairs, defined as follows. Suppose lattice link ci (i > 0) starts from the lattice node at position r = (x,y) and crosses the x = x0 boundary and that ρ(r) = ρ(r + ci ) + , i.e., ρ(r) > ρ(r + ci ). LB distribution functions fi (r,t) and fj (r + ci ,t), with j = i ⊕ (Q − 1)/2, form a pair which counterpropagate along link ci , between times t and t + 1. Here, ⊕ denotes addition, modulo Q (Q = 9 here, of course). The distribution function pair fj (r + ci ,t) and fi (r,t) is characterized by the two values of r and i [j = i ⊕ (Q − 1)/2 is determined by the choice of i]. The “opposed” distribution function pair fi (r,t) and fj (r + ci ,t), chosen in this way, has equilibrium components determined solely by the density (pressure) and velocity of the fluid, which we denote fi(0) (ρ(r),u(r)), fj(0) (ρ(r + ci ),u(r + ci )),
(25)
still j = i ⊕ (Q − 1)/2, note. As we have indicated, to ensure that the required pressure difference cs2 exists between the fluid occupying lattice nodes r and r + ci it is necessary to enforce the pressure step condition: ρ(r) − ρ(r + ci ) = .
(26)
This can be achieved, Kim and Pitsch showed, by that notional exchange, during the LB propagate step, of a fixed quantum of density ti , implicit in the following replacements, imposed at every time step: fi (r,t) → fi (r,t) −
fi(0) (,u(r)),
fj (r + ci ,t) → fj (r + ci ,t) + fj(0) (,u(r + ci )).
(27) (28)
Note that the fi(1) component of the distribution is left intact. However, the total mass of the system is invariant only if fi(0) (,u(r)) = fj(0) (,u(r + ci )). We return to this point in Sec. IV. The above replacements are readily computed, using Eq. (13), for a chosen value of . To produce interfacial physics appropriate to a discontinuous hydrodynamic boundary between two immiscible fluids, designated red and blue, one can apply a version of the steps defined in Eq. (27), above. We present below two related algorithms, each of which does essentially this. C. Algorithm
Let us consider the interface of a two component MCLB model. The latter is henceforth defined as the closed, sublattice contour ρ N = 0 of the continuous phase field and is denoted C. Let us suppose that the red fluid is enclosed by the blue fluid and therefore subject to a Laplace pressure, that C has
PHYSICAL REVIEW E 88, 063305 (2013)
known local curvature K0 , that lattice node r lies inside C [so ρ N (r) > 0], and that lattice node r ≡ r + ci lies outside C [so ρ N (r + ci ) < 0]. The density of the red fluid differs from that of the blue fluid by an amount = α0 k0 /cs2 , determined by the Laplace pressure, and (for the moment) no viscosity contrast exists. 1. Mass conserving discontinuous hydrodynamic boundary
We use the following extension to the MCLB evolution rule set, at every time step, on all the distribution function component pairs “connected” by the link ci (which intersects C, recall): Ri (r,t) → Ri (r,t) − 12 fi(0) (,0) ,
(29)
Bi (r,t) → Bi (r,t) − 12 fi(0) (,0) ,
(30)
Rj (r ,t) → Rj (r ,t) + 12 fj(0) (,0) ,
(31)
Bj (r ,t) → Bj (r ,t) + 12 fj(0) (,0) .
(32)
Here, = α0 K0 /cs2 . K0 , the curvature on contour C, is interpolated using Eq. (10) from the curvature field defined by Eqs. (4) and (5). In Eqs. (4) and (5) all derivatives of ρ N and nˆ are evaluated using the stencil in Eq. (6). To a first approximation, it is appropriate that additional mass is allocated equally between red and blue components in (29)–(32) because, on C, the value of ρ N ≈ 0, hence the fluid there will be composed of red and blue components equally. The sum Ri (r,t) + Rj (r ,t) is invariant under replacements (29)–(32) which therefore conserve colored mass. There is, however, local momentum introduced between the pair of sites at r and r = r + ci . Replacements (29)–(32) introduce a net momentum: Pij (r,r ) = fi(0) (,0,t)ci + fj(0) (−,0,t)cj .
(33)
Now ci = −cj , weights ti = tj , and so, from Eq. (13): Pij (r,r ) = 2ti ci , where = a drop is therefore
α0 K0 /cs2 ,
(34)
recall. The total momentum imparted to
P = 2
ti ci .
(35)
r
Whatever its shape, a closed contour must cross any lattice line an even number of times. Therefore, the contributions to the local momentum r ti ci cancel in pairs and the total momentum change P(r,r ). Hence, our first, most computationally economical, colored mass conserving discontinuous hydrodynamic boundary method also conserves global (though not local) momentum. 2. Low microcurrent discontinuous hydrodynamic boundary
The more costly method of this section approximately conserves color mass and is based upon a more considered approach to allocating the transferred mass in the interfacial region. The interfacial microcurrent is a small unphysical velocity which arises in the region of the interface. This spurious signal is substantially reduced when the following, modified
063305-5
T. J. SPENCER AND I. HALLIDAY
PHYSICAL REVIEW E 88, 063305 (2013)
extension to the MCLB evolution rule set is used at every time step, to modify the values of the “opposed,” counterpropagating distribution function component pairs, connected by the link ci : Ri (r,t) → Ri (r,t) −
g(i,r)fi(0)
(,u0 ) ,
(36)
Bi (r,t) → Bi (r,t) − g (i,r)fi(0) (,u0 ) ,
(37)
Rj (r ,t) → Rj (r ,t) + g(j,r )fj(0) (,u0 ) ,
(38)
Bj (r ,t) → Bj (r ,t) + g (j,r )fj(0) (,u0 ) .
(39)
The factors g(i,r) and g (i,r) which appear above will be discussed shortly. In these replacements, = α0 K0 /cs2 . K0 , the curvature on contour C, is interpolated using Eq. (10). u0 denotes the velocity interpolated onto the ρ N contour C, again using Eq. (10). It is important to note that the underlying curvature field is now calculated from Eqs. (4) and (5) but now the derivatives of ρ N and nˆ are evaluated using the modified stencil in Eq. (7). The factors g(r) and g (r) appearing in replacements (36)–(39) are defined Ri (r,t) , g(i,r) ≡ Ri (r,t) + Bi (r,t)
Bi (r,t) g (i,r) ≡ . Ri (r,t) + Bi (r,t) (40)
They are designed to conserve the value of the phase field ρ N in the adjusted link populations, as follows. Consider the node at r . It sends information in the form of Rj (r ,t) and Bj (r ,t) along link cj , from node r , where the value of the phase field is ρ N (r ). However, Rj (r ,t) and Bj (r ,t) are adjusted by replacements (36)–(39). Temporarily drop the variables t and r . In order to preserve the perceived value, at node r, of the phase field ρ N (r ), the changes in the population, δRj , δBj , should satisfy the condition ρjN ≡
Rj − B j Rj + δRj − (Bj + δBj ) = . Rj + B j Rj + δRj + Bj + δBj
(41)
Taking a binomial expansion, we obtain δRi − δBj δRj + δBj = + o δRj2 , Rj − B j Rj + B j
(42)
which may be solved δBj =
Bj δRj , Rj
(43)
and since δRj + δBj = fj(0) (,u0 ) it is straightforward to obtain Rj fj(0) (,u0 ), δRj = Rj + B j Bj δBj = fj(0) (,u0 ), (44) Rj + B j which explains the definition of the factor g(j,r ) appearing in replacements (36)–(39). Clearly, colored mass is now not conserved, and small fluctuations in mass arise. We return to this issue in Sec. IV. The net momentum introduced between the pair of sites at r and r = r + ci by the replacements (36)–(39) is given by
the expression Pij (r) = fi(0) (,u0 )ci + fj(0) (−,u0 )cj .
(45)
Again, ci = −cj and the link weights equate, ti = tj , so, taking a linear approximation from Eq. (13), (u0 − u0 ) · ci ci + o u20 , Pij (r,r ) = 2ti ci + ti 2 cs (46) where = α0 K0 /cs2 , recall. The total momentum imparted to a drop is therefore (47) ti ci + o u20 . P = 2 r
Whatever its shape, a closed contour must cross any lattice line an even number of times. Therefore, the contributions to r ti ci again cancel in pairs and the principal contribution to the total momentum, above, is zero to o(u20 ). D. Viscosity contrasts
Thus far, we have considered immiscible fluids with identical viscosity, separated by an interface located on a sublattice contour C, i.e., on the closed curve ρ N = 0. Let us now allow a change in fluid viscosity across contour C. Thus, the red (blue) fluid occupying lattice node r (r = r + ci ) has a value of kinematic viscosity determined by τr (τr = τb ). Physical fluid with a shear viscosity η, which varies, possibly in response to local heating, is described by Eq. (21) [2], which shows that viscous stress must be continuous (i.e., differentiable) in the fluid, given that velocity and density are also continuous. Likewise, in LBGK dynamics, viscous stress must be continuous in the presence of inhomogeneities in τ if the velocity and density are continuous. In the present context, velocity across the region of the discontinuous interface of Sec. III A is continuous by design [21] and since we only effect a step change in pressure (density) at the sublattice C, long wavelength variations in density are also continuous. It follows that the viscous stress is continuous throughout the interfacial region. Put another way, a rapid change in shear viscosity η will produce compensating changes in the shear rate sαβ but the viscous stress σαβ will remain continuous. We return to this point in Sec. IV shortly. Given the continuity of the viscous stress [20], 1 i fi(1) (r)ciα ciβ , (48) σαβ (r) = − 1 − 2τ in the presence of a sudden change in viscosity, i.e., LBGK collision parameter τ , interfacial viscous stress contractions ˆ σαβ tα , with tˆ the assumed smoothly varying unit tangent vector, must also be continuous. Clearly, this is a necessary condition for the emergence of the appropriate hydrodynamic boundary (red) (blue) conditions, such as σαβ tˆβ − σαβ tˆβ = 0. E. Density contrasts
The additional algorithmic extension we have used previously to amplify density contrasts [22] can be applied for the present model. Of course, compounding algorithmic complexity must have a range of adverse consequences.
063305-6
MULTICOMPONENT LATTICE BOLTZMANN EQUATION . . .
However, the addition of the stress insulating forces postulated previously is perfectly feasible. For the sake of brevity we reserve this investigation for a future publication. F. Analysis
Consider a flat MCLB fluid-fluid interface orientated perpendicular to the x axis, maintained by the segregation rule in Eq. (16), embedded in a uniform fluid of constant density ρ0 . By Taylor expanding the discrete kinetic equation of the MCLB and considering the steady state, we have shown elsewhere [4] that such an interface when centered upon x = x0 [i.e., with ρ N (x) < 0,x < x0 ], will be described by the ordinary differential equation, β d 2 d2 β d R+2 R = 0, R−2 dx 2 h dx hρ0 dx
(49)
in which h is the lattice spacing. Note that this equation has a first integral. Throughout this work we have taken h = 1, with which simplification the solution is [4,19] 1 + tanh (β(x − x0 )) R(x) = ρ0 , 2 1 − tanh (β(x − x0 )) . (50) B(x) = ρ0 2 Here, it has been assumed that the region x < x0 is initially occupied by blue fluid. Simulation data support this solution. We seek now to determine the corresponding red fluid density variation in the presence of a density discontinuity: ρ(x) = ρ0 + δ[1 − H (x − x0 )],
(51)
where H (x) is the Heaviside step function. Noting that the point x0 will never correspond to a point lattice location, one can follow an identical approach to that used to derive Eq. (49). An additional assumption is that for x < x0 (x > x0 ) ρ(x) = ρ0 + (ρ0 ) is used to encapsulate the boundary condition implicit in the rule set of Eqs. (29)–(32), straightforwardly to show that the introduction of a density discontinuity at x = x0 modifies only one term in the ordinary differential equation (49): β d 2 d2 β d R+2 R = 0. R−2 2 dx h dx hρ(x) dx
that the phase field does not contain any discontinuity, ρ N (x) ≡
Note that this solution meets those density boundary conditions implicit in our replacement expressions (29)–(32), namely Eq. (51), (55)
R(x) − B(x) = tanh (β(x − x0 )), R(x) + B(x)
(56)
and that the latter will therefore remain differentiable, provided certain straightforward conditions are met. In practice, to obtain curvatures K, the phase-field profile ρ N must be numerically differentiable, without taking elaborate precautions (which would undermine the simplicity of the method). Implicit in the numerical differentiability of ρ N is the differentiability of R(x) and B(x). For R(x) to be differentiable, the step change in R(x), which occurs at x = x0 and is of size /2, must be small compared with the change in the continuous contribution tanh (β(x − x0 )) in the region of x0 . This leads directly to the condition d ρ0 − tanh (β(x − x0 )) (57) . dx 2 2 x=x0 Performing the differentiation and noting that cosh2 (β(x − x0 )) = 1 when x = x0 we obtain β(ρ0 − ) which may be re-arranged to provide the following stability criterion:
βρ0 . 1+β
(58)
Figure 1 shows the steady-state color density profiles associated with a flat interface y = const, centered upon the sublattice position x0 = 0.5 lattice units, for a D2Q9 lattice, with density ρ0 = 2.0, a range of values of the density step , segregation parameter β = 0.67, and LBGK collision parameter τ = 1. The fluids are at rest. From inequality (58), the approximate value at which the phase field ρ N ceases to become numerically differentiable is 0.8, which agrees with the trend in the figure. The existence of this instability must not be overemphasized; a density step = 0.4 [Fig. 1(b)] corresponds to a 20% density step and hence to a large Laplace law pressure step and interfacial tension. Moreover, Eq. (58) shows that it may be controlled by an appropriate choice of ρ0 and β. The latter is not now subject to restrictions placed upon it using immersed boundary techniques, as we shall now see.
(52)
For x = x0 , ρ(x) is constant and it is possible to write down a solution to Eq. (52) from Eq. (49), by replacing ρ0 with ρ(x). Using Eq. (51), we obtain ρ0 + R(x) = − H (x − x0 ) [1 + tanh (β(x − x0 ))] , 2 2 (53) ρ0 + B(x) = − H (x − x0 ) [1 − tanh (β(x − x0 ))] . 2 2 (54)
ρ(x) ≡ R(x) + B(x) = ρ0 + [1 − H (x)],
PHYSICAL REVIEW E 88, 063305 (2013)
IV. RESULTS AND DISCUSSION
The microcurrent, or spurious velocity, associated with MCLB interfaces is a small but unphysical circulation in the near-interfacial region occurring even when the separated fluids are otherwise at rest. See Fig. 2(d). We shall denote the maximum microcurrent velocity U . The vorticity created at an interface will diffuse, at a rate controlled by the fluid viscosity, into the bulk fluids where it may accumulate, unless balanced by dissipation. Lee et al. have successfully reduced microcurrents in free energy MCLB methods [23] and Shan [17] has argued that they are, in the Shan-Chen model [15], caused by lattice anisotropy and that they are therefore expensive to eliminate. It is not our aim here to carefully consider the origin of this essentially numerical artifact. We adopt a more practical approach. Strong flow obscures the microcurrent. Conversely, microcurrent activity obscures physical flow when the latter is
063305-7
T. J. SPENCER AND I. HALLIDAY
PHYSICAL REVIEW E 88, 063305 (2013)
Min(Ca) ≡
ηU (2τ − 1)U = . α0 6α0
(59)
By measuring a steady state U , then evaluating Min(Ca) defined above, it is possible to estimate the value of Ca at which, for a given level of resolution, the physical flow will cease to be distinguishable at all points in a simulation. A variety of other factors including lattice resolution and the value of the collision parameter τ also affect U . We fix lattice resolution for all the results presented in this section. We consider throughout a red drop of constant radius R = 35 lattice units, at rest at the center of a lattice of size 125 × 125 lattice units, subject to periodic boundary conditions. The systems simulated are therefore infinite in extent and contain periodic images of the red drop. Accordingly, the velocity induced by the microcurrent vorticity can increase for a considerable time. Throughout this section data obtained after 1.0 × 106 lattice time steps is designated steady state. The initial density of both red and blue fluids was set to a value 2.0. The hydrodynamically discontinuous interfaces developed here may be assessed by comparing their Min(Ca) with the corresponding quantity obtained from otherwise identical simulations using Lishchuk’s method. Let us consider first the case of zero viscosity contrast, setting τr = τb = 1. For the steady state of a drop, using segregation parameter β = 0.67 and a large range of interfacial tensions α0 , Fig. 3 shows the variation in Min(Ca), for the mass conserving discontinuous interface method of Sec. III C1 compared with Lishchuk’s method [3,4](continuous line). Figure 4 shows the equivalent variation in Min(Ca), for the low microcurrent discontinuous interface method of Sec. III C2, again compared with Lishchuk’s method [3,4] (continuous line). See the figure captions for a summary of the simulation parametrizations. From Fig. 3, it is clear that for α0 0.02 Lishchuk’s method [3,4] has the practical advantage, However, for smaller values of α0 a discontinuous interface may be used without apparent penalty. Indeed, it may be used to advantage. In these data, as α0 approaches zero, Min(Ca) approaches a constant, which indicates the existence of a regime where U ∝ α0 . More unexpected is the existence of an optimum value of α0 for the mass conserving hydrodynamically discontinuous
-2
10
min (Ca)
10-3
10-4
-5
10
-6
10
-5
10
-4
10
-3
-2
10 10 α0 (lattice units)
-1
10
0
10
FIG. 3. For all the data in this figure, the drop was assigned an initial radius of 35 lattice units, it was simulated on a lattice of 125 × 125 lattice units, without viscosity or density contrast, using LBGK collision parameter τ = 1.0 for both fluids, separation parameter β = 0.67, and a large range of interfacial tensions α0 . The solid squares show the variation of Min(Ca) for the mass conserving discontinuous hydrodynamic interface method, in Sec. III C1. For reference, the continuous line shows the value of Min(Ca) for Lishchuk’s method [3,4], which is an immersed boundary technique.
interface method. The data in Fig. 4 shows that the effect of the microcurrent is indeed reduced and that our low microcurrent discontinuous hydrodynamic interface method actually performs better than Lishchuk’s method over the range of interfacial tension parameter investigated. The value of segregation parameter β may be decreased, to make the phase field boundary more diffuse and, hence, to reduce the microcurrent. This exacerbates the problem of unphysical scales in the case of Lishchuk’s method. However, while reducing the value of β makes the phase field vary over a greater distance, the ρ N = 0 contour remains sharp and the -2
10
10-3 min (Ca)
weak. The length scale associated with the interface has been eliminated from the hydrodynamic signals’ spatial variation by the methods related above. Therefore, where flow is strong, any penalty in the form of an increase in U will be unimportant and the discontinuous hydrodynamic interface is clearly the preferred method. However, in weak flow the practical utility of the discontinuous interface method may be restricted by an increased in U . Hence, our current methods can be succinctly evaluated by assessing their maximum microcurrent velocity compared with an equivalent implementation of Lishchuk’s method, in the context of weak flow. While the relationship may not be linear, as the interfacial Laplace pressure step increases (in response to an increase in interfacial tension α0 ) so does the spurious velocity U . To take account of this effect, we characterize a given MCLB method by the dimensionless minimum capillary number accessible to simulation using that MCLB interface. Our minimum capillary number is defined as follows:
-4
10
10-5 -6
10
-7
10
-5
10
-4
10
-3
-2
10 10 α0 (lattice units)
-1
10
0
10
FIG. 4. For all the data in this figure, the drop was assigned an initial radius of 35 lattice units, it was simulated on a lattice of 125 × 125 lattice units, without viscosity or density contrast, using LBGK collision parameter τ = 1.0 for both fluids, separation parameter β = 0.67, and a large range of interfacial tensions α0 . The solid squares show the variation of Min(Ca) for the low microcurrent discontinuous hydrodynamic interface method of Sec. III C2. For reference, the continuous line shows the value of Min(Ca) for Lishchuk’s method [3,4].
063305-8
MULTICOMPONENT LATTICE BOLTZMANN EQUATION . . .
PHYSICAL REVIEW E 88, 063305 (2013)
-2
10-2
10
-3
10-3
10-4
10
-4
min (Ca)
min (Ca)
10
-5
10
10-5
10-6
10-6
10-7
10
-7
-8
10
10-5
10-4
10-3 10-2 α0 (lattice units)
10-1
10-8 -5 10
100
FIG. 5. For all the data in this figure, the drop was assigned an initial radius of 35 lattice units, it was simulated on a lattice of 125 × 125 lattice units, without viscosity or density contrast, using LBGK collision parameter τ = 1.0 for both fluids, separation parameter β = 0.4, and a large range of interfacial tensions α0 . The solid squares show the variation of Min(Ca) for the low microcurrent discontinuous hydrodynamic interface method of Sec. III C2. For reference, the continuous line shows the value of Min(Ca) for Lishchuk’s method [3,4] with β = 0.4. Note that using width of the interface in the latter is now considerable; these data should be considered only as reference data.
discontinuous hydrodynamic interface method may be used with β = 0.4 without penalty. Figure 5 compares Lishchuk’s method with our low microcurrent method, which again outperforms the latter, over the range of interfacial tension parameter investigated. See the figure caption for a summary of the simulation parametrization. It is very important to note that Figs. 3–5 show that the low microcurrent version of our algorithm with β = 0.4 considerably outperforms Lishchuk’s method with β = 0.69. Figure 6 shows data similar to that in Fig. 4, obtained with our low microcurrent method. The sole difference is the value of LBGK collision parameter τ = 5.56, corresponding to a tenfold increase of both fluids’ kinematic viscosities. Again, the solid squares show the variation of Min(Ca) for the low microcurrent discontinuous interface and, as in Fig. 3, the broken line shows the value of Min(Ca) for Lishchuk’s method [3,4]. Color masses are not conserved for the systems of Figs. 4 and 5. It is important to quantify that mass fluctuation. Suppose a red drop of initial red mass R0 has, after T time steps, a red mass RT . Let L denote the length of its interface, to which quantity the mass fluctuation must be proportional. The quantity , defined below, measures the fractional rate of change in red (drop) mass, expressed as a percentage of the original mass, R0 , per unit length of interface, RT − R0 100 . (60) ≡ T R0 L By fitting simulation data, over the range of interfacial tensions, α0 < 0.01, investigated the following empirical model for was obtained: = aα02 + bα0 + c,
(61)
-4
10
-3
-2
10 10 α0 (lattice units)
-1
10
0
10
FIG. 6. For these data, the drop was assigned an initial radius of 35 lattice units; it was simulated on a lattice of 125 × 125 lattice units, again without viscosity or density contrast, using an increased LBGK collision parameter τ = 5.56 (viscosities) for both fluids. The separation parameter β = 0.67 and a large range of interfacial tensions α have been used. The solid squares show the variation of Min(Ca) for the discontinuous hydrodynamic interface and, as in Fig. 3, the broken line shows the value of Min(Ca) for Lishchuk’s method [3,4]. The data for that method are little changed from those in Fig. 3. By contrast, those obtained from a discontinuous interface show an increase in microcurrent activity over the whole range of α0 , especially for larger values.
with coefficients a = 3.70 × 10−6 , b = 4.90 × 10−7 , and c = 5.12 × 10−9 , obtained by least squares fitting and characteristic of LB model and dimensionality. Clearly these values are very small indeed. For a small drop with radius 20 lattice units (and, proportionally, a long interface) and large interfacial tension α = 0.1 lattice units, a two-dimensional (2D) simulation using the discontinuous hydrodynamic interface method of Sec. III C2 would need to run for approximately 1 × 105 time steps before its mass would have changed by 1%. Let us next consider a viscosity contrast between the red and blue fluids. For simplicity, we consider a flat interface x = 0 and shear the fluids in the y direction, in the tangential direction. The red fluid occupies the region x 0 and has a kinematic viscosity νr = 10 (i.e., τr = 30.5). The blue fluid occupies the region x > 0 with νb = 1 (τb = 3.5). Hence, a viscosity contrast of 10 exists between the fluids. Both fluids are initially at rest. In Fig. 7, which concerns the near interfacial region and not the entire flow domain, note that the variation of the computed MCLB model’s viscous stress, evaluated from Eq. (20), is expressed as a fraction of its steady state value (solid circles). Solid diamonds correspond to t = 250 time steps, solid squares correspond to t = 500 time steps, and the solid triangles correspond to t = 1000 time steps. The viscous stress is apparently continuous at all times. Figure 8 shows an excerpt of the steady state velocity profile of the system. The steady state shear rates in the different fluids are in the expected ratio. Numerical expressions for gradients in ρ N are used in the ˆ present method to determine the interface normal vector n. Hence, the arguments of Shan [17] regarding the origin of the microcurrent do apply, albeit indirectly, for the present method appears to avoid use of explicit force. An effective
063305-9
T. J. SPENCER AND I. HALLIDAY
PHYSICAL REVIEW E 88, 063305 (2013)
two methods was observed to arise in two distinct contributions. First, the introduction of the factors g(i,r) and g (i,r), defined in Eq. (40), does not affect qualitative behavior but it does dampen the numerical instability shown in Fig. 1. That figure was obtained with our low microcurrent model; the equivalent data obtained for the mass conserving variant show the instability occurring for smaller values of δ. Second, the change in template used for the calculation of the derivatives of continuous phase field ρ N , from that in Eq. (6) to that defined in Eq. (7), produces a further, independent reduction in microcurrent activity. Note that more extensive use of this stencil was observed to produce only small benefits. We reserve further comment for our conclusions.
2.5
σxy(t)/σxy0
2 1.5 1 0.5 0 -6
-4
-2 0 2 y (lattice units)
4
6
FIG. 7. Variation of the viscous stress, evaluated from Eq. (20), expressed as a fraction of its steady value (solid circles) about a flat interface, at which the kinematic viscosity changes. For x 0, ν = 10 (τr = 30.5), x > 0, ν = 1 (τb = 3.5), so that a viscosity contrast of 10 exists. Diamonds correspond to t = 250 time steps, solid squares correspond to t = 500 time steps, and the solid triangles correspond to t = 1000 time steps. The viscous stress is apparently continuous at all times.
force, or momentum change [note Eq. (46), above], clearly accompanies any process that does not explicitly conserve momentum. This fact, and similarity between the level of microcurrent activity measured with the present method and that observed when using Lishcuk’s method [3,4], across a range of α0 , adds further support to the relevance of Shan’s theory to both models considered in this work. In summary, Figs. 3–5 demonstrate that the mass conserving discontinuous hydrodynamic variant, defined in Sec. III C1, to be precise in steps (29)–(32), outperforms our previous models, until large interfacial tensions create Laplace pressure, or density, steps which are associated with that numerical instability identified in Fig. 1; we return to this point shortly. The low microcurrent variant presented in Sec. III C2 and summarized in steps (36)–(39) performs still better and, furthermore, is stable for a larger range of interfacial tensions. The improvement between the performance of our 1 0.98 0.96 ux/uxmax
0.94 0.92 0.9 0.88 0.86 0.84 0.82 -6
-4
-2
0
2
4
6
x (lattice units)
FIG. 8. Steady-state shear rate corresponding to the system used to obtain the data for Fig. 7. The steady-state shear rates are in the expected ratio.
V. CONCLUSIONS AND FURTHER WORK
Flow calculations in the arrested coalescence regime of continuum hydrodynamics are facilitated by the multiple immiscible fluids lattice Boltzmann equation method introduced here. The method removes unphysical interfacial scales from the hydrodynamic signals and simplifies both the computational algorithm and theoretical framework. Our “discontinuous hydrodynamic interface” method still relies upon an auxiliary phase field which changes continuously which limits “direct” computational cost savings, as we discuss below. However, it still provides indirect cost savings in applications where, e.g., computational mesh resolution is limited and comparisons with boundary integral or analytical calculation are sought by removing the need to simulate at larger than necessary system sizes to achieve the narrow interface limit. Put another way, by using the discontinuous interface method advanced here, the simulator is always in the narrow interface, hydrodynamic limit. The discontinuous hydrodynamic interface method is simple and easy to implement; it needs only low order interpolation to be stable and accurate over a range of interfacial tensions and it can be generalized to three dimensions and to multiple relaxation time schemes easily, since it relies upon no external force distribution whatsoever. We have assessed it against our existing approach (designated here as Lishchuk’s method [3,4]) which is also adapted to the continuum regime. When considering two fluids of identical density (but possibly different viscosity) separated by an interface which is subject to the Laplace law and no-traction conditions on interfacial stresses, our method introduces a small gain when assessed in terms of a simulation’s minimum accessible capillary number over a range of interfacial tension parameters when flow is weak. For very large interfacial tensions however, our methods produce microcurrent activities in excess of those produced in Lishchuk’s method, which is important when the simulated flow is challengingly weak. Matters are less ambiguous where flow is strong compared with the microcurrent activity; here the present method will be seen to have the clear advantage of physical accuracy over Lishchuk’s method. A natural question arises regarding further work to incorporate, in additional to a Laplace pressure step, a fluid density step at the interface. While such a density step is simple to directly incorporate, only small factors of density of order unity are possible, due to the numerical instability shown in Fig. 1 and considered in Sec. III F. To facilitate careful,
063305-10
MULTICOMPONENT LATTICE BOLTZMANN EQUATION . . .
PHYSICAL REVIEW E 88, 063305 (2013)
fully 3D comparisons with similar techniques [24], it will be necessary to address the matter of large density difference by adapting previous work [22] and to extend the model’s dimensionality. Assessed directly, the discontinuous interface method is only a little more computationally efficient than a comparable force-density based scheme and its principal advantages remain the removal of the unphysical scales from continuum multicomponent flow calculations and the simplification of the multicomponent lattice Boltzmann equation method algorithm. The introduction of a sharp physical interface is achieved but a diffuse auxiliary phase field and the associated processing cost remain to be found. This raises another question in respect of the extent of any direct saving in computational cost. Direct cost savings are realized because (i) a distributed immersed boundary force density is not, in the present method, impressed throughout the interfacial
region and (ii) no associated velocity field corrections are required [16] (our models are amenable to the simplest form of Chapman-Enskog analysis [5] used for lattice Boltzmann models [20]). Put another way, all direct computational costs associated with the processing of the phase field remain. Hence we must be cautious in making claims for the computational savings associated with the algorithmic simplifications presented here. Where high resolutions are essential, an extended phase field boundary region will be required and the computational savings are, in relative terms, small but they improve as resolution decreases, for one is always in the narrow interface limit. Stressing that, in general, direct computational cost savings must depend upon the extent of interface in a particular calculation, as well as other factors like geometry, we can tentatively quantify the computational cost savings associated with the results of Sec. IV at about 10% reduction in execution time.
[1] C. S. Peskin, J. Comput. Phys. 220, 20 (1977). [2] L. Landau and E. M. Lifshitz, Fluid Mechanics, 2nd ed. (Pergamon, New York, 1987). [3] S. V. Lishchuk, C. M. Care, and I. Halliday, Phys. Rev. E 67, 036701 (2003). [4] I. Halliday, A. P. Hollis, and C. M. Care, Phys. Rev. E 76, 026708 (2007). [5] S. Succi, The Lattice Boltzmann Equation for Fluid Mechanics and Beyond (Oxford University Press, Oxford, 2001). [6] T. J. Spencer, I. Halliday, and C. M. Care, Phys. Rev. E 82, 066701 (2010). [7] Y. H. Qian, D. d’Humi`eres, and P. Lallemand, Europhys. Lett. 17, 479 (1992). [8] U. D’Ortona, D. Salin, M. Cieplak, R. Rybka, and J. Banavar, Phys. Rev. E 51, 3718 (1995). [9] A. K. Gunstensen, D. H. Rothman, S. Zaleski, and G. Zanetti, Phys. Rev. A 43, 4320 (1991). [10] M. R. Swift, W. R. Osborn, and J. M. Yeomans, Phys. Rev. Lett. 75, 830 (1995). [11] M. R. Swift, E. Orlandini, W. R. Osborn, and J. M. Yeomans, Phys. Rev. E 54, 5041 (1996).
[12] A. J. Wagner, Phys. Rev. E 74, 056703 (2006). [13] Q. Li and A. J. Wagner, Phys. Rev. E 76, 036701 (2007), and references therein. [14] A. J. Wagner and C. M. Pooley, Phys. Rev. E 76, 045702(R) (2007), and references therein. [15] X. W. Shan and H. D. Chen, Phys. Rev. E 49, 2941 (1994). [16] Z. Guo, C. Zheng, and B. Shi, Phys. Rev. E 65, 046308 (2002). [17] X. Shan, Phys. Rev. E 73, 047701 (2006) [18] X. He, S. Chen, and G. D. Doolen, J. Comput. Phys. 146, 282 (1998). [19] I. Halliday, S. V. Lishchuk, T. J. Spencer, G. Pontrelli, and C. M. Care, Phys. Rev. E 87, 023307 (2013). [20] S. Hou, Q. Zou, S. Chen, G. D. Doolen, and A. C. Cogley, J. Comput. Phys. 118, 329 (1995). [21] S. H. Kim and H. Pitsch, Phys. Fluids 19, 108101 (2007). [22] S. V. Lishchuk, I. Halliday, and C. M. Care, Phys. Rev. E 77, 036702 (2008). [23] T. Lee and P. F. Fischer, Phys. Rev. E 74, 046709 (2006). [24] H. Liu, A. J. Valocchi, and Q. Kang, Phys. Rev. E 85, 046309 (2012).
063305-11
Multicomponent lattice Boltzmann equation method with a discontinuous hydrodynamic interface T. J. Spencer and I. Halliday* Materials and Engineering Research Institute, Sheffield Hallam University, Sheffield S1 1WB, United Kingdom (Received 4 October 2012; revised manuscript received 5 August 2013; published 6 December 2013) In the multicomponent lattice Boltzmann equation simulation method (MCLB), applied to the continuum regime of fluid flow, the finite width of the fluid-fluid interface introduces unphysical scales. We present a practical, robust, computationally efficient, and easy to implement solution to this problem which needs only low order interpolation to be stable and accurate and is applicable to any MCLB variant which uses a continuous phase field to distinguish between immiscible fluids with arrested coalescence. Our method extends the ideas of Kim and Pitsch, [Phys. Fluids 19, 108101 (2007)] and uses no external force distribution whatsoever to generate continuum interfacial physics, i.e., the Laplace law and no traction conditions on interfacial stresses. As such, it is amenable to the simplest form of Chapman-Enskog analysis used for lattice Boltzmann models. We assess our method and proceed to compare key results obtained with it against other equivalent data, obtained using the established continuum regime MCLB technique based upon the work of Lishchuk, Care, and Halliday, [Phys. Rev. E 67, 036701 (2003)] and Halliday, Hollis, and Care, [Phys. Rev. E 76, 026708 (2007)], quantifying performance in terms of the minimum feasible capillary available to simulation using that technique. DOI: 10.1103/PhysRevE.88.063305
PACS number(s): 47.11.−j
I. INTRODUCTION
Lattice Boltzmann method (LB) has advantages over conventional computational fluid dynamics in only a few areas, one being multicomponent flow. A facility with multicomponent flow is, alone, sufficient to justify the existence of the LB technique, however, for the flow of immiscible fluids with arrested coalescence, in the continuum hydrodynamic regime, it is of significance in a range of biological, microfluidic, and general engineering applications. The multicomponent lattice Boltzmann simulation method (MCLB) is a simple Eulerian numerical technique, in which a fluid-fluid interface is represented by a phase field, or order parameter, throughout denoted ρ N . Here we are concerned with the particular class of MCLB aimed at the continuum regime of hydrodynamics which uses a form of immersed boundary force [1], the action of which is directed by ρ N , to embed appropriate interfacial physics. The fact that ρ N cannot change discontinuously imposes extraneous length scales on MCLB simulations, characteristic of the distance over which the phase field changes between the values which represent the separated components. In mesoscopic applications, this is often acceptable but in large scale applications to continuum hydrodynamics it is inconvenient (since computational mesh is a scarce resource) and inaccurate (since the fluid-fluid interface is viewed as an infinitely thin surface, separating two fluids, upon which certain boundary conditions must be imposed [2]). In short, a diffuse hydrodynamic boundary obstructs detailed comparison with analytic calculations and calculations obtained using other numerical methods (such as boundary integral calculations), it introduces unphysical scales into simulations, and it is a drain on those computational resources which are associated with resolution. To address this problem we recount a means of imposing “discontinuous” interfaces between completely immiscible
* Corresponding author: [email protected]; FAX: (+44)114 225 3501.
1539-3755/2013/88(6)/063305(11)
LB fluids. For present purposes, we define a “discontinuous hydrodynamic interface” to be a closed, reducible, sublattice contour at which the separated fluids’ physical properties, such as pressure (or density) and viscosity but not, note, the phase field ρ N (which must remain differentiable) undergo a step change. The MCLB method we devise to achieve this is simple; it needs only low order interpolation to be stable and, judged in terms of the interfacial microcurrent (see below), the solution we advance introduces no penalty over the range of envisaged capillary numbers Ca to which MCLB may be accurately applied. Furthermore, our solution is computationally efficient and theoretically straightforward, since it requires no external force distribution, such as has been used previously [3,4] to generate appropriate continuum interfacial physics; the interface method advanced here more closely approximates a true boundary condition, which can be treated with the simplest form of Chapman-Enskog analysis [5]. In Sec. II selective background information is presented, upon which we shall build in Sec. III, where we introduce two closely related discontinuous hydrodynamic interface methods. In Sec. IV we will present results which we proceed to discuss and in Sec. V we present our conclusions and plans for further work. For the sake of simplicity, and to maintain parity with previous, related work [3,6], we consider the LB variant designated D2Q9 LBGK [7], which has lattice Q = 9 links denoted ci , i = 0–8, each with an associated weight denoted ti [5]. However, the work presented here will generalize very transparently to other multiple relaxation time LB variants, since no force distribution is used. We shall use common notation throughout. See Fig. 1 of Ref. [4] for the lattice link numbering convention in use. For simplicity, units of the lattice are used, so that x = t = 1 throughout. II. MCLB USING LISHCHUK’S METHOD
In this work, we rely upon certain aspects of our previous methods, i.e., that advanced by Lishchuk et al. in Ref. [3], combined with a distributed interface method after d’Ortona et al. [4,8]. We shall refer to this combination as “Lishchuk’s
063305-1
©2013 American Physical Society
T. J. SPENCER AND I. HALLIDAY
PHYSICAL REVIEW E 88, 063305 (2013)
method” [3,4] for brevity. In this section we review those of its features which are directly relevant and introduce related extensions required by the discontinuous hydrodynamic interface method. In the wake of the pioneering work of Gunstensen and Rothman [9], several two-component MCLB methods have arisen, each distinguished by the physical regime of its fluid-fluid interface. At the mesoscale, where the kinematics of phase separation must be considered, free energy methods [10,11] and the thermodynamically consistent extensions of Wagner et al. [12–14] (based, as they are, on the Cahn-Hilliard theory) are appropriate; for workers with a background in molecular simulation, the Shan-Chen method [15] is a natural choice. At the continuum scale, addressing completely immiscible fluids, the MCLB interface of Lishchuk et al. [3,4] provides a robust MCLB with assigned surface tension, appropriate dynamics, and, when used with an appropriate segregation [4], continuum interfacial kinematics. It is important to note that all the methods cited above contain distributed fluid-fluid interfaces, which extend over about >5 lattice sites (sites having −0.9 ρ N 0.9), with a variable body force applied on all of them, to produce a Laplace pressure step. Lishchuk’s method [3,4] imposes Laplace law and “no-traction” effects by impressing a curvature-dependent, immersed-boundary external force density, primarily in regions where fluid component phase parameter ρ N varies. The calculation of this force and the necessary correction to the model’s velocity [16] complicates the algorithm and is expensive. Let two fluid distributions which occupy lattice link i at position r and time t be described by distribution functions, Ri (r,t) and Bi (r,t), with the overall distribution function
3 (a)
ρ (lattice units)
2.5 2 1.5 1 0.5 0 -5
-4
-3
-2
-1 0 1 x (lattice units)
2
3
4
5
3 (b)
ρ (lattice units)
2.5 2 1.5 1 0.5 0 -5
-4
-3
-2
-1 0 1 x (lattice units)
2
3
4
5
3 (c)
ρ (lattice units)
2.5 2 1.5 1
fi (r,t) = Ri (r,t) + Bi (r,t) .
0.5
The nodal density of the red and blue fluid is R (r,t) = Ri (r,t) , B (r,t) = Bi (r,t) ,
0 -5
-4
-3
-2
-1 0 1 x (lattice units)
2
3
4
5
3
i
(d)
(2)
i
and the local phase field is [4]
2.5
ρ (lattice units)
(1)
ρ N (r,t) =
2 1.5
0.5 0 -4
-3
-2
-1 0 1 x (lattice units)
2
3
4
(3)
This definition is used without modification in the present work, to define an extended, smoothly varying phase field. Closed contour ρ N = 0 will define the discontinuous hydrodynamic interface, which, having no transverse extension, is always “sublattice.” Throughout the interfacial region, where |ρ N | < 1, a local interface normal is easily derived from derivative of ρ N :
1
-5
R (r,t) − B (r,t) . R (r,t) + B (r,t)
5
FIG. 1. The development of numerical instability as the interfacial pressure, or density step, , increases. Rest interface density profiles for red (triangles) and blue (squares) fluids obtained with various values of pressure step . The lines are included only to guide the eye. The interface is centered upon the sublattice position x0 = 0.5 lattice units, on a D2Q9 lattice, with density ρ0 = 2.0, interface segregation parameter β = 0.67, and LBGK collision parameter τ = 1. The fluids are at rest. Numerical differentiability degrades as increases. It is important to note that the stable system with = 0.4 depicted in [Fig. 1(b)] corresponds to a density step of 0%, which corresponds to a large Laplace law pressure step.
nˆ = −
∇ρ N . |∇ρ N |
(4)
The extended interfacial curvature field is obtained from the surface gradient of nˆ = (nˆ x ,nˆ y ), which, in two dimensions, is given by [3] ∂ nˆ y ∂ nˆ y ∂ nˆ x ∂ nˆ x − nˆ 2y + − nˆ 2x . (5) K(r) ≡ nˆ x nˆ y ∂x ∂y ∂x ∂y This field will be used to determine the curvature on the hydrodynamic interface ρ N = 0.
063305-2
MULTICOMPONENT LATTICE BOLTZMANN EQUATION . . .
PHYSICAL REVIEW E 88, 063305 (2013)
All the derivatives in Eqs. (4) and (5) may be computed to good accuracy by using the stencil ∂φ = 3 i ti φ(r + ci )ciα + o(x 3 ), (6) ∂xα r
(a)
cs2 ρ (lattice units)
in which φ(r) denotes ρ N , nˆ x , or nˆ y . The above result applies to the D2Q9 lattice. Shan has pointed out that the above stencil contains anisotropic lattice tensors in its error term [17]. Hence, normally directed forces calculated from it may contain tangential components which drive an interfacial microcurrent; we return to this issue in Sec. IV. The largest such anisotropic contribution in the error term of Eq. (6) may be eliminated, resulting in the following alternative stencil (again for the D2Q9 lattice): ∂φ = 3.2 i ti φ(r + ci )ciα ∂xα r (7)
u0 =
u(r )ρ N (r) − u(r)ρ N (r ) . ρ N (r) − ρ N (r )
ρ -1 -1.5
(c)
K0 (lattice units )
-1
|u| x10 (lattice units)
6
Figure 2(c) shows the value of curvature K0 interpolated using Eq. (10) on the contour ρ N = 0, for a circular drop of initial radius 35 lattice units. In Lishchuk’s method [3,4], continuum regime interfacial effects emerge by applying an external force density F = 1 ∇ρ N α0 K, where 12 ∇ρ N is the appropriate weight [3] and α0 2 is the interfacial tension. Impressing this force, one recovers a Laplace law pressure step [2] and the no-traction condition [3]. Inserting external forces, however, complicates LB’s underlying theory, its implementation, and reduces execution
0.03 0.025 0.02 0.015 0.01 0.005 0 -0.005 -0.01
(d)
(10)
(11)
0 -0.5
which, using Eq. (8) gives
A similar interpolation may be used to obtain the fluid velocity on the ρ N = 0 contour:
60
0.5
The local curvature K0 of the ρ N = 0 boundary, or fluid-fluid interface, at its intersection with ci may then be obtained from the interpolation: K(r ) − K(r) K0 = K(r) + d , (9) |ci | K(r )ρ N (r) − K(r)ρ N (r ) . ρ N (r) − ρ N (r )
40
1
We note that the phase field generates no hydrodynamic behavior whatsoever until some interaction between it and the LB distribution function fi is defined. It is for this essential reason that a continuous phase field may be used to produce a discontinuous hydrodynamic interface. In the present work, it will be necessary to evaluate the curvature at points on the contour ρ N = 0. This is achieved with sufficient accuracy, using a low order interpolation as follows. Identify all pairs of lattice sites r [with ρ N (r) > 0] and r [with ρ N (r ) < 0] using the criterion ρ N (r)ρ N (r ) < 0. Every such pair defines a lattice link ci which intersects the interface contour ρ N = 0 at distance d from the node at r, where |ci | d = N . (8) ρ N (r) ρ (r) − ρ N (r )
K0 =
20
140 120 100 80 60 40 20 80 100 120 140 0
(b)
N
− 0.1 i ti φ(r + 2ci )ciα + o(x ). 3
0.6668 0.66678 0.66676 0.66674 0.66672 0.6667 0.66668 0.66666 0.66664 0.66662 0.6666 0
3 2.5 2 1.5 1 0.5 0 -0.5
FIG. 2. For all these data the two-dimensional (2D) drop has initial radius of 35 lattice units. It was simulated on a lattice of 125 × 125 lattice units, without viscosity contrast, using LBGK collision parameter τ = 1.0, separation parameter β = 0.67, and an interfacial tension parameter α = 5.0 × 10−3 lattice units. All parts of this figure show nodal values, hence the smooth ρ N = 0 contour cannot be distinguished. The data represent the nodal value of steady state (a) discontinuous pressure field, (b) corresponding continuous phase field, ρ N (c) interpolated curvature value displayed at the lattice node closest to the ρ N = 0 contour, (d) interfacial microcurrent [note the scale factor in (d)]. The data of this figure are typical of both model variants described in Secs. III C1 and III C2.
speed, for to describe fluid subject to an external force density requires both the computation of that force and the insertion
063305-3
T. J. SPENCER AND I. HALLIDAY
PHYSICAL REVIEW E 88, 063305 (2013)
of a “source term” Fi into the evolution equation [16,18]: 1 † fi (r,t) = fi(0) (ρ,u) + 1 − fi(1) (∇ρ, . . . ∂x uy . . . ,F) τ (12) + Fi (τ,F,u) , 1 )( cic−u 2 2τ s
+ · F. Here we have conwith Fi = ti (1 − sidered the LBGK variant [7]; the dagger superscript indicates a postcollision, prepropagate quantity, fi(0) (ρ,u) denotes the equilibrium distribution function [5,7]:
uα uβ ciα ciβ − cs2 δαβ uα ciα (0) fi (ρ,u) = ti ρ 1 + 2 + , (13) cs 2cs4 (ci .u)(ci ) ) cs4
and density and velocity are defined as follows: 1 ρ≡ fi , ρu ≡ fi ci + F, 2 i i
τ together determine the viscous stress tensor, 1 i fi(1) (r)ciα ciβ , (r) = − 1 − σαβ 2τ
(20)
and that its weakly compressible hydrodynamics has a NavierStokes equation [20], D ∂ 2 ∂ ρuα = − cs ρ + σ . Dt ∂xα ∂xβ αβ
(21)
Note that lattice units have been used in Eqs. (19)–(21). III. A DISCONTINUOUS HYDRODYNAMIC INTERFACE FOR MCLB
(17)
Laplace law compliance and stress continuity boundary conditions require that two separated fluids’ hydrostatic pressure, p = cs2 ρ, evaluated on either side of a short section of interface, should differ by an amount α0 K0 , where K0 is the local interface curvature and α0 is the interfacial tension while tangential and normal viscous stress contractions should equate [2]. The kinematic condition of mutual impenetrability requires that the two fluids should move with the same velocity. Let us consider first the interface between two immiscible fluids, still designated red and blue, having identical shear viscosity, separated by an effectively discontinuous interface. Later we will consider viscosity contrast, i.e., a change in the value of collision parameter τ (which actually introduces no further complications). How are such dynamic and kinematic boundary conditions to be incorporated into a simple, efficient, and accurate MCLB simulation? It is appropriate to recall that there is, in fact, only a single “background” lattice fluid in all MCLB methods and that only physical observables like pressure and density need to change discontinuously.
in which parameter β (β 0.69, note) regulates phase field width. Furthermore, the red color density (say) of a locally flat interface with normal nˆ embedded between fluids in uniform motion u is described by [19] β ρ0 ˆ − ut) ; 1 − tanh R(x,t) = n.(r (18) 2 x
To comply with continuum physics, a LB distribution function which streams, or propagates, across a sublattice boundary ρ N = 0 must see only its equilibrium part, fi(0) (ρ,u), change by an amount which is consistent with the lattice fluid density change associated with the Laplace pressure step:
(14)
all other symbols having their usual meaning [4]. Streaming may be written †
fi (r + ci ,t + 1) = fi (r,t).
(15)
In the present work, we will obtain the same physics without recourse to an external force density F. Note that for F = 0, Fi = 0 and Eq. (12) simplifies considerably. Physical interfacial kinematics and isotropy emerge from carefully chosen postcollision color segregation rules. We have previously shown [4] that correct motion of the phase field boundary accrues from the use of d’Ortona’s segregation [8], R(r,t)B(r,t) R(r,t) † † ˆ Ri (r,t) = (16) f (r,t) + βti ci .n, ρ(r,t) i ρ(r,t) †
†
†
Bi (r,t) = fi (r,t) − Ri (r,t),
here x is lattice spacing. Recall that in this work x = 1. This result shows that the time development of a flat interface of the type used in this work corresponds to advection at the same speed as the underlying fluid, which is consistent with the kinematic condition of mutual impenetrability, and demonstrates the Galilean invariance of the color algorithm. Note that while the variation of R and B may be influenced by β, in Sec. III F the variation of ρ = R + B will be shown to be independent of β. It is relevant to recall that the LBGK method’s kinematic viscosity is determined by the collision parameter: 2τ − 1 η , (19) ν= = ρ 6 that the principal part of its distribution function, fi = Ri + Bi , is the equilibrium contribution fi(0) (ρ,u) which, note, is linear in the density ρ, as is clear in Eq. (13), that its nonequilibrium contribution to fi , fi(1) and collision parameter
A. Kinematic and dynamic compliance
δP = α0 K0 .
(22)
It will often be convenient to think of the corresponding density step: ≡ R0 − B0 ,
(23)
where R0 and B0 are the density of the red and blue fluids at rest, at a large distance from any interface. Clearly, =
α0 K0 , cs2
(24)
by the Laplace law. Any distribution function which streams across a sublattice boundary should see its equilibrium part change by an amount fi(0) (,u), where = α0 K0 /cs2 , while in the absence of viscosity and density contrasts, its nonequilibrium part should not change (so that viscous stresses and, hence, their normal and tangential contractions, will be continuous).
063305-4
MULTICOMPONENT LATTICE BOLTZMANN EQUATION . . . B. Background
cs2 ρ
The aim of adjusting the value of in accord with the Laplace law while maintaining f (1) continuous across our MCLB interface is facilitated by previous work, which we now summarize. To simulate pressure driven internal flows, Kim and Pitsch devised a periodic boundary closure which produces a spatially and temporally fixed pressure step of constant magnitude cs2 at, say, x = x0 , with the velocity and fluid stresses continuous across that boundary [21]. This closure may be thought of as the exchange of a constant “quantum” of density ti between appropriately chosen LB distribution function component pairs, defined as follows. Suppose lattice link ci (i > 0) starts from the lattice node at position r = (x,y) and crosses the x = x0 boundary and that ρ(r) = ρ(r + ci ) + , i.e., ρ(r) > ρ(r + ci ). LB distribution functions fi (r,t) and fj (r + ci ,t), with j = i ⊕ (Q − 1)/2, form a pair which counterpropagate along link ci , between times t and t + 1. Here, ⊕ denotes addition, modulo Q (Q = 9 here, of course). The distribution function pair fj (r + ci ,t) and fi (r,t) is characterized by the two values of r and i [j = i ⊕ (Q − 1)/2 is determined by the choice of i]. The “opposed” distribution function pair fi (r,t) and fj (r + ci ,t), chosen in this way, has equilibrium components determined solely by the density (pressure) and velocity of the fluid, which we denote fi(0) (ρ(r),u(r)), fj(0) (ρ(r + ci ),u(r + ci )),
(25)
still j = i ⊕ (Q − 1)/2, note. As we have indicated, to ensure that the required pressure difference cs2 exists between the fluid occupying lattice nodes r and r + ci it is necessary to enforce the pressure step condition: ρ(r) − ρ(r + ci ) = .
(26)
This can be achieved, Kim and Pitsch showed, by that notional exchange, during the LB propagate step, of a fixed quantum of density ti , implicit in the following replacements, imposed at every time step: fi (r,t) → fi (r,t) −
fi(0) (,u(r)),
fj (r + ci ,t) → fj (r + ci ,t) + fj(0) (,u(r + ci )).
(27) (28)
Note that the fi(1) component of the distribution is left intact. However, the total mass of the system is invariant only if fi(0) (,u(r)) = fj(0) (,u(r + ci )). We return to this point in Sec. IV. The above replacements are readily computed, using Eq. (13), for a chosen value of . To produce interfacial physics appropriate to a discontinuous hydrodynamic boundary between two immiscible fluids, designated red and blue, one can apply a version of the steps defined in Eq. (27), above. We present below two related algorithms, each of which does essentially this. C. Algorithm
Let us consider the interface of a two component MCLB model. The latter is henceforth defined as the closed, sublattice contour ρ N = 0 of the continuous phase field and is denoted C. Let us suppose that the red fluid is enclosed by the blue fluid and therefore subject to a Laplace pressure, that C has
PHYSICAL REVIEW E 88, 063305 (2013)
known local curvature K0 , that lattice node r lies inside C [so ρ N (r) > 0], and that lattice node r ≡ r + ci lies outside C [so ρ N (r + ci ) < 0]. The density of the red fluid differs from that of the blue fluid by an amount = α0 k0 /cs2 , determined by the Laplace pressure, and (for the moment) no viscosity contrast exists. 1. Mass conserving discontinuous hydrodynamic boundary
We use the following extension to the MCLB evolution rule set, at every time step, on all the distribution function component pairs “connected” by the link ci (which intersects C, recall): Ri (r,t) → Ri (r,t) − 12 fi(0) (,0) ,
(29)
Bi (r,t) → Bi (r,t) − 12 fi(0) (,0) ,
(30)
Rj (r ,t) → Rj (r ,t) + 12 fj(0) (,0) ,
(31)
Bj (r ,t) → Bj (r ,t) + 12 fj(0) (,0) .
(32)
Here, = α0 K0 /cs2 . K0 , the curvature on contour C, is interpolated using Eq. (10) from the curvature field defined by Eqs. (4) and (5). In Eqs. (4) and (5) all derivatives of ρ N and nˆ are evaluated using the stencil in Eq. (6). To a first approximation, it is appropriate that additional mass is allocated equally between red and blue components in (29)–(32) because, on C, the value of ρ N ≈ 0, hence the fluid there will be composed of red and blue components equally. The sum Ri (r,t) + Rj (r ,t) is invariant under replacements (29)–(32) which therefore conserve colored mass. There is, however, local momentum introduced between the pair of sites at r and r = r + ci . Replacements (29)–(32) introduce a net momentum: Pij (r,r ) = fi(0) (,0,t)ci + fj(0) (−,0,t)cj .
(33)
Now ci = −cj , weights ti = tj , and so, from Eq. (13): Pij (r,r ) = 2ti ci , where = a drop is therefore
α0 K0 /cs2 ,
(34)
recall. The total momentum imparted to
P = 2
ti ci .
(35)
r
Whatever its shape, a closed contour must cross any lattice line an even number of times. Therefore, the contributions to the local momentum r ti ci cancel in pairs and the total momentum change P(r,r ). Hence, our first, most computationally economical, colored mass conserving discontinuous hydrodynamic boundary method also conserves global (though not local) momentum. 2. Low microcurrent discontinuous hydrodynamic boundary
The more costly method of this section approximately conserves color mass and is based upon a more considered approach to allocating the transferred mass in the interfacial region. The interfacial microcurrent is a small unphysical velocity which arises in the region of the interface. This spurious signal is substantially reduced when the following, modified
063305-5
T. J. SPENCER AND I. HALLIDAY
PHYSICAL REVIEW E 88, 063305 (2013)
extension to the MCLB evolution rule set is used at every time step, to modify the values of the “opposed,” counterpropagating distribution function component pairs, connected by the link ci : Ri (r,t) → Ri (r,t) −
g(i,r)fi(0)
(,u0 ) ,
(36)
Bi (r,t) → Bi (r,t) − g (i,r)fi(0) (,u0 ) ,
(37)
Rj (r ,t) → Rj (r ,t) + g(j,r )fj(0) (,u0 ) ,
(38)
Bj (r ,t) → Bj (r ,t) + g (j,r )fj(0) (,u0 ) .
(39)
The factors g(i,r) and g (i,r) which appear above will be discussed shortly. In these replacements, = α0 K0 /cs2 . K0 , the curvature on contour C, is interpolated using Eq. (10). u0 denotes the velocity interpolated onto the ρ N contour C, again using Eq. (10). It is important to note that the underlying curvature field is now calculated from Eqs. (4) and (5) but now the derivatives of ρ N and nˆ are evaluated using the modified stencil in Eq. (7). The factors g(r) and g (r) appearing in replacements (36)–(39) are defined Ri (r,t) , g(i,r) ≡ Ri (r,t) + Bi (r,t)
Bi (r,t) g (i,r) ≡ . Ri (r,t) + Bi (r,t) (40)
They are designed to conserve the value of the phase field ρ N in the adjusted link populations, as follows. Consider the node at r . It sends information in the form of Rj (r ,t) and Bj (r ,t) along link cj , from node r , where the value of the phase field is ρ N (r ). However, Rj (r ,t) and Bj (r ,t) are adjusted by replacements (36)–(39). Temporarily drop the variables t and r . In order to preserve the perceived value, at node r, of the phase field ρ N (r ), the changes in the population, δRj , δBj , should satisfy the condition ρjN ≡
Rj − B j Rj + δRj − (Bj + δBj ) = . Rj + B j Rj + δRj + Bj + δBj
(41)
Taking a binomial expansion, we obtain δRi − δBj δRj + δBj = + o δRj2 , Rj − B j Rj + B j
(42)
which may be solved δBj =
Bj δRj , Rj
(43)
and since δRj + δBj = fj(0) (,u0 ) it is straightforward to obtain Rj fj(0) (,u0 ), δRj = Rj + B j Bj δBj = fj(0) (,u0 ), (44) Rj + B j which explains the definition of the factor g(j,r ) appearing in replacements (36)–(39). Clearly, colored mass is now not conserved, and small fluctuations in mass arise. We return to this issue in Sec. IV. The net momentum introduced between the pair of sites at r and r = r + ci by the replacements (36)–(39) is given by
the expression Pij (r) = fi(0) (,u0 )ci + fj(0) (−,u0 )cj .
(45)
Again, ci = −cj and the link weights equate, ti = tj , so, taking a linear approximation from Eq. (13), (u0 − u0 ) · ci ci + o u20 , Pij (r,r ) = 2ti ci + ti 2 cs (46) where = α0 K0 /cs2 , recall. The total momentum imparted to a drop is therefore (47) ti ci + o u20 . P = 2 r
Whatever its shape, a closed contour must cross any lattice line an even number of times. Therefore, the contributions to r ti ci again cancel in pairs and the principal contribution to the total momentum, above, is zero to o(u20 ). D. Viscosity contrasts
Thus far, we have considered immiscible fluids with identical viscosity, separated by an interface located on a sublattice contour C, i.e., on the closed curve ρ N = 0. Let us now allow a change in fluid viscosity across contour C. Thus, the red (blue) fluid occupying lattice node r (r = r + ci ) has a value of kinematic viscosity determined by τr (τr = τb ). Physical fluid with a shear viscosity η, which varies, possibly in response to local heating, is described by Eq. (21) [2], which shows that viscous stress must be continuous (i.e., differentiable) in the fluid, given that velocity and density are also continuous. Likewise, in LBGK dynamics, viscous stress must be continuous in the presence of inhomogeneities in τ if the velocity and density are continuous. In the present context, velocity across the region of the discontinuous interface of Sec. III A is continuous by design [21] and since we only effect a step change in pressure (density) at the sublattice C, long wavelength variations in density are also continuous. It follows that the viscous stress is continuous throughout the interfacial region. Put another way, a rapid change in shear viscosity η will produce compensating changes in the shear rate sαβ but the viscous stress σαβ will remain continuous. We return to this point in Sec. IV shortly. Given the continuity of the viscous stress [20], 1 i fi(1) (r)ciα ciβ , (48) σαβ (r) = − 1 − 2τ in the presence of a sudden change in viscosity, i.e., LBGK collision parameter τ , interfacial viscous stress contractions ˆ σαβ tα , with tˆ the assumed smoothly varying unit tangent vector, must also be continuous. Clearly, this is a necessary condition for the emergence of the appropriate hydrodynamic boundary (red) (blue) conditions, such as σαβ tˆβ − σαβ tˆβ = 0. E. Density contrasts
The additional algorithmic extension we have used previously to amplify density contrasts [22] can be applied for the present model. Of course, compounding algorithmic complexity must have a range of adverse consequences.
063305-6
MULTICOMPONENT LATTICE BOLTZMANN EQUATION . . .
However, the addition of the stress insulating forces postulated previously is perfectly feasible. For the sake of brevity we reserve this investigation for a future publication. F. Analysis
Consider a flat MCLB fluid-fluid interface orientated perpendicular to the x axis, maintained by the segregation rule in Eq. (16), embedded in a uniform fluid of constant density ρ0 . By Taylor expanding the discrete kinetic equation of the MCLB and considering the steady state, we have shown elsewhere [4] that such an interface when centered upon x = x0 [i.e., with ρ N (x) < 0,x < x0 ], will be described by the ordinary differential equation, β d 2 d2 β d R+2 R = 0, R−2 dx 2 h dx hρ0 dx
(49)
in which h is the lattice spacing. Note that this equation has a first integral. Throughout this work we have taken h = 1, with which simplification the solution is [4,19] 1 + tanh (β(x − x0 )) R(x) = ρ0 , 2 1 − tanh (β(x − x0 )) . (50) B(x) = ρ0 2 Here, it has been assumed that the region x < x0 is initially occupied by blue fluid. Simulation data support this solution. We seek now to determine the corresponding red fluid density variation in the presence of a density discontinuity: ρ(x) = ρ0 + δ[1 − H (x − x0 )],
(51)
where H (x) is the Heaviside step function. Noting that the point x0 will never correspond to a point lattice location, one can follow an identical approach to that used to derive Eq. (49). An additional assumption is that for x < x0 (x > x0 ) ρ(x) = ρ0 + (ρ0 ) is used to encapsulate the boundary condition implicit in the rule set of Eqs. (29)–(32), straightforwardly to show that the introduction of a density discontinuity at x = x0 modifies only one term in the ordinary differential equation (49): β d 2 d2 β d R+2 R = 0. R−2 2 dx h dx hρ(x) dx
that the phase field does not contain any discontinuity, ρ N (x) ≡
Note that this solution meets those density boundary conditions implicit in our replacement expressions (29)–(32), namely Eq. (51), (55)
R(x) − B(x) = tanh (β(x − x0 )), R(x) + B(x)
(56)
and that the latter will therefore remain differentiable, provided certain straightforward conditions are met. In practice, to obtain curvatures K, the phase-field profile ρ N must be numerically differentiable, without taking elaborate precautions (which would undermine the simplicity of the method). Implicit in the numerical differentiability of ρ N is the differentiability of R(x) and B(x). For R(x) to be differentiable, the step change in R(x), which occurs at x = x0 and is of size /2, must be small compared with the change in the continuous contribution tanh (β(x − x0 )) in the region of x0 . This leads directly to the condition d ρ0 − tanh (β(x − x0 )) (57) . dx 2 2 x=x0 Performing the differentiation and noting that cosh2 (β(x − x0 )) = 1 when x = x0 we obtain β(ρ0 − ) which may be re-arranged to provide the following stability criterion:
βρ0 . 1+β
(58)
Figure 1 shows the steady-state color density profiles associated with a flat interface y = const, centered upon the sublattice position x0 = 0.5 lattice units, for a D2Q9 lattice, with density ρ0 = 2.0, a range of values of the density step , segregation parameter β = 0.67, and LBGK collision parameter τ = 1. The fluids are at rest. From inequality (58), the approximate value at which the phase field ρ N ceases to become numerically differentiable is 0.8, which agrees with the trend in the figure. The existence of this instability must not be overemphasized; a density step = 0.4 [Fig. 1(b)] corresponds to a 20% density step and hence to a large Laplace law pressure step and interfacial tension. Moreover, Eq. (58) shows that it may be controlled by an appropriate choice of ρ0 and β. The latter is not now subject to restrictions placed upon it using immersed boundary techniques, as we shall now see.
(52)
For x = x0 , ρ(x) is constant and it is possible to write down a solution to Eq. (52) from Eq. (49), by replacing ρ0 with ρ(x). Using Eq. (51), we obtain ρ0 + R(x) = − H (x − x0 ) [1 + tanh (β(x − x0 ))] , 2 2 (53) ρ0 + B(x) = − H (x − x0 ) [1 − tanh (β(x − x0 ))] . 2 2 (54)
ρ(x) ≡ R(x) + B(x) = ρ0 + [1 − H (x)],
PHYSICAL REVIEW E 88, 063305 (2013)
IV. RESULTS AND DISCUSSION
The microcurrent, or spurious velocity, associated with MCLB interfaces is a small but unphysical circulation in the near-interfacial region occurring even when the separated fluids are otherwise at rest. See Fig. 2(d). We shall denote the maximum microcurrent velocity U . The vorticity created at an interface will diffuse, at a rate controlled by the fluid viscosity, into the bulk fluids where it may accumulate, unless balanced by dissipation. Lee et al. have successfully reduced microcurrents in free energy MCLB methods [23] and Shan [17] has argued that they are, in the Shan-Chen model [15], caused by lattice anisotropy and that they are therefore expensive to eliminate. It is not our aim here to carefully consider the origin of this essentially numerical artifact. We adopt a more practical approach. Strong flow obscures the microcurrent. Conversely, microcurrent activity obscures physical flow when the latter is
063305-7
T. J. SPENCER AND I. HALLIDAY
PHYSICAL REVIEW E 88, 063305 (2013)
Min(Ca) ≡
ηU (2τ − 1)U = . α0 6α0
(59)
By measuring a steady state U , then evaluating Min(Ca) defined above, it is possible to estimate the value of Ca at which, for a given level of resolution, the physical flow will cease to be distinguishable at all points in a simulation. A variety of other factors including lattice resolution and the value of the collision parameter τ also affect U . We fix lattice resolution for all the results presented in this section. We consider throughout a red drop of constant radius R = 35 lattice units, at rest at the center of a lattice of size 125 × 125 lattice units, subject to periodic boundary conditions. The systems simulated are therefore infinite in extent and contain periodic images of the red drop. Accordingly, the velocity induced by the microcurrent vorticity can increase for a considerable time. Throughout this section data obtained after 1.0 × 106 lattice time steps is designated steady state. The initial density of both red and blue fluids was set to a value 2.0. The hydrodynamically discontinuous interfaces developed here may be assessed by comparing their Min(Ca) with the corresponding quantity obtained from otherwise identical simulations using Lishchuk’s method. Let us consider first the case of zero viscosity contrast, setting τr = τb = 1. For the steady state of a drop, using segregation parameter β = 0.67 and a large range of interfacial tensions α0 , Fig. 3 shows the variation in Min(Ca), for the mass conserving discontinuous interface method of Sec. III C1 compared with Lishchuk’s method [3,4](continuous line). Figure 4 shows the equivalent variation in Min(Ca), for the low microcurrent discontinuous interface method of Sec. III C2, again compared with Lishchuk’s method [3,4] (continuous line). See the figure captions for a summary of the simulation parametrizations. From Fig. 3, it is clear that for α0 0.02 Lishchuk’s method [3,4] has the practical advantage, However, for smaller values of α0 a discontinuous interface may be used without apparent penalty. Indeed, it may be used to advantage. In these data, as α0 approaches zero, Min(Ca) approaches a constant, which indicates the existence of a regime where U ∝ α0 . More unexpected is the existence of an optimum value of α0 for the mass conserving hydrodynamically discontinuous
-2
10
min (Ca)
10-3
10-4
-5
10
-6
10
-5
10
-4
10
-3
-2
10 10 α0 (lattice units)
-1
10
0
10
FIG. 3. For all the data in this figure, the drop was assigned an initial radius of 35 lattice units, it was simulated on a lattice of 125 × 125 lattice units, without viscosity or density contrast, using LBGK collision parameter τ = 1.0 for both fluids, separation parameter β = 0.67, and a large range of interfacial tensions α0 . The solid squares show the variation of Min(Ca) for the mass conserving discontinuous hydrodynamic interface method, in Sec. III C1. For reference, the continuous line shows the value of Min(Ca) for Lishchuk’s method [3,4], which is an immersed boundary technique.
interface method. The data in Fig. 4 shows that the effect of the microcurrent is indeed reduced and that our low microcurrent discontinuous hydrodynamic interface method actually performs better than Lishchuk’s method over the range of interfacial tension parameter investigated. The value of segregation parameter β may be decreased, to make the phase field boundary more diffuse and, hence, to reduce the microcurrent. This exacerbates the problem of unphysical scales in the case of Lishchuk’s method. However, while reducing the value of β makes the phase field vary over a greater distance, the ρ N = 0 contour remains sharp and the -2
10
10-3 min (Ca)
weak. The length scale associated with the interface has been eliminated from the hydrodynamic signals’ spatial variation by the methods related above. Therefore, where flow is strong, any penalty in the form of an increase in U will be unimportant and the discontinuous hydrodynamic interface is clearly the preferred method. However, in weak flow the practical utility of the discontinuous interface method may be restricted by an increased in U . Hence, our current methods can be succinctly evaluated by assessing their maximum microcurrent velocity compared with an equivalent implementation of Lishchuk’s method, in the context of weak flow. While the relationship may not be linear, as the interfacial Laplace pressure step increases (in response to an increase in interfacial tension α0 ) so does the spurious velocity U . To take account of this effect, we characterize a given MCLB method by the dimensionless minimum capillary number accessible to simulation using that MCLB interface. Our minimum capillary number is defined as follows:
-4
10
10-5 -6
10
-7
10
-5
10
-4
10
-3
-2
10 10 α0 (lattice units)
-1
10
0
10
FIG. 4. For all the data in this figure, the drop was assigned an initial radius of 35 lattice units, it was simulated on a lattice of 125 × 125 lattice units, without viscosity or density contrast, using LBGK collision parameter τ = 1.0 for both fluids, separation parameter β = 0.67, and a large range of interfacial tensions α0 . The solid squares show the variation of Min(Ca) for the low microcurrent discontinuous hydrodynamic interface method of Sec. III C2. For reference, the continuous line shows the value of Min(Ca) for Lishchuk’s method [3,4].
063305-8
MULTICOMPONENT LATTICE BOLTZMANN EQUATION . . .
PHYSICAL REVIEW E 88, 063305 (2013)
-2
10-2
10
-3
10-3
10-4
10
-4
min (Ca)
min (Ca)
10
-5
10
10-5
10-6
10-6
10-7
10
-7
-8
10
10-5
10-4
10-3 10-2 α0 (lattice units)
10-1
10-8 -5 10
100
FIG. 5. For all the data in this figure, the drop was assigned an initial radius of 35 lattice units, it was simulated on a lattice of 125 × 125 lattice units, without viscosity or density contrast, using LBGK collision parameter τ = 1.0 for both fluids, separation parameter β = 0.4, and a large range of interfacial tensions α0 . The solid squares show the variation of Min(Ca) for the low microcurrent discontinuous hydrodynamic interface method of Sec. III C2. For reference, the continuous line shows the value of Min(Ca) for Lishchuk’s method [3,4] with β = 0.4. Note that using width of the interface in the latter is now considerable; these data should be considered only as reference data.
discontinuous hydrodynamic interface method may be used with β = 0.4 without penalty. Figure 5 compares Lishchuk’s method with our low microcurrent method, which again outperforms the latter, over the range of interfacial tension parameter investigated. See the figure caption for a summary of the simulation parametrization. It is very important to note that Figs. 3–5 show that the low microcurrent version of our algorithm with β = 0.4 considerably outperforms Lishchuk’s method with β = 0.69. Figure 6 shows data similar to that in Fig. 4, obtained with our low microcurrent method. The sole difference is the value of LBGK collision parameter τ = 5.56, corresponding to a tenfold increase of both fluids’ kinematic viscosities. Again, the solid squares show the variation of Min(Ca) for the low microcurrent discontinuous interface and, as in Fig. 3, the broken line shows the value of Min(Ca) for Lishchuk’s method [3,4]. Color masses are not conserved for the systems of Figs. 4 and 5. It is important to quantify that mass fluctuation. Suppose a red drop of initial red mass R0 has, after T time steps, a red mass RT . Let L denote the length of its interface, to which quantity the mass fluctuation must be proportional. The quantity , defined below, measures the fractional rate of change in red (drop) mass, expressed as a percentage of the original mass, R0 , per unit length of interface, RT − R0 100 . (60) ≡ T R0 L By fitting simulation data, over the range of interfacial tensions, α0 < 0.01, investigated the following empirical model for was obtained: = aα02 + bα0 + c,
(61)
-4
10
-3
-2
10 10 α0 (lattice units)
-1
10
0
10
FIG. 6. For these data, the drop was assigned an initial radius of 35 lattice units; it was simulated on a lattice of 125 × 125 lattice units, again without viscosity or density contrast, using an increased LBGK collision parameter τ = 5.56 (viscosities) for both fluids. The separation parameter β = 0.67 and a large range of interfacial tensions α have been used. The solid squares show the variation of Min(Ca) for the discontinuous hydrodynamic interface and, as in Fig. 3, the broken line shows the value of Min(Ca) for Lishchuk’s method [3,4]. The data for that method are little changed from those in Fig. 3. By contrast, those obtained from a discontinuous interface show an increase in microcurrent activity over the whole range of α0 , especially for larger values.
with coefficients a = 3.70 × 10−6 , b = 4.90 × 10−7 , and c = 5.12 × 10−9 , obtained by least squares fitting and characteristic of LB model and dimensionality. Clearly these values are very small indeed. For a small drop with radius 20 lattice units (and, proportionally, a long interface) and large interfacial tension α = 0.1 lattice units, a two-dimensional (2D) simulation using the discontinuous hydrodynamic interface method of Sec. III C2 would need to run for approximately 1 × 105 time steps before its mass would have changed by 1%. Let us next consider a viscosity contrast between the red and blue fluids. For simplicity, we consider a flat interface x = 0 and shear the fluids in the y direction, in the tangential direction. The red fluid occupies the region x 0 and has a kinematic viscosity νr = 10 (i.e., τr = 30.5). The blue fluid occupies the region x > 0 with νb = 1 (τb = 3.5). Hence, a viscosity contrast of 10 exists between the fluids. Both fluids are initially at rest. In Fig. 7, which concerns the near interfacial region and not the entire flow domain, note that the variation of the computed MCLB model’s viscous stress, evaluated from Eq. (20), is expressed as a fraction of its steady state value (solid circles). Solid diamonds correspond to t = 250 time steps, solid squares correspond to t = 500 time steps, and the solid triangles correspond to t = 1000 time steps. The viscous stress is apparently continuous at all times. Figure 8 shows an excerpt of the steady state velocity profile of the system. The steady state shear rates in the different fluids are in the expected ratio. Numerical expressions for gradients in ρ N are used in the ˆ present method to determine the interface normal vector n. Hence, the arguments of Shan [17] regarding the origin of the microcurrent do apply, albeit indirectly, for the present method appears to avoid use of explicit force. An effective
063305-9
T. J. SPENCER AND I. HALLIDAY
PHYSICAL REVIEW E 88, 063305 (2013)
two methods was observed to arise in two distinct contributions. First, the introduction of the factors g(i,r) and g (i,r), defined in Eq. (40), does not affect qualitative behavior but it does dampen the numerical instability shown in Fig. 1. That figure was obtained with our low microcurrent model; the equivalent data obtained for the mass conserving variant show the instability occurring for smaller values of δ. Second, the change in template used for the calculation of the derivatives of continuous phase field ρ N , from that in Eq. (6) to that defined in Eq. (7), produces a further, independent reduction in microcurrent activity. Note that more extensive use of this stencil was observed to produce only small benefits. We reserve further comment for our conclusions.
2.5
σxy(t)/σxy0
2 1.5 1 0.5 0 -6
-4
-2 0 2 y (lattice units)
4
6
FIG. 7. Variation of the viscous stress, evaluated from Eq. (20), expressed as a fraction of its steady value (solid circles) about a flat interface, at which the kinematic viscosity changes. For x 0, ν = 10 (τr = 30.5), x > 0, ν = 1 (τb = 3.5), so that a viscosity contrast of 10 exists. Diamonds correspond to t = 250 time steps, solid squares correspond to t = 500 time steps, and the solid triangles correspond to t = 1000 time steps. The viscous stress is apparently continuous at all times.
force, or momentum change [note Eq. (46), above], clearly accompanies any process that does not explicitly conserve momentum. This fact, and similarity between the level of microcurrent activity measured with the present method and that observed when using Lishcuk’s method [3,4], across a range of α0 , adds further support to the relevance of Shan’s theory to both models considered in this work. In summary, Figs. 3–5 demonstrate that the mass conserving discontinuous hydrodynamic variant, defined in Sec. III C1, to be precise in steps (29)–(32), outperforms our previous models, until large interfacial tensions create Laplace pressure, or density, steps which are associated with that numerical instability identified in Fig. 1; we return to this point shortly. The low microcurrent variant presented in Sec. III C2 and summarized in steps (36)–(39) performs still better and, furthermore, is stable for a larger range of interfacial tensions. The improvement between the performance of our 1 0.98 0.96 ux/uxmax
0.94 0.92 0.9 0.88 0.86 0.84 0.82 -6
-4
-2
0
2
4
6
x (lattice units)
FIG. 8. Steady-state shear rate corresponding to the system used to obtain the data for Fig. 7. The steady-state shear rates are in the expected ratio.
V. CONCLUSIONS AND FURTHER WORK
Flow calculations in the arrested coalescence regime of continuum hydrodynamics are facilitated by the multiple immiscible fluids lattice Boltzmann equation method introduced here. The method removes unphysical interfacial scales from the hydrodynamic signals and simplifies both the computational algorithm and theoretical framework. Our “discontinuous hydrodynamic interface” method still relies upon an auxiliary phase field which changes continuously which limits “direct” computational cost savings, as we discuss below. However, it still provides indirect cost savings in applications where, e.g., computational mesh resolution is limited and comparisons with boundary integral or analytical calculation are sought by removing the need to simulate at larger than necessary system sizes to achieve the narrow interface limit. Put another way, by using the discontinuous interface method advanced here, the simulator is always in the narrow interface, hydrodynamic limit. The discontinuous hydrodynamic interface method is simple and easy to implement; it needs only low order interpolation to be stable and accurate over a range of interfacial tensions and it can be generalized to three dimensions and to multiple relaxation time schemes easily, since it relies upon no external force distribution whatsoever. We have assessed it against our existing approach (designated here as Lishchuk’s method [3,4]) which is also adapted to the continuum regime. When considering two fluids of identical density (but possibly different viscosity) separated by an interface which is subject to the Laplace law and no-traction conditions on interfacial stresses, our method introduces a small gain when assessed in terms of a simulation’s minimum accessible capillary number over a range of interfacial tension parameters when flow is weak. For very large interfacial tensions however, our methods produce microcurrent activities in excess of those produced in Lishchuk’s method, which is important when the simulated flow is challengingly weak. Matters are less ambiguous where flow is strong compared with the microcurrent activity; here the present method will be seen to have the clear advantage of physical accuracy over Lishchuk’s method. A natural question arises regarding further work to incorporate, in additional to a Laplace pressure step, a fluid density step at the interface. While such a density step is simple to directly incorporate, only small factors of density of order unity are possible, due to the numerical instability shown in Fig. 1 and considered in Sec. III F. To facilitate careful,
063305-10
MULTICOMPONENT LATTICE BOLTZMANN EQUATION . . .
PHYSICAL REVIEW E 88, 063305 (2013)
fully 3D comparisons with similar techniques [24], it will be necessary to address the matter of large density difference by adapting previous work [22] and to extend the model’s dimensionality. Assessed directly, the discontinuous interface method is only a little more computationally efficient than a comparable force-density based scheme and its principal advantages remain the removal of the unphysical scales from continuum multicomponent flow calculations and the simplification of the multicomponent lattice Boltzmann equation method algorithm. The introduction of a sharp physical interface is achieved but a diffuse auxiliary phase field and the associated processing cost remain to be found. This raises another question in respect of the extent of any direct saving in computational cost. Direct cost savings are realized because (i) a distributed immersed boundary force density is not, in the present method, impressed throughout the interfacial
region and (ii) no associated velocity field corrections are required [16] (our models are amenable to the simplest form of Chapman-Enskog analysis [5] used for lattice Boltzmann models [20]). Put another way, all direct computational costs associated with the processing of the phase field remain. Hence we must be cautious in making claims for the computational savings associated with the algorithmic simplifications presented here. Where high resolutions are essential, an extended phase field boundary region will be required and the computational savings are, in relative terms, small but they improve as resolution decreases, for one is always in the narrow interface limit. Stressing that, in general, direct computational cost savings must depend upon the extent of interface in a particular calculation, as well as other factors like geometry, we can tentatively quantify the computational cost savings associated with the results of Sec. IV at about 10% reduction in execution time.
[1] C. S. Peskin, J. Comput. Phys. 220, 20 (1977). [2] L. Landau and E. M. Lifshitz, Fluid Mechanics, 2nd ed. (Pergamon, New York, 1987). [3] S. V. Lishchuk, C. M. Care, and I. Halliday, Phys. Rev. E 67, 036701 (2003). [4] I. Halliday, A. P. Hollis, and C. M. Care, Phys. Rev. E 76, 026708 (2007). [5] S. Succi, The Lattice Boltzmann Equation for Fluid Mechanics and Beyond (Oxford University Press, Oxford, 2001). [6] T. J. Spencer, I. Halliday, and C. M. Care, Phys. Rev. E 82, 066701 (2010). [7] Y. H. Qian, D. d’Humi`eres, and P. Lallemand, Europhys. Lett. 17, 479 (1992). [8] U. D’Ortona, D. Salin, M. Cieplak, R. Rybka, and J. Banavar, Phys. Rev. E 51, 3718 (1995). [9] A. K. Gunstensen, D. H. Rothman, S. Zaleski, and G. Zanetti, Phys. Rev. A 43, 4320 (1991). [10] M. R. Swift, W. R. Osborn, and J. M. Yeomans, Phys. Rev. Lett. 75, 830 (1995). [11] M. R. Swift, E. Orlandini, W. R. Osborn, and J. M. Yeomans, Phys. Rev. E 54, 5041 (1996).
[12] A. J. Wagner, Phys. Rev. E 74, 056703 (2006). [13] Q. Li and A. J. Wagner, Phys. Rev. E 76, 036701 (2007), and references therein. [14] A. J. Wagner and C. M. Pooley, Phys. Rev. E 76, 045702(R) (2007), and references therein. [15] X. W. Shan and H. D. Chen, Phys. Rev. E 49, 2941 (1994). [16] Z. Guo, C. Zheng, and B. Shi, Phys. Rev. E 65, 046308 (2002). [17] X. Shan, Phys. Rev. E 73, 047701 (2006) [18] X. He, S. Chen, and G. D. Doolen, J. Comput. Phys. 146, 282 (1998). [19] I. Halliday, S. V. Lishchuk, T. J. Spencer, G. Pontrelli, and C. M. Care, Phys. Rev. E 87, 023307 (2013). [20] S. Hou, Q. Zou, S. Chen, G. D. Doolen, and A. C. Cogley, J. Comput. Phys. 118, 329 (1995). [21] S. H. Kim and H. Pitsch, Phys. Fluids 19, 108101 (2007). [22] S. V. Lishchuk, I. Halliday, and C. M. Care, Phys. Rev. E 77, 036702 (2008). [23] T. Lee and P. F. Fischer, Phys. Rev. E 74, 046709 (2006). [24] H. Liu, A. J. Valocchi, and Q. Kang, Phys. Rev. E 85, 046309 (2012).
063305-11
