Critical properties and phase separation in lattice Boltzmann fluid mixtures

Nicos S. Martys1 and Jack F. Douglas2

1

Building Materials Division, Building and Fire Research Laboratory, National Institute of Standards and Technology,

Gaithersburg, Maryland 20899 2

Polymers Division, Material Science Engineering Laboratory, National Institute of Standards and Technology,

Gaithersburg, Maryland 20899

Received 4 May 2000; published 27 February 2001

Basic equilibrium properties of lattice Boltzmann LBfluid mixtures coexistence curve, surface tension,

interfacial profile, correlation lengthare calculated to characterize the critical phenomena occurring in these

model liquids and to establish a reduced variable description allowing a comparison with real fluid mixtures.

We observe mean-field critical exponents and amplitudes so that the LB model may be useful for modeling

high molecular weight polymer blends and other fluid mixtures approximated over a wide temperature range by

mean-field theory. We also briefly consider phase separation under quiescent and shearing conditions and point

out the strong influence of interacting boundaries on the qualitative form of the late-stage phase-separation

morphology.

DOI: 10.1103/PhysRevE.63.031205 PACS numbers: 47.11.j, 05.20.Dd, 02.70.c

I. INTRODUCTION

In many applications involving materials processing and

development, it is necessary to understand and control the

morphology of multiphase fluid mixtures and particulate dispersions subject to a complex flow history. These applications often involve free liquid-air boundaries that can respond to flow and phase-separation processes, solid

boundaries that can be preferentially wet by certain liquid

components, thin-film geometries, complex solid substrate

geometries, and high Reynolds number flows in which fluid

inertia is important. The development of computational

methods of sufficient flexibility and generality to treat such

realistic fluid dynamics problems is a basic theoretical challenge.

The lattice Boltzmann LBmethod and other related

computational methods based on cellular-automata ideas

e.g., lattice gas1 have emerged as powerful tools for

modeling complex fluid dynamics problems. These methods

are developing rapidly in response to recent theoretical advances and the availability of resources for large-scale computation. Applications of LB to modeling high Reynolds

flow 1, the dynamics of fluid phase separation 2,3, and

multicomponent fluid flow in porous media 4–6 have

proven the potential of LB as a general purpose computational scheme for modeling complex fluid dynamics problems. Many of these exploratory studies have emphasized the

development of the LB methodology and have not considered a direct comparison to the properties of real liquids.

Therefore, basic characteristics of these computational models of liquids are still largely unknown. In this paper, we

characterize the type of critical phenomena observed in two

LB models of multicomponent liquids. The first model, due

to Shan and Chen 2, modifies the Boltzmann equilibrium

distribution to account for fluid-fluid interactions. The second approach, considered in the present paper, incorporates

the forcing between two fluids directly into the body-forcing

term of the Boltzmann equation 7. In Appendix A, we

show how the second method is related to a density-gradient

expansion of a BBGKY collision operator 8–10, which

should facilitate comparison with the Cahn-Hilliard theory of

phase separation 11–16 and the Ginzburg-Landau theory of

critical phenomena 17,18.

We calculate basic equilibrium properties coexistence

curve, interfacial width and correlation length, surface tensionand express our results in terms of a reduced variables

description that allows comparison with real fluid measurements. Asymmetry in the mass and volume of the fluid components is considered in this comparison since this property

is characteristic of real liquids. The effect of flow and interacting boundaries on the phase separation process is briefly

explored to identify some basic phenomena of experimental

interest.

II. BRIEF REVIEW OF LATTICE BOLTZMANN MODEL

OF FLUID PHASE SEPARATION

The LB method of modeling fluid dynamics is actually a

family of models with varying degrees of faithfulness to the

properties of real liquids. These methods are currently in a

state of evolution as the models become better understood

and corrected for various deficiencies. In the present work,

we utilize two LB models of complex fluids. The first and

primarily studied method in this paper was proposed by Shan

and Chen 2,5. It is particularly simple in form and adaptable to complex flow conditions such as the presence of

solid-fluid and air-fluid boundaries. For comparison, a second approach is studied 7 that directly incorporates the

fluid/fluid interaction into a body-force term. This approach

removes second-order contributions with respect to the forcing see below. The physical basis of such an approach is

further described in Appendix A. A third approach, not studied in this paper, is due to Yeomans et al. 3,16 and is

strongly rooted in the Cahn-Hilliard model of binary mixtures 12–16. A major criticism of this approach is the lack

of Galilean invariance, but recent work suggests that the errors involved can be controlled in the description of the

phase separation of fluids in the absence of shear 16. It

PHYSICAL REVIEW E, VOLUME 63, 031205

63 031205-1
remains to be seen whether these higher-order effects cause

problems under conditions of fluid flow. All of these models,

however, are not strictly energy conserving, which should be

important under conditions of deep quenches and high rates

of flow. Recent work 10 has introduced a framework for

overcoming this technical problem and allows for systematic

improvements of the LB calculation of fluid properties. In

the present paper, we focus on isothermal flows and will

examine the effect of energy conservation on the two-phase

region of LB fluid mixtures in a future work.

We now present a brief description of the LB methods

used in this study. A considerable computational advantage

in modeling the fluid is obtained by restricting the particle

positions to the sites of a lattice. The LB method extends the

standard ‘‘quasicrystalline’’ fluid model 19–22 of classical

statistical mechanics by specifying the particle velocity distribution at each particle position so that both equilibrium

and dynamical properties of the fluid can be calculated. Macroscopic variables such as density and fluid velocity are obtained by taking suitable moments of the velocity distribution function. The velocity distribution function, na

i

(x,t),

where the superscript i labels the fluid component and the

subscript a indicates the velocity direction, is the number

density of particles at node x, time t with velocity ea , where

a1, . . . ,b. In this study, the particle velocities are represented in terms of a basis set defined on a cubic lattice. The

velocity vectors are directed towards points that are nearest

and next-nearest neighbors about a central lattice site and are

in units of the lattice spacing divided by the time step. In the

literature, this lattice is called a D3Q19 lattice, where the 19

corresponds to the basis set size, b19 23. Time is counted

in discrete time steps, and the fluid particles can collide with

each other as they move under applied forces surface tension, applied shear, etc.. The directions of the particle velocities are discretized in such a way that the isotropy of the

pressure tensor and other properties are preserved 1. The

discretization of the time and spatial coordinates greatly reduces the information required to specify the fluid dynamics,

thus extending the scale of the system that can be modeled

by a given computational resource. This is the computational

virtue of a lattice formulation. The price is that care must be

taken to avoid artifacts that arise from the lattice model.

Macroscopic quantities such as density, i

(x,t), and fluid

velocity, ui

, of each fluid component, i, are obtained by the

following moment sums:

i

x,tmi

ni

x,tmi

a

na

i

x,t1

and

ui

x,t

a

na

i

x,tea

ni

x,t, 2

where mi is the molecular mass of the ith component. While

the distribution function is defined only over a discrete set of

velocities, the actual macroscopic velocity field of the fluid is

nearly continuous.

The time evolution of the particle velocity distribution

function satisfies the following LB equation:

na

i

xea ,t1na

i

x,ta

i

x,t, 3

where a

i is the collision operator representing the rate of

change of the particle distribution due to collisions. The collision operator is greatly simplified by use of the single time

relaxation approximation 23,24,

a

i

x,t1

ina

i

x,tna

i(eq)x,t, 4

where na

i(eq)(x,t) is the equilibrium distribution at (x,t) and

 i is the relaxation time that controls the rate of approach to

equilibrium. The equilibrium distribution can be represented

in the following form for particles of each type 23:

na

i(eq)xtani

x 3

2 1do

i

3ea•v 3

2 3eaea :vvv2,

5

n19

i(eq)xt19ni

x3do3

2 v2, 6

where

v



i

S

mi

a

na

i

ea / i



i

S

mini

x/ i

, 7

and the weights are ta 1

36 for 1	a	12, ta 136 for 13	a

	18, and t19 1

3 . The parameter do

i can be related by selfconsistency to an effective temperature, T, by the following

moment of the equilibrium distribution:

Tx,t



i

mi

a

na

i(eq)x,teav

2

3

i

ni

x,t

. 8

In order that both fluid components have the same temperature, do

i may be defined by the relation doi

12(T/mi),

where we choose units here such that Boltzmann’s constant

kB equals 1.

It has been shown that the above formalism leads to a

velocity field that is a solution of the Navier-Stokes 24

equation with the kinematic viscosity 
 2,5,


c2



i

S

ci i1

2

6 , 9

where ci is the concentration of each component and the

lattice constant c& for this model.

NICOS S. MARTYS AND JACK F. DOUGLAS PHYSICAL REVIEW E 63 031205

031205-2
A. Fluid interaction

An interaction force Fi for each fluid is needed to drive

the phase-separation process. We use a simple interaction,

suggested by Shan and Chen the physical basis of this forcing is given in Appendix A, that depends on the density of

each fluid:

Fi

ni

x

i

S

a

Gii

a nixeaea , 10

where Gii

a 2G, G, and 0 for the cases ea

1, ea2,

and, ii, respectively. G is a coupling constant controlling

the interaction strength. This term is analogous to a nearestneighbor interaction in lattice models of interacting fluids.

The forcing term has been shown to drive the phase separation and to produce an interfacial surface tension effect consistent with the Laplace law 5, which states that there is a

pressure drop proportional to the local curvature at the interface boundary between two fluids.

In the LB model of Shan and Chen, phase separation

takes place when the mutual diffusivity of the binary mixture

becomes negative, providing a condition determining the

critical coupling Gc for phase separation. An analytical expression for the mutual diffusivity has been determined 25.

For a viscosity matched binary mixture in which the particle

masses are also matched ‘‘symmetric fluid mixture’’, phase

separation occurs when the critical coupling equals

Gc T

48

c1n1c2n2c1n1c2n2

28n1n2

n1n2 .

11

It is not ordinarily possible to exactly calculate the critical

coupling for phase separation in three-dimensional liquids,

and this condition for the critical coupling, Gc , is evidently

a clue to the nature of the phase-separation process. Since the

LB method neglects thermal fluctuations that renormalize the

critical coupling constant G, this method is a mean-field

model of fluid mixtures. This observation, which has basic

ramifications for the applicability of the model in comparison with real fluid mixtures, is established numerically in the

next section, where the critical properties are examined to

establish the nature of the model.

Once the forcing is described, it must be properly incorporated into a LB model. Shan-Chen introduced the forcing

by modifying the equilibrium velocity v 2:

ni

xvxni

vx i

Fi

i

, 12

where v is the new velocity used in Eqs. 5and 6. This

approach introduces a momentum transfer between fluids

that preserves momentum globally. The main criticism of the

Shan-Chen model is that when shifting the velocity in the

equilibrium distribution, additional corrections in the pressure tensor will appear that are of order F2.

Instead of shifting the velocity in the equilibrium distribution as in the Shan-Chen model, the forcing between two

fluids can be directly included in the body-force term of the

Boltzmann equation. In the continuum Boltzmann equation,

the body-force term is written Bi

Fi/i•eni(x,e), where

Fi

/i is an acceleration field due to a body force. An expression of this body-force term, to second order in Hermite

polynomials 7, in the discrete velocity space of the D3Q19

lattice is given by

Ba

i

3tani

xeaÀv•

Fi

i

3ea"v ea•

Fi

i  . 13

One can think of this acceleration field as being due to a

‘‘mean field’’ produced by the surrounding molecules see

Appendix A. The main difference between this approach

and the Shan-Chen model is that it avoids terms of order F2

that result from the shift of the velocity in the equilibrium

distribution, so that the linearity of the forcing is preserved.

The effect of this modification of the LB model is investigated below.

Finally, Eq. 10can be modified to mimic an interaction

between the solid surface and fluid 5. Here ni(x¿ea) is

given the value 1 or 0 depending on whether xea resides

on a point in the solid or fluid, respectively, and the value of

Gii

a is then set to allow the solid to attract a fluid wetting

or to repulse a fluid nonwetting.

III. EQUILIBRIUM CRITICAL PHENOMENA

IN A LB FLUID MIXTURE

An understanding of the equilibrium critical phenomena

of fluid mixtures is necessary for modeling the flow of multiphase materials. This understanding is not only required for

estimating phase stability and the type of phase-separation

process droplet growth or bicontinuous fluid pattern formation 12–16, but also for transport properties collective

diffusion, self-diffusion, viscositythat depend sensitively on

the nature of the critical phenomena occurring in the liquid.

This is natural given the existence of order-parameter fluid

compositionfluctuations that cause a mode coupling between momentum and mass transport processes 26–32.

Our first task in developing the LB model for these many

practical applications involving phase-separating fluid mixtures under flow conditionsis to establish the type of equilibrium critical phenomena exhibited by this model fluid

mixture. While the LB models do have hydrodynamic interactions, they do not treat hydrodynamic couplings associated

with non-mean-field contributions to the mode-coupling between compositional and velocity fluctuations 26–29. The

theory is thus a mean-field theory also in the sense of fluid

dynamics, and the implications of this constraint require further investigation.

IV. COEXISTENCE CURVE FOR LB MIXTURE

For a sufficiently large fluid interaction coupling, G, the

LB mixtures of components A and B phase-separate into liquids having coexisting compositions A and B at equilibrium. The composition variables A and B denote the relative volume fractions of the two fluid components. These

CRITICAL PROPERTIES AND PHASE SEPARATION IN . . . PHYSICAL REVIEW E 63 031205

031205-3
dimensionless concentration units are normalized so that

AB1. Compressibility effects on the fluid mixture can

be treated through the addition of an additional vacancy

component c such that ABc1, as in equilibrium

lattice model calculations of compressible mixtures 33, but

this complication is not considered in the present paper.

Increasing the LB coupling G makes the coexisting compositions more enriched in the pure components and thus has

the same qualitative effect as lowering the temperature in

systems exhibiting an ‘‘upper critical temperature’’ type

phase separation i.e., phase separation upon cooling. The

parameter G thus plays a role similar to the binary interaction parameter 
 in the lattice model of fluid mixtures 19–

22,34,



AB
AA
BB/2, 14

where 
AB and 
AA ,
BB denote the mutual and self-nearestneighbor interactions of the fluid components. G is also

analogous to the well depth parameter 
LJ in off-lattice models of phase separation based on Lenard-Jones or related potentials 35. In all these models, it is the relative value of the

‘‘interaction strength’’ (G,
,
LJ) to the temperature that is

the dimensionless coupling constant defining the tendency

toward phase separation. For example, the dimensionless interaction in the lattice model of fluid mixtures is conventionally defined as 19–22,34

q
/kbT, 15

where q and kbT denote the lattice coordination number and

thermal energy, respectively. The inclusion of the q factor is

made to weight the number of possible nearest-neighbor interactions. In magnetic phase transitions, we have the same

form of dimensionless couplings as Eq. 13, where 
 is

replaced by the ‘‘exchange interaction’’ J modeling the

short-range magnetic interparticle interaction 18,36. From

these measures of interaction, we see that lowering the temperature has basically the same effect as increasing the interaction coupling (
,
LJ) for the usual case where ordering

occurs upon cooling. Phase separation in the LB liquid also

occurs when the temperature is lowered with G fixed, and we

similarly define a dimensionless coupling constant,

 GG/kbT. 16

A reduced variable temperature  may then be defined from

the interaction coupling constant,  G ,

 G T/GT/Gc

T/Gc

17

for our simulation performed at fixed temperature T and variable G. For a particular fluid mixture it is natural to fix G and

to vary T so that the reduced temperature variable is defined

as

 G

  T/GTc /G

Tc /G , 18

where Tc is the critical temperature for a fixed value of G.

The absolute value definition in Eqs. 17and 18ensures

that the reduced temperature variable is positive for notational simplicity, but this requires that we must carefully

distinguish between the one-phase and two-phase regions.

All of the computations of the present paper are performed in

the two-phase region.

In Fig. 1, we present our results for the coexistence curve

of a symmetric LB fluid mixture both mass and viscosity

ratios of fluid components are equal. The y axis denotes the

ratio of critical dimensionless coupling to the dimensionless

coupling,  Gc / G , defined in Eq. 16and the x axis denotes

the composition A of the A fluid. We observe that the critical composition c,A of the A component equals c 1

2 a

‘‘symmetric mixture’’, as required by the symmetry of exchange of the fluid components. This exchange symmetry is

well known in lattice models of fluid phase separation 37–

41. The composition difference A

(1)A(2) between

the coexisting phases defines an order parameter for the fluid

phase-separation process. The relation of  to the reduced

temperature is indicative of the type of critical phenomena

‘‘universality class’’under discussion. In a mean-field

model of fluid phase separation,  is described by the general relation 17,18,34

2B , TTc /Tc , TTc, 19

where the order-parameter exponent  and critical amplitude

B for a symmetric incompressible fluid mixture equal 34

B3/2,  1

2 . 20

The dashed line in Fig. 1 is the predicted value. Our data are

consistent with the mean-field prediction as  G→0. Note

that the Shan-Chen model deviates more from the mean-field

prediction than does the simple body-forcing model. In general, it was found that the linear body forcing was somewhat

FIG. 1. Phase diagram of LB fluid mixture. Normalized coupling  Gc / G versus the composition A of fluid A. The solid

circles represent data from the Shan-Chen model and the triangles

represent data from the body-forcing model.  Gc / G also corresponds to the temperature ratio T/Tc .

NICOS S. MARTYS AND JACK F. DOUGLAS PHYSICAL REVIEW E 63 031205

031205-4
more stable. The mean-field theory prediction Eq. 19 is

further examined in Fig. 2, where we plot log10() versus

log10(G) for the lattice data shown in Fig. 1. It is apparent

that a power-law scaling of  on  G is observed over an

appreciable temperature range. The solid line denotes the

prediction of Eq. 19with no free parameters, where  is

equated with  G . Note that the critical temperature is not

adjustable in this comparison, in contrast to most simulations

and experiments in which this quantity is not known exactly.

Of course, the solution of the two-dimensional Ising spin

model and its lattice-gas analog is an exception to this general situation 17,18. Sengers gives an excellent review of

the critical properties of fluids and fluid mixtures that provides much further information about mean-field and nonmean-field critical properties and the ‘‘crossover’’ between

these property scaling regimes 42.

Figures 1 and 2 not only verify that the phase-separation

process in LB fluids is described well by mean-field theory,

but they also establish the utility of our definition of reduced

temperature scale,  G , which is required for other applications involving LB fluid mixtures. For example, we can

quantify the quench depth of our phase-separation measurements by specifying the  G value. These simulations can be

compared to experiments on real fluids at the corresponding

 value. Quantitative agreement with the properties of real

liquids can only be expected for liquids that can be modeled

by mean-field theory over a broad temperature range see the

discussion below. This identification between computational

and real fluids is generally restricted to a temperature range

over which mean-field critical behavior is exhibited to a

good approximation. Strictly speaking, no real fluids are described by mean-field critical behavior, but for many fluids

the approximation should be reasonable provided  is sufficiently far from the critical point defined by the limit 

→0. The Ginzburg criterion defines the temperature range

over which mean-field theory is a reasonable approximation

33,43–50.

Real fluid mixtures are characterized by differences in the

molecular shapes and volumes of the fluid molecules and

asymmetries in the intermolecular interaction potentials that

destroy the symmetry of exchange between the fluid components 38–40. This symmetry breaking is evident in the

shape of the coexistence curve. The graph of  versus 

becomes ‘‘skewed’’ so that the critical composition c no

longer equals 1

2 38–40. The molecular asymmetry effect is

particularly evident in polymer fluid mixtures where the ratio

of the molecular weights and the backbone chain structure

can be adjusted to ‘‘tune’’ the asymmetry of the coexistence

curve 21,22,34,51. The asymmetry becomes extreme in the

case of high molecular weight polymers dissolved in low

molecular weight solvents where A,c of the high molecular

weight component approaches zero with increasing molecular weight 21,52. It is also possible to modify the molecular

weights of a blend to achieve an almost perfect symmetry as

in Fig. 1 53. This symmetry is not usually observed in fluid

mixtures or in single-component fluid phase transitions, although the degree of asymmetry is usually modest in comparison with polymer solutions.

The breaking of the particle exchange symmetry arising

from differences in molecular shape, rigidity, mass, and

other molecular parameters is difficult to describe in a mesoscale fluid model of phase separation. We can obtain a

simple model of this symmetry-breaking phenomenon, however, by considering the idealized Flory-Huggins FHmeanfield model of polymer blend phase separation 21, which

accounts minimally for the molecular mass asymmetry of the

fluid components actually, the model accounts for a volume

asymmetry since this incompressible polymer blend model

assumes all lattice sites are occupied and have equal density.

Notably, the FH model completely ignores polymer topology, monomer asymmetry, polydispersity in the size and

monomer-monomer interactions, and other factors that surely

influence polymer blend stability, but the mass ratio in the

FH model does provide a parameter that allows the asymmetry of the coexistence curve to be ‘‘tuned’’ to fit observations

on real blends. The recently developed lattice cluster meanfield theory generalizes the FH model by incorporating

leading-order correlations associated with chemical connectivity and monomer structure 33.We first consider the

case in which the particle masses are ‘‘asymmetric’’ in the

LB fluid model in the same spirit of approximation. Figure 3

shows the coexistence curve for a LB fluid mixture having a

mass ratio  MMA /MB3, where the concentration difference  between the coexisting phases is given in number

density concentration units rather than the volume fraction

units of Fig. 1. We examine the scaling of  on the quench

depth parameter  G in Fig. 4, where we find a mean-field

scaling exponent 1

2 as in Figs. 1 and 2 and a shift of the

critical coupling to the value Gc0.0135. The asymmetry of

the coexistence curve is quantified by calculating the dependence of the average composition A  ̄(A

(1)A(2))/2 in the

coexistence curve shown in Fig. 3, where A

(1) and A(2) are

compositions of the coexisting phases. According to the

‘‘law of rectilinear diameter’’ of Cailletet and Mathias 54,

 ̄ is linear function of . This linearity is found to a good

approximation in the asymmetric fluid phase-separation coexistence curve shown in Fig. 3. The average composition

A is shown in Fig. 5, where the line denotes the rectilinear

diameter fit,

FIG. 2. Order parameter  versus quench depth parameter

 G .

CRITICAL PROPERTIES AND PHASE SEPARATION IN . . . PHYSICAL REVIEW E 63 031205

031205-5
A  ̄A,cA G , 21

where A is a constant, A0.11, and A  ̄ and A,c are given in

number fraction units. This type of plot is an effective way to

determine the critical composition of an asymmetric fluid

mixture 55. Fluctuation corrections to mean-field theory

can lead to deviations from Eq. 21in real fluids that are

important near the critical point 56.

We can also break the symmetry of interparticle exchange

and thus distort the shape of the phase boundary by varying

the relative volumes of the fluid particles see Appendix B.

In Fig. 6, we show the phase boundary calculated for a range

of values of , the ratio of the volumes of components A and

B, VA /VB . We assume spherically shaped particles so

that  scales as the cube of the ratio of the particle radii, 

(RA /RB)3. By convention, we take the A fluid component

to have the largest molecular volume. Increasing the particle

size asymmetry strongly increases the asymmetry of the

phase boundary, in qualitative agreement with the FloryHuggins theory of polymer phase separation 21,22,34.

Note that the data in Fig. 6 are given in volume fraction

units. The simple Flory mean-field treatment of the FloryHuggins lattice model indicates that c1/11/2, where

 is the relative chain molecular volume see below. The

true critical composition seems to be approximated reasonably well by a similar expression, c1/1, and the

arrows in the figure show the result in comparison with the

data. We note that the phase diagrams of micelle and protein

solutions, where there is a large asymmetry in the size of the

phase-separating species, tend to be asymmetric as in Fig. 6

57.

In applications of the LB model to real measurements, we

can phenomenologically adjust the relative mass and thus

the critical composition in Figs. 3 and 4or relative particle

volume  and identify the LB order-parameter variable A

number densityor volume fraction units, respectively with

the experimentally determined order-parameter concentration unit. While this is generally an approximation, we expect it to provide a reasonable mimic of the critical properties of asymmetric fluid mixtures as in previous experience

with the FH model 58.

A. Interfacial composition profile and correlation length

The interface between phase-separated liquids becomes

diffuse near a critical point where the interfacial tension becomes relatively low. The width of this interface can be

quantified through the determination of the composition interfacial profile A(z), which measures the local composition along a coordinate, z, normal to an interface between the

coexisting phases.

An interface in near-critical fluid mixtures can be probed

by optical reflectivity 59–61 or ellipsometry 62–64 to

determine its width, but accurate measurements of composition gradients across the interface are difficult. Direct measurement of A(z) has recently become possible in thin films

by neutron reflection 65,66, but the broadening of these

profiles by capillary waves and surface wetting effects complicates the interpretation of those measurements so that

it is hard to quantitatively evaluate theory in this important

area 67–70. The LB model allows the determination of

A(z) for an ideal mean-field theory fluid. Some insight into

the fluctuation contribution to A(z) can be obtained by

comparing these calculations to Monte Carlo calculations of

A(z) 69,70. An important property that derives from the

determination of A(z) is the correlation length in the

two-phase region, which governs the average width of the

interface see below. This definition of the correlation

length is more involved in asymmetric fluid mixtures since

the composition profile, A(z), is asymmetric about the center (z0) of the fluid interface.

Figure 7 shows an equilibrium interfacial composition

profile A(z) for a symmetric LB mixture having a quench

depth in the two-phase region,  G0.1. The numerically determined profile A(z) is fit well by the mean-field theory

prediction 11,63,70–72,

Az ̄ /2tanhz/w22

for all  G considered in our paper. The dependence of the

interfacial width w on reduced temperature  G is indicated in

Fig. 8. The mean-field correlation length of the fluid mixture in the two-phase mixture is related to w by 64,70,72

2w 23

so that the determination A(z) affords a means of determining the basic property . A fit to the w data nearest the

critical point in Fig. 8 gives

0.960.05 G

 , 24

where the mean-field value of the critical exponent 
 1

2 is

assumed. Far away from the critical point, where  G

O(1), the correlation length becomes comparable to the

lattice spacing, as in Monte Carlo simulation of phase separation in small molecule liquids 70. The lattice spacing in

the LB model should be interpreted as being comparable to

the average range of the interparticle interaction potential.

This scale is typically comparable to the average molecular

dimensions of the molecule involved 73, and for polymers

this scale can be fairly large 74. Particle clustering can also

increase the magnitude of this scale in small molecule liquids

75.

In a mean-field model of phase separation, the correlation

length  has the same singular dependence on reduced temperature,  G , in the one- and two-phase regions ( and ,

respectively76,

 o

 G1/2 ,  o G1/2 , 25

but the correlation length amplitudes  o

 , oare related by a

constant ratio in mean-field theory 76,

 o

/ o21/2. 26

This ‘‘universal’’ ratio is closer to 2 in real fluid mixtures

exhibiting Ising-type criticality 77,78. The discrepancy between Eq. 26and measurement is illustrative of the large

NICOS S. MARTYS AND JACK F. DOUGLAS PHYSICAL REVIEW E 63 031205

031205-6
property changes that critical fluctuations can induce. Moreover, the critical temperature, in three dimensions, can be

shifted from its mean-field value by as much as 25% by

fluctuations, so that the interaction parameters must be

treated as phenomenological parameters in comparison to experiments in order to ‘‘absorb’’ these discrepancies 34.

The interfacial composition profile A(z) of asymmetric

fluid mixtures composed of mixtures of dissimilar molecules

is also asymmetric 79. We briefly illustrate this effect in

Fig. 9 for a fluid having a mass ratio  M3 and a quench

depth  G0.1. Again, we present our concentration data for

the asymmetric fluid in terms of number density units to

simply realize the effect of fluid asymmetry 80. While the

equilibrium profile may appear to be like the tanh profile

Eq. 22 found in the symmetric case, we could not obtain

a good fit to this function. A good empirical description of

this profile is found by first taking a derivative of the profile

in Figure 9 to find its inflection point and by then fitting to a

tanh profile on each side of the inflection point 65,81. Fig.

9 shows a fit to the data where the characteristic widths wL

5.75 and wR5.0 correspond to the left and right sides of

the inflection point. The average w ̄ (wLwR)/2 provides a

good measure of the average interfacial width. We plan to

discuss the properties of the LB model with a molecular

volume asymmetry further in a separate publication so that

this case is not discussed here. In all the discussion below,

we restrict ourselves to 1.

It has sometimes been reported that two correlation

lengths exist in the two-phase region of fluid mixtures having

asymmetric coexistence curves 75. These measurements

are made by performing light or neutron scattering on coexisting phases in macroscopically phase-separated samples.

The scale of the composition fluctuations appears to occur at

different scales in the A-rich and A-poor coexisting phases

75. Apparently, the measurement process can give rise to

unequal weighting in the different phases to the two sides of

the A(z) profile, leading to different estimates. Although such measurements provide some insight into the

FIG. 3. Phase diagram of an asymmetric LB fluid mixture.

FIG. 4. Order parameter  vs quench depth parameter  G for

an asymmetric LB mixture.

FIG. 5. Rectilinear diameter (AB)/2 versus quench

depth  G .

FIG. 6. Influence of particle size on phase boundary asymmetry.

1.0 filled circles, 4.63 filled squares, and 125.0 open

circles. Dashed lines are included to help guide the eye.

CRITICAL PROPERTIES AND PHASE SEPARATION IN . . . PHYSICAL REVIEW E 63 031205

031205-7
asymmetry of the A(z) profile, they should not be interpreted as implying the existence of two distinct correlation

lengths in the two-phase region.

B. Interfacial tension

Interfacial tension measurements provide a direct means

of probing the interaction between fluids. This property is

crucial in an industrial context for controlling the size and

phase stability of mechanically dispersed droplets and other

transient structures formed in the course of phase separation.

In principle, the interfacial tension, , provides a conceptually simple means of determining the reduced temperature

variable (TTc)/Tc needed to characterize the phase stability of fluid mixtures, but experimental complications

82,83 e.g., high viscosity in polymeric systemshave limited somewhat the application of this method to the critical

phenomena of fluid mixtures. Part of the difficulty is the

need for a more predictive theory of interfacial tension on

which reliable thermodynamic interaction ( or G) measurements can be based. Recently, there has been great effort in

modeling the interfacial tension of polymeric blends by

Monte Carlo simulation methods as a guide to improving

analytic theory in this important area of technological application 69,70.

We calculate the LB interfacial tension  through an integration of the interfacial composition profile,

  Pzz1

2 PxxPyy dz, 27

where Pzz and 1

2 PxxPyy  are the normal and tangential

parts of the pressure tensor, respectively. The numerical values of the interfacial tension for the symmetric LB fluid mixture, shown in Fig. 10, are consistent with a power law,

o G

1.5 , 04.2 28

over a broad temperature range. The exponent 1.5 is an established result for the interfacial tension in mean-field

FIG. 7. Interface composition profile A for a quench depth

 G0.08. z is in the direction normal to the center of the fluid

interface. Solid line is fit to Eq. 22.

FIG. 8. Interfacial width w as a function of the reduced interaction  G .

FIG. 9. The interfacial composition number densityprofile

A(z) of an asymmetric fluid mixture having a mass ratio  M3

and quench depth  G0.08.

FIG. 10. Interfacial tension versus quench depth,  G .

NICOS S. MARTYS AND JACK F. DOUGLAS PHYSICAL REVIEW E 63 031205

031205-8
theory 12,70, and is found to agree reasonably well with

observations outside the critical region 84–86. Fluctuation

effects modify the exponent to a value 1.25 87. This

fluctuation-modified exponent is often found to be quite accurate near the critical point for phase separation 77,78.

The amplitude of the interfacial tension, 0, is a nonuniversal quantity that depends on the interparticle potential range,

interparticle spacing, and molar volume of the liquid. Further

discussion of the origin of interfacial tension in the LB

model is given in Appendix A.

V. COMPARISON OF THE LB FLUID MIXTURE MODEL

TO POLYMER BLENDS

It is apparent from the examination of the LB mixture

model above that the critical properties of this fluid are described by mean-field theory and that the critical exponents

predicted by this model are inconsistent with those measured

for many real mixtures. This limits, of course, the comparison of LB model calculations to certain qualitative trends in

the equilibrium properties of near-critical fluid mixtures.

The theory should become more reliable, however, away

from the critical point.Such inconsistencies can also be expected for certain dynamic properties near the critical point.

For example, the shear viscosity of a near critical Ising-type

fluid mixture diverges near the critical point while no divergence occurs in a model mean-field mixture 88. Modecoupling effects due to compositional flucuations have an

even larger effect on the collective diffusion coefficient

30,31. Although mean-field models of fluid mixtures are

idealized, there is a class of real liquids whose behavior approaches this ideal type critical behavior. The phase separation of polymer blends in the theoretical limit of infinite

molecular weight of the homopolymer components has been

argued 45,46 to be described by mean-field theory, so that

the phase separation of high molecular weight blends should

be reasonably approximated by this idealized model. Monte

Carlo calculations support these theoretical arguments, although the chain length must generally be rather high for this

approximation to apply 89,90. There have also been recent

reports of apparent mean-field critical behavior in fluid mixtures with salts 90 and in ionic fluid mixtures 91.

It is common practice in polymer science to fit the critical

properties of polymer blends to the Flory-Huggins FH

mean-field lattice model of phase separation where all inaccuracies of the modeling monomer structure, chain architecture, compressibility, . . .are absorbed into the phenomenological ‘‘ parameter’’ in Eq. 1692. A virtue of the FH

lattice model is that it often allows the prediction of qualitative trends in the scattering properties of polymer blends. We

can retain this advantage and avoid the conceptual pitfalls of

interpreting  too literally as a ‘‘molecular’’ parameter by

establishing a formal correspondence between the parameters

of this model and the LB fluid mixture model.

In the FH model of polymer blend phase separation, the

reduced temperature variable is 46,92

 FHc/TcT/Tc , 29

where the critical interaction c is defined by the condition

c2NAc 

12NB1c 1 30

and NA and NB are homopolymer polymerization indices.

Symmetric blends are defined by the idealized condition NA

NBN so that c 1

2 and c2/N. The critical composition c of the blend no longer equals 1

2 when the blends are

not symmetric (NANB) 21,46,74,

c NB

1/2

NA

1/2NB1/2 . 31

As mentioned above, the incorporation of asymmetry into

the LB model requires adjusting the mass asymmetry or volume asymmetry to give a variation in the critical composition. We can then mimic the asymmetric phase of the phase

boundary of polymer blends by varying the mass asymmetry

 M or  and formally replacing  G by  FH .

The correlation length  of the FH model in conjunction

with the random-phase approximation 46,74 yields a scaling relation for  in the two-phase region for a symmetric

blend

 0

 FH1/2,  0Rg /6, 32

where Rg is the chain radius of gyration. We see from a

comparison of Eqs. 32to the LB expression Eq. 23that

the lattice spacing in a coarse-grained model of polymer

blends must be large since the lattice spacing is on the order

of Rg . This implies that the lattice spacing must be taken to

depend on chain molecular weight in comparison with measurements. Moreover, the predictions of the LB model must

be considered with caution when physically relevant scales

in physical problems become smaller than this coarsegraining scale lattice spacingof the LB lattice model. This

limitation is natural since the LB model is a mesoscopic

description of a fluid rather than a microscopic model.

There are a number of points to be drawn from our discussion of polymer blend critical properties in comparison

with the LB model of fluid phase separation.

iPolymer blends are reasonable candidates for comparison with the LB mixture model.

iiThe mean-field model gives rise to universal scaling

relationships that should allow fixing the parameters of the

LB model according to the blend molecular characteristics.

This gives some insight into the qualitative variation of the

LB parameters with the variation of molecular structure.

iiiComparison of the LB model with parameters fixed

by the FH model with Monte Carlo calculations of the lattice

model of polymer blends should provide some insight into

the mean-field approximation in the case of properties not

tractable using analytic mean-field theory. For example, we

can compare LB calculations of the interfacial tension to

Monte Carlo calculations for polymer blends that avoid the

mean-field approximation.

ivFixing the LB model parameters through ‘‘matching’’

to the FH model then allows a comparison with dynamical

properties transport coefficientsand processes phase separation, wetting, dewettingof blends calculated using the

mean-field approximation.

CRITICAL PROPERTIES AND PHASE SEPARATION IN . . . PHYSICAL REVIEW E 63 031205

031205-9
vThe expression of the results of measurements in

terms of general and universal scaling relations when they

existoffers advantages to representations involving the phenomenological  interaction parameter. Expressions between

large-scale observable properties deduced from mean-field

theory often have greater applicability than expressions between observables and temperature-like variables such as 

and .

VI. SOME ILLUSTRATIVE CALCULATIONS OF PHASE

SEPARATION WITH AND WITHOUT SHEAR AND

THE INFLUENCE OF INTERACTION BOUNDARIES ON

PHASE SEPARATION

Now that we have established the type of critical phenomena exhibited by the LB model of fluid mixtures and a reduced variable description for some of the basic thermodynamic properties of this model fluid mixture, we can apply

the LB model to the description of phase separation under a

wide range of conditions. In this section, we will illustrate

some phenomena we have investigated in connection with

recent measurements.

The comparison of nonequilibrium phenomena such as

fluid phase separation to LB model calculations requires the

introduction of a dimensionless time unit that is common

between the experimental and computational fluids. For fluid

phase separation, it is conventional to express reduced time

in terms of the mutual diffusion coefficient Dm and the correlation length, 93–95. We thus divide our computational time t by the average initial rate of phase separation,

tps2()2/Dm, deduced from Cahn-Hilliard theory

11,93–95. The mutual diffusion coefficient obtained for the

Shan-Chen model studied in this paper for a viscosity

matched binary mixture is given by

DmT 1G

2n1n2

1Gc1n2c2n1

1

2 , 33

where G12G/T. All times below are reported in the dimensionless time t

 ̄t/tps in our discussion of the LB model

of fluid phase separation. Results are presented from simulations corresponding to a sytem size of 803 in units of lattice

spacing cubed.

A. Phase separation without shear

In Fig. 11, we illustrate the case of a critical composition

(A 1

2 ,mAmB). The value of the ‘‘quench depth’’ equals

 G0.537 and the reduced time, t

 ̄3.6. Periodic boundary

conditions are employed in this LB calculation to minimize

wall effects. Figure 12 shows separation after a later time,

t

 ̄17.8. The pattern is similar in geometric form to Fig. 11,

but the characteristic scale of the pattern is larger after longer

FIG. 11. ‘‘Spinodal’’ phase separation morphology in critical

composition 50-50 relative compositionfluid mixture. The quench

depth equals  G0.537 and t

 ̄3.6.

FIG. 12. Spinodal phase separation in critical composition fluid

mixture. Note the similarity of the structure in Fig. 12 to Fig. 11,

apart from scale. This observation reflects the existence of dynamic

scaling in the mixture coarsening. The quench depth equals  G

0.537 and t

 ̄17.8.

FIG. 13. Off-critical 10-90 relative compositionphase separation showing droplet formation and coarsening. The quench depth

 G0.133 and t

 ̄8.4.

NICOS S. MARTYS AND JACK F. DOUGLAS PHYSICAL REVIEW E 63 031205

031205-10
times, geometrically illustrating the notion of dynamic similarity in the phase-separation coarsening process. A study of

the time dependence of the growth of the phase separation

shows that the pattern scale grows slowly in the early stage

of the phase-separation process, as the local composition

builds up to its coexisting composition value one of the

coexisting phases, B , has been rendered transparent in Fig.

11. At a later stage of phase separation, the pattern scale

growth is roughly linear in time 96,97. Other LB studies

have recently focused on modeling the kinetics of phase

separation, so we do not dwell on this well known phenomFIG. 14. ColorSheared critical composition blend. The

‘‘string’’ structures are observed along shearing planes where the

strings are oriented in the direction of the fluid flow. The quench

depth  G0.537, t

 ̄8.5, and  ̇ t p0.18.

FIG. 15. ColorEarly-stage off-critical 15-85phase separation under shear. For Figs. 15–17, orange represents the regions of

high localized phase fraction of fluid A (A0.15). The green

regions represent the transition to a high localized phase fraction of

fluid B, B0.85. Note the incipient Taylor-Tomitaka instability in

some of the fluid ‘‘strings.’’ The quench depth  G0.287, t

 ̄

16.4, and the dimensionless shear rate  ̇ t p0.56.

FIG. 16. ColorIntermediate stage off-critical 15-85phase

separation (t

 ̄27.6). Droplets form after the fluid threads in Fig.

15 break up by the Taylor-Tomotika instability.

FIG. 17. ColorLate-stage off-critical 15-85phase separation

(t

 ̄56.8). Droplets shown in Fig. 15 reconnect into stringlike

structures that appear to persist and coarsen in cross-sectional dimension with time. The strings in Fig. 17 are contrasted with the

‘‘stringslike’’ structures near the mixture surface in Fig. 14, representative of a critical composition fluid mixture.

CRITICAL PROPERTIES AND PHASE SEPARATION IN . . . PHYSICAL REVIEW E 63 031205

031205-11
FIG. 18. ColorInitial density of three phase systems to undergo surface-driven phase separation. The blue and green regions

correspond to the location of the two-phase mixture. A third phase,

lying above, is rendered as the translucent red region. The quench

depth  G0.7.

FIG. 19. Color‘‘Surface-directed’’ phase separation. Phase

separation of a fluid mixture between interacting solid and air

boundaries. The layered morphology corresponds to ‘‘surfacedirected spinodal decomposition.’’ This stage of the phase separation corresponds to t

 ̄1.0.

FIG. 20. ColorPhase separation of a fluid mixture between

interacting solid and air boundaries. This image indicates the development of an instability that disrupts the layers in late-stage phase

separation. The reduced time equals t

 ̄1.6.

FIG. 21. Color‘‘Disrupted’’ surface-directed phase separation. The disruption effect of ‘‘pinching’’ of the layers allows for

further coarsening and leads to a collapse of the layered structure

shown in Figs. 19 and 20. The reduced time equals t

 ̄4.0.

NICOS S. MARTYS AND JACK F. DOUGLAS PHYSICAL REVIEW E 63 031205

031205-12
enon in the present paper 11,97–100. We next illustrate the

qualitative change in the phase-separation morphology that

occurs under off-critical conditions. Figure 13 shows the

phase-separation morphology for a quench depth,  G

0.133. The rendered composition A is taken to have the

off-critical value A0.1. The phase separation then occurs

through droplet formation rather than the formation of a bicontinuous ‘‘spinodal’’ phase-separation pattern. At a later

time, we observe the droplets to coarsen by coalescence in a

normal manner for off-critical fluids 101. In future work,

we plan to explore the conditions e.g., quench depth, viscosity mismatch, etc.that determine the crossover between

the droplet and bicontinuous phase-separation patterns observed in the early stages of phase separation. Even this basic

aspect of fluid phase separation remains poorly understood

so that materials scientists must rely on engineering correlations 102.

B. Phase separation under a steady shearing flow

We next illustrate a nontrivial application of the LB

method to a situation in which fluid flow is crucially important. Figure 14 shows the phase separation of a critical composition (AB , MAMB) blend for the same quench as

shown in Fig. 11. The upper and lower boundaries in the

figure are energetically ‘‘neutral’’ neither fluid preferentially wets the surface.A hydrodynamic ‘‘stick’’ boundary

condition is imposed at the walls. The top and bottom walls

move at velocities uw and uw such that the dimensionless

shear rate equals  ̇ tps(2uw /d)tps0.18, where d is the

spacing between the walls. The boundary condition is periodic in the direction parallel to the translating planes. We

observe in the top view that the phase-separation pattern appears to have a ‘‘stringlike’’ form at the boundary, but that

the structure within the film is actually more complicated.

The phase separation in the plane perpendicular to the flow is

remarkably undisturbed by the flow and closely resembles a

two-dimensional phase-separation pattern in the absence of

shear. As time proceeds, the ‘‘penetration depth’’ of the

surface-induced ‘‘stringlike’’ structures in Fig. 14 increases

and the phase-separation morphology ultimately coarsens to

a state where the fluid interface lies perpendicular to the

shearing planes and parallel to the flow direction. The

surface-induced ‘‘structuring’’ of the blend morphology

would be accentuated in the early stages of phase separation

if one of the blend components had a preferential attraction

to the shearing boundaries. This surface interaction evidently

has a symmetry-breaking effect, since the fluid interface was

found to lie parallel to the interface of the shearing boundaries at long times, when the boundary interaction was modified in this fashion. This illustrative calculation shows that

the surface interaction can have a large influence on the ultimate alignment of a phase-separated fluid under shear. It is

also apparent that the interpretation of optical and scattering

data on sheared phase-separated fluids is complicated by the

existence of gradients in the composition and structure of the

fluid. Clearly, these observations warrant a thorough investigation of the many parameters that seem relevant to the phenomena quench depth, surface energy, viscosity, molecular

weight mismatch, roughness of shearing surface, steady and

oscillatory shear, etc.. We note that many experiments have

recently reported ‘‘stringlike’’ structures in sheared phaseseparating fluids using both light scattering and optical microscopy techniques 48,103,104, and ‘‘stringlike’’ structures have also been reported in two-dimensional LB phaseseparation simulations 105. The case of two dimensions is

somewhat special, however, since the Taylor-Tomotika

106 instability is suppressed in two dimensions 107.

An illustration of phase-separation in an off-critical blend

under steady shear also provides important insights into the

kinds of phase-separation morphologies that can be expected

experimentally. Figure 15 shows an off-critical blend 15-85

for a dimensionless shear rate,  ̇ tps(2uw /d)tps0.58. Initially, long narrow filaments formed and eventually they

broke apart into droplets due to the well-known TaylorTomotika instability 106. The droplets Fig. 16then became elongated and tilted at approximately 45° relative to

the shear plane, as predicted in the limit of a low concentration of dispersed droplets 108. However, at a later stage of

phase separation Fig. 17, the droplets coalesce to form undulating string structures that seem to persist indefinitely in a

‘‘dynamic string state’’ a video of our simulation can be

found in Ref. 109. The correlated motions of the strings

suggest that the hydrodynamic interactions between the

strings and/or between the strings and the boundary of the

sheared fluid seem to be playing an important role in the

stabilization and formation of the extended string structures.

Subsequent experiments have indicated a similar string formation phenomenon in an off-critical blend sheared at low

shear rates in a parallel plate geometry having a narrow gap

relative to the droplet size 110.

C. Surface-directed phase separation

As a final illustrative application of the LB method, we

consider an example of blend phase separation where one

boundary is solid and the other interface is fluid Figs.

18–21. Figure 18 shows the initial fluid composition at a

quench depth of  G0.7. The dark liquid phase has a preferential interaction 111 with both the solid substrate bottom boundaryand a third fluid ‘‘air’’. This image illustrates the well known phenomena of ‘‘surface-directed

spinodal decomposition,’’ in which the compositional waves

of phase separation are brought into registry with the

symmetry-breaking walls 95,112–116. The coarsening of

the layer structure at early and intermediate times occurs

much like a bulk blend see Fig. 12, but the continued coarsening at long times requires the intermittent loss of fluid

layers. At some point, the undulations within the layers grow

sufficiently large perhaps associated with the rupture of inner layers as required by coarseningto induce perforations

in the outer surface of the blend film at the polymer-air

boundary Fig. 20. This undulation phenomenon then causes

the layered structure to break up into a structure that superficially resembles a spinodal decomposition pattern when

seen from above Fig. 21. A number of studies have indicated the presence of a ‘‘fast mode’’ 117 in layered blend

films, corresponding to a rapidly growing length scale conCRITICAL PROPERTIES AND PHASE SEPARATION IN . . . PHYSICAL REVIEW E 63 031205

031205-13
sistent with a hydrodynamic instability. The instability

shown in Figs. 18–21 provides a possible explanation for the

geometrical nature of this transition. Further computational

and experimental studies of the late-stage coarsening instabilities in layered blends would clearly be interesting to

check this novel picture of phase separation in thin quiescent

films in which the surface exerts a strong perturbing influence on the phase-separating blend film structure.

VII. CONCLUSIONS

The development of the lattice Boltzmann methods of

simulating flows in multiphase liquids has developed rapidly

in recent years. The time has come to evaluate the critical

phenomena that characterize basic thermodynamic and hydrodynamic properties of this type of model. We performed

numerical experiments on a LB fluid model to determine the

equilibrium critical properties that are most important for

comparison with real fluids. The results of those simulations

are represented in a reduced variable description that is

largely independent of the particulars of the model facilitating comparisons with other models of fluid mixtures and

with experiment. This type of representation should also be

advantageous in expressing experimental measurements in a

model-independent form. Our observations indicate that the

critical properties coexistence curves, correlation length, interfacial profile, surface tensionof the LB fluid correspond

to an ideal mean-field fluid over a broad range of temperatures. This makes comparisons of the model to experiment

particularly appropriate to the high molecular weight polymer blends and other fluid mixtures perhaps also including

some ionic fluids and fluid mixtures containing dissolved

salts 90,91that can be reasonably modeled by mean-field

theory.

Now that we have established the equilibrium critical

properties of LB fluid mixtures, we are in a position to study

much more complicated problems involving fluid flow,

phase separation, and interacting complex boundaries. We

illustrated this type of problem in the case of phase separation in critical and off-critical fluid mixtures with and without shear. We also considered the perturbing influence of

boundaries on quiescent phase separation and the flow instabilities that can occur in late stages of phase separation.

Our bulk blend phase-separation studies show that the

morphology of the phase-separation process in its early and

intermediate stages depends on the fluid composition. The

early stage of phase separation corresponds to the growth of

the local composition to the value of the coexisting composition. Coarsening proceeds at a later stage. The phaseseparation morphology had a bicontinuous form under offcritical conditions and the minority phase had a droplet

morphology in a far off-critical blend. Previous LB calculations have emphasized the kinetics of phase separation,

which is not the emphasized phenomena here 11,98–100.

In future work, we plan to study the crossover between the

droplet to bicontinuous phase-separation morphology as a

function of viscosity mismatch, composition, and quench

depth.

Next, we considered the more challenging problem of

phase separation under steady shear. Again we considered

on- and off-critical blends and found the morphologies to be

qualitatively different. Shear had the effect of causing the

phase-separation morphology to ‘‘streak’’ into a stringlike

morphology near the boundary of the critical composition

phase-separating blend, leading to a complex gradient structure within the blend. The ‘‘penetration’’ depth of the

surface-induced strings seemed to grow with time in the

course of phase separation. The ultimate configuration of the

phase-separated blend, alignment parallel or perpendicular to

the flow direction, depends on the polymer surface interaction. These observations of fluid heterogeneity on intermediate time scales clearly raise questions about the proper interpretation of light scattering and optical microscopy studies of

blends under shear, since these methods often involve an

averaging over the gradient structure or are limited to observations of the near-surface properties of the mixture, respectively. The off-critical sheared blend simulations revealed a

tendency toward droplet distortion and tilted alignment with

respect to the shear flow direction. At a latter stage, we observed droplet alignment and the droplets subsequently coalescenced into a stringlike morphology. These strings seem

to be very stable under shear, which we expect to arise from

the strong hydrodynamic interactions between the strings and

the shearing boundaries in these highly confined phaseseparating fluids see Fig. 17. A similar phenomenon in offcritical blends sheared at low rates in a confined geometry

has recently been observed experimentally 110.

In our final illustrative example, we considered the perturbing influence of solid and air boundaries on the phase

separation of a blend. The existence of a ‘‘free’’ deformable

boundary polymer-air interfacemakes this a particularly

instructive example of some of the advantages of the LB

method. We observe the development of composition waves

in the phase-separating blend, as observed in many previous

experimental and simulation studies with a preferential interaction between one of the blend components and the boundaries 110,112–116. The simulation illustrates the process

by which layers are lost in the course of phase coarsening.

These film coarsening processes apparently lead to a destabilization of the layer structure in a late stage of phase separation. The fluctuations within the film associated with successive film rupture processes cause the layer structure to

collapse like a disturbed ‘‘house of cards,’’ leading to a polymer blend morphology superficially resembling a bicontinuous ‘‘spinodal’’ pattern. These observations emphasize the

importance of time-dependent studies of blend film morphologies in measurements on real blend films to properly

interpret their origin.

Our illustrations of LB calculations of blend phase separation were purposely restricted to relatively simple geometries and flows that are under current study for their potential relevance to processing applications. It is also possible to

incorporate many other fluid properties of interest shear dependence of fluid viscosityand other important effects temperature gradients and time-dependent temperature variations, density mismatch of fluid components and segregation

with gravitational and centrifugal fields, fluid wetting and

dewetting on heterogeneous substrates, phase separation of

NICOS S. MARTYS AND JACK F. DOUGLAS PHYSICAL REVIEW E 63 031205

031205-14
blend films on patterned substrates, phase separation in electric fields, phase separation at high rates of flow where inertial effects are important, flow in complex geometries and

with the presence of filler inclusions, etc.. There are many

possibilities for further application. An important challenge

for future theoretical work is the incorporation of fluctuation

effects to better describe fluid properties near the critical

point for phase separation.

ACKNOWLEDGMENTS

The authors would like to acknowledge useful conversations with Hudong Chen and Xiayio He. We thank Newell

Washburn for his constructive remarks related to phase separation in asymmetric fluid mixtures. We also appreciate the

assistance of John Hagedorn and Terry Griffin of the Information and Technology Laboratory at the National Institute

of Standards and Technology in the flow visualization.

N.S.M. would like to acknowledge partial support from the

National Institute of Standards and Technology Center for

Theoretical and Computational Materials Science.

APPENDIX A: CONTRIBUTION TO PRESSURE TENSOR

FROM FLUIDÕFLUID INTERACTION

Here it is useful to discuss the relation between the LB

model and other mean-field models of phase separation.

First, consider the equation for a single-particle distribution

function P1 based on the continuum BBGKY 8,9 formalism, which is extended to the case of multiple species,

 tP1

i

k 1•“P1

i

F •“kP1

i

i, A1

where k is the microscopic momentum, F is the acceleration

due to a body force, and  is a collision operator. It can be

shown, when making a molecular chaos approximation 9,

that the collision operator can be written as

i



j1

s P1

i

r

 1 ,k 1 ,t

mk 1

  d3r2j

r

 2 ,tgi jr 1 ,r 2 ,t

Vi jr12

r

 1

. A2

This approximation of the collision operator is of the form of

a body-force term, F•“kP1. For r

 1r 2d, where d is of

order a few ‘‘effective’’ hard-sphere diameters, gi j(r

 1 ,r 2)

1. After expanding j about r1, the contribution to the

collision operator associated with the attractive intermolecular interaction, i can be approximated by i

j1 s “Vm

i j

•“kP1

i , where Vmi j2ai jj

(r

 1)i j“2j

(r

 1) with ai j

 1

2 dd3rVi j(r) and i j 16 d d3r r2Vi j(r). Vm

i j can be

thought of as a mean-field potential produced by neighboring

particles and “Vm

i j is the associated mean-field force.

The pressure tensor can be determined for this system,

PJa11  1

2

a22  2

2

2a12 12

11 1

2 “1

21“2122 1

2 “2

22“22

12“1•“21“222“21JI 11“1“1

22“2“2212“1“2“2“1.

This expansion is a counterpart of the Cahn-Hilliard or Landau free-energy expansion 12,17. The forcing used in this

paper is for the case in which all terms in the above pressure

tensor are zero except that with the coefficient a12 , which is

proportional to the coupling constant, G, described earlier in

the paper. While there is no explicit inclusion of a surface

tension term in the model studied in this paper, an effective

surface tension force results in the Shan-Chen model due to

how the forcing between fluid components is incorporated.

This can be seen from the leading term in the expansion of F,

Fxxxx 

x , A3

where F2 corrections to the pressure tensor, a feature in the

Shan-Chen model 7, scale with the surface tension as in

standard Cahn-Hilliard models. Higher-order terms, from

Chapman-Enskog analysis, also contribute to the effective

surface tension, along with contributions that arise from

finite-difference approximations of the contiuum equations.

Quantitative predictions of the surface tension require an understanding of all these terms.

APPENDIX B: HARD-SPHERE CORRECTION

TO INTERACTION TERM

The usual lattice Boltzmann method assumes that the fluid

is composed of point particles. To include a volume exclusion interaction and in effect obtain a relative volume of the

fluid particles, we utilize an Enskog hard-sphere model. The

relative volume fraction can be determined from the sphere

radius and the number density. The application of Enskog

theory to multicomponent fluid mixtures is described by Lo ́-

pez de Haro et al. 118 and in references cited in this work.

While there are different formulations of hard-sphere models, such as standard Enskog theory SETand revised Enskog theory RET, we will utilize a form of forcing, arising

from hard-sphere interactions that are treated to lowest order

in density. In this case, the two theories are identical. Further

details are described in Ref. 118. In the isothermal regime,

the additional correction to the forcing due to hard-sphere

collisions, Bi(HS), is

CRITICAL PROPERTIES AND PHASE SEPARATION IN . . . PHYSICAL REVIEW E 63 031205

031205-15
Ba

i(HS)2bi jijc

nj na

i(eq)eav•

nj

x , B1

where ijc is the equilibrium value of the pair correlation

function for spheres of species i and j at contact with the

equilibrium density replaced by the total local equilibrium

density at the point x, bi j 2

3 nji j

3 /, with i j equal to the

distance between sphere centers in contact and  is the local

density. The total forcing on a fluid component i is then Ba

i

Ba

i(HS) , where Bai is defined in Eq. 13.

1 For a review of applications of lattice Boltzmann and lattice

gas methods, see S. Chen and G. Doolen, Annu. Rev. Fluid

Mech. 30, 329 1998; D. H. Rothman and S. Zaleski, Rev.

Mod. Phys. 66, 1417 1994.

2 X. Shan and H. Chen, Phys. Rev. E 47, 1815 1992; 49,

2941 1994.

3 M. R. Swift, W. R. Osborn, and J. M. Yeomans, Phys. Rev.

Lett. 75, 830 1995; W. R. Osborn, E. Orlandini, M. R.

Swift, J. M. Yeomans, and J. R. Banavar, ibid. 75, 4031

1995.

4 A. K. Gunstensen, D. H. Rothman, S. Zaleski, and G. Zanetti,

Phys. Rev. A 43, 4320 1991; A. K. Gunstensen, Ph.D. thesis, MIT, 1992; A. K. Gunstensen and D. H. Rothman, Europhys. Lett. 18, 157 1992.

5 N. S. Martys and H. Chen, Phys. Rev. E 53, 743 1996.

6 F. M. Auzerais, J. Dunsmuir, B. B. Ferre ́ol, N. Martys, J.

Olson, T. S. Ramakrishnan, D. H. Rothman, and L. M.

Schwartz, J. Geophys. Res. B 23, 705 1996.

7 N. S. Martys, X. Shan, and H. Chen, Phys. Rev. E 58, 6855

1998.

8 S. Chapman and T. G. Cowling, The Mathematical Theory of

Non-Uniform Gases, 3rd ed. Cambridge University Press,

London, 1970. BBGKY is an abbreviation for Bogoliubov,

Born, Green, Kirkwood, and Yvon.

9 R. L. Liboff, Kinetic Theory, 2nd ed. John Wiley & Sons,

New York, 1998.

10 N. S. Martys, Int. J. Mod. Phys. C 10, 1367 1999.

11 J. W. Cahn and J. E. Hilliard, J. Chem. Phys. 28, 258 1958;

31, 688 1959.

12 J. W. Cahn, Trans. Metall. Soc. AIME 242, 166 1968.

13 A. J. Bray, Adv. Phys. 43, 357 1994.

14 J. D. Gunton, M. Miguel, and P. S. Sahni, Phase Transitions

and Critical Phenomena, edited by C. Domb and J. L. Lebowicz Academic Press, New York, 1983, Vol. 9.

15 S. C. Glotzer, Annu. Rev. Comput. Phys. 2, 1 1995.

16 E. Orlandini, M. R. Swift, and J. M. Yeomans, Europhys.

Lett. 32, 463 1995.

17 H. E. Stanley, Introduction to Phase Transitions and Critical

Phenomena Oxford University Press, New York, 1971.

18 R. J. Baxter, Exactly Solved Models in Statistical Mechanics

Academic Press, New York, 1982.

19 A. R. Miller, The Theory of Solutions of High Polymers

Clarendon Press, Oxford, 1948.

20 J. H. Hildebrand and R. L. Scott, The Solubility of NonElectrolytes Reinhold, New York, 1950.

21 P. J. Flory, Principles of Polymer Chemistry Cornell University Press, Ithaca, New York, 1953; P. G. de Gennes, Scaling

Concepts in Polymer Science Cornell University Press, Ithaca, 1979.

22 K. F. Freed and J. Dudowicz, Theor. Chim. Acta 82, 357

1992.

23 Y. H. Qian, D. d’Humie`res, and P. Lallemand, Europhys.

Lett. 17, 479 1992.

24 H. Chen, S. Y. Chen, and W. H. Matthaeus, Phys. Rev. A 45,

R5339 1992.

25 X. Shan and G. Doolen, Phys. Rev. E 54, 3614 1996.

26 R. Perl and R. A. Ferrel, Phys. Rev. A 6, 2358 1972; R. A.

Ferrel, Phys. Rev. Lett. 24, 1169 1970; see also J. F. Douglas, Macromolecules 25, 1468 1992.

27 K. Kawasaki, Ann. Phys. N.Y.61, 1 1970.

28 O. Onuki and K. Kawasaki, Ann. Phys. N.Y.121, 456

1979.

29 O. Onuki, K. Yamazaki, and K. Kawasaki, Ann. Phys. N.Y.

131, 217 1981.

30 J. V. Sengers, Int. J. Thermophys. 6, 203 1985; H. C.

Burstyn and J. V. Sengers, Phys. Rev. Lett. 45, 259 1980;

Phys. Rev. A 25, 448 1982.

31 P. Stepanek, T. P. Lodge, C. Kedrowski, and F. S. Bater, J.

Chem. Phys. 94, 8289 1991.

32 J. C. Nieuwoudt and J. V. Sengers, J. Chem. Phys. 90, 457

1989.

33 J. Dudowicz, M. Lifschitz, K. F. Freed, and J. F. Douglas, J.

Chem. Phys. 99, 4804 1993; M. Lifschitz, J. Dudowicz, and

K. F. Freed, ibid. 100, 3957 1994.

34 A. Sariban and K. Binder, J. Chem. Phys. 86, 5859 1987.

35 D. A. McQuarrie, Statistical Mechanics Harper and Row,

New York, 1976.

36 D. C. Mattis, Theory of Magnetism Harper and Row, New

York, 1985.

37 T. Lee and C. Yang, Phys. Rev. 87, 410 1952.

38 N. D. Mermin and J. J. Rehr, Phys. Rev. A 4, 2408 1971; N.

D. Mermin, Phys. Rev. Lett. 26, 169 1971.

39 B. Widom and J. Rowlinson, J. Chem. Phys. 52, 1670 1970.

40 J. A. Zollweg and G. W. Mulholland, J. Chem. Phys. 57,

1021 1972.

41 G. Stell and P. C. Hemmer, J. Chem. Phys. 56, 4274 1972.

42 J. V. Sengers, in Supercritical Fluids: Fundamentals for Applications, edited by E. Kiran and J. M. H. Levelt Sengers

Kluwer, Dordrecht, 1994, pp. 231–271.

43 V. L. Ginzburg, Fiz. Tverd. Tela Leningrad2, 2031 1960

Sov. Phys. Solid State 2, 1824 1960.

44 M. A. Anisimov, S. B. Kiselev, J. V. Sengers, and S. Tan,

Physica A 188, 487 1992.

45 P. G. de Gennes, J. Phys. FranceLett. 38, L441 1977.

46 K. Binder, J. Chem. Phys. 79, 6378 1983; Phys. Rev. A 29,

341 1984.

47 F. S. Bates, J. Rosedale, P. Stepanek, T. P. Lodge, P. Wiltzius, G. H. Fredrickson, and R. P. Hjelm, Phys. Rev. Lett. 65,

1893 1990.

NICOS S. MARTYS AND JACK F. DOUGLAS PHYSICAL REVIEW E 63 031205

031205-16
48 E. K. Hobbie, L. Reed, C. C. Huang, and C. C. Han, Phys.

Rev. E 48, 1579 1993.

49 D. Schwahn, T. Schmackers, and M. Mortensen, Phys. Rev.

E 52, R1288 1995.

50 J. H. M. Levelt Sengers and J.V. Sengers, Perspectives in

Statistical Physics, edited by H. J. Raveche ́ North Holland

Publ. Co., New York, 1981.

51 S. V. Kazakov and N. I. Chernnova, Russ. J. Phys. Chem. 69,

1010 1995.

52 A. A. Povodyrev, M. A. Anisimov, and J. V. Sengers,

Physica A 264, 345 1999.

53 I. C. Sanchez, J. Phys. Chem. 93, 6983 1989; P. Damay and

F. Leclercq, J. Chem. Phys. 95, 590 1991.

54 L. Cailletet and E. C. Mathias, C. R. Hebd. Seances Acad.

Sci. 102, 1202 1886; 104, 1563 1887; E. A. Guggenheim,

J. Chem. Phys. 13, 253 1945.

55 J. W. Yu, J. F. Douglas, E. K. Hobbie, S. Kim, and C. C.

Han, Phys. Rev. Lett. 78, 2664 1997.

56 M. W. Pestak, R. R. Goldstein, M. H. W. Chan, J. R. de

Bruyn, D. A. Balzarini, and N. W. Ashcroft, Phys. Rev. B 35,

599 1987; J. S. Rowlinson, Nature London319, 362

1985.

57 D. Blankschtein, G. M. Thurston, and G. B. Benedek, Phys.

Rev. Lett. 54, 995 1985; X.-H. Guo and S.-H. Chen, ibid.

64, 1979 1990.

58 The dissimilarity in molecular size in a fluid mixture can be

modeled in a more physically realistic way by including different size parameters for the interacting species in the specification of the collision operator in an Enskog-type model of

the particle interaction. See Appendix B and X. He, H. Shan,

and G. Doolen, Phys. Rev. E 57, R13 1998.

59 G. H. Gilmer, W. Gilmore, J. Huang, and W. W. Web, Phys.

Rev. Lett. 14, 491 1965.

60 D. Beaglehole, Physica B & C 112, 330 1982.

61 D. Beysens and M. Robert, J. Chem. Phys. 87, 3056 1987.

62 D. G. Miles, Jr. and J. W. Schmidt, J. Chem. Phys. 92, 3881

1990.

63 J. W. Schmidt, Physica A 172, 40 1991.

64 J. W. Schmidt and M. R. Moldover, J. Chem. Phys. 99, 582

1993.

65 A. Karim, G. P. Felcher, and T. P. Russel, Macromolecules

27, 6973 1994.

66 T. Karle, J. Klein, and K. Binder, Phys. Rev. Lett. 77, 1318

1998.

67 F. P. Buff, R. A. Lovett, and F. H. Stillinger, Phys. Rev. Lett.

15, 621 1965.

68 J. V. Sengers, J. M. J. Van Leeuwen, and J. W. Schmidt,

Physica A 172, 20 1991.

69 A. Werner, F. Schmid, M. Mu ̈ller, and K. Binder, J. Chem.

Phys. 107, 8175 1997.

70 M. Mu ̈ller, K. Binder, and W. Oed, J. Chem. Soc., Faraday

Trans. 91, 2369 1995.

71 J. D. van der Waals, Z. Phys. Chem., Stoechiom. Verwandtschaftsl. 13, 657 1894.

72 S. Fisk and B. Widom, J. Chem. Phys. 50, 3219 1969.

73 P. Debye, J. Chem. Phys. 31, 680 1959.

74 K. Binder, Adv. Polym. Sci. 112, 181 1994.

75 P. Chieux, J.-F. Jal, L. Hily, J. Dupuy, P. Leclercq, and P.

Damay, J. Phys. IV C-5, 3 1991.

76 V. G. Vaks, A. I. Larkin, and S. A. Pikin, Zh. E ́ ksp. Teor.

Fiz. 51, 361 1966Sov. Phys. JETP 24, 240 1967.

77 M. R. Moldover, Phys. Rev. A 31, 1022 1985.

78 H. Chaar, M. R. Moldover, and J. W. Schmidt, J. Chem.

Phys. 85, 418 1986; B. Widom, J. Stat. Phys. 52, 1343

1988; I. Szleifer and B. Widom, J. Chem. Phys. 90, 7524

1989.

79 I. Szleifer, J. Chem. Phys. 92, 6940 1990.

80 The coexistence curve for  M1 is actually fairly symmetric

in volume fraction units.

81 A. Karim, Ph.D. thesis, Northwestern University, 1990 unpublished.

82 R. L. Sammler, R. P. Dion, C. J. Carriere, and A. Cohen,

Rheol. Acta 31, 554 1992.

83 S. H. Anastasiadis, I. Gencarz, and J. T. Koberstin, Macromolecules 21, 2980 1988.

84 D. Atack and O. K. Rice, Discuss. Faraday Soc. 15, 210

1953. See Ref. 1 for a discussion of these data.

85 A. Stammer and B. A. Wolf, Polymer 39, 2065 1998.

86 M. Heinrich and B. A. Wolf, Polymer 33, 1927 1992.

87 B. Widom, J. Chem. Phys. 43, 3892 1965.

88 J. F. Douglas, Macromolecules 25, 1468 1992.

89 H. P. Deutsch and K. Binder, J. Phys. II 3, 1049 1993; M.

Mu ̈ller and K. Binder, Macromolecules 28, 1825 1995.

90 J. Jacob, A. Kumar, M. A. Anisimov, A. A. Povodyrev, and

J. V. Sengers, Phys. Rev. E 58, 2188 1998; B. M. Jaffer Ali

and A. Kumar, Phys. Lett. A 237, 257 1998.

91 K. S. Pitzer, M. Conceicao, P. de Lima, and D. R. Schreiber,

J. Chem. Phys. 89, 1854 1985; D. R. Schreiber, M. Conceicao, P. de Lima, and K. S. Pitzer, ibid. 91, 4087 1987; R. R.

Singh and K. S. Pitzer, ibid. 92, 6775 1990.

92 A. Sariban, K. Binder, and D. W. Heerman, Phys. Rev. B 35,

6873 1987.

93 Y. C. Chou and W. I. Golburg, Phys. Rev. A 23, 2105

1979.

94 M. Takenaka, T. Izumitani, and T. Hashimoto, Macromolecules 20, 2257 1987.

95 B. P. Lee, J. F. Douglas, and S. C. Glotzer, Phys. Rev. E 60,

5812 1999.

96 L. P. McMaster, Adv. Chem. Ser. 142, 43 1975; E. D. Siggia, Phys. Rev. A 20, 595 1970.

97 S. Chen and T. Lookman, J. Stat. Phys. 81, 223 1995.

98 M. R. Swift, E. Orlandini, W. R. Osborn, and J. M. Yeomans,

Phys. Rev. E 54, 5041 1996.

99 J. Olson and D. H. Rothman, J. Stat. Phys. 81, 199 1995.

100 F. J. Alexander, S. Chen, and D. W. Grunau, Phys. Rev. B

48, 634 1993.

101 F. Perrot, P. Guenoun, T. Baumberg, D. Beysens, Y. Garrabos, and B. LeNeindre, Phys. Rev. Lett. 73, 688 1994.

102 I. S. Miles and A. Zurek, Polym. Eng. Sci. 28, 1796 1988;

G. M. Jordhamoo, J. A. Manson, and L. H. Sperling, ibid. 26,

517 1986.

103 T. Hashimoto, K. Matsuzaka, E. Moses, and A. Onuku, Phys.

Rev. Lett. 74, 126 1995; E. Moses, T. Kume, and T. Hashimoto, ibid. 72, 2037 1994; E. K. Hobbie, S. Kim, and C. C.

Han, Phys. Rev. E 54, 5909 1996; S. Kim, E. K. Hobbie, J.

Yu, and C. C. Han, Macromolecules 30, 8245 1997.

104 A. H. Krall, J. V. Sengers, and K. Hamano, Phys. Rev. Lett.

69, 1963 1992.

CRITICAL PROPERTIES AND PHASE SEPARATION IN . . . PHYSICAL REVIEW E 63 031205

031205-17
105 A. J. Wagner and J. M. Yeomans, Phys. Rev. E 59, 4366

1999.

106 G. I. Taylor, Proc. R. Soc. London, Ser. A 138, 41 1932; S.

Tomotika, ibid. 150, 322 1935.

107 M. San Miguel, M. Grant, and J. D. Gunton, Phys. Rev. A 31,

1001 1985; L. Sung, A. Karim, J. F. Douglas, and C. C.

Han, Phys. Rev. Lett. 76, 4368 1996.

108 J. M. Rallison, J. Fluid Mech. 209, 465 1981; H. A. Stone,

Annu. Rev. Fluid Mech. 26, 65 1994.

109 http://www.itl.nist.gov/div895/savg/example/index.html.

110 K. B. Migler, Phys. Rev. Lett. to be published.

111 The incorporation of component-specific surface interactions

in multicomponent fluids requires a modification to the equilibrium distribution to maintain its positive-definite property.

This technical point will be discussed in a separate publication.

112 R. A. L. Jones, L. J. Norton, E. J. Kramer, F. S. Bates, and P.

Wiltzius, Phys. Rev. Lett. 66, 1326 1991.

113 G. Krausch, C. A. Dai, E. J. Kramer, and F. Bates, Phys. Rev.

Lett. 71, 3669 1993.

114 R. C. Ball and R. L. H. Essery, J. Phys.: Condens. Matter 2,

10 303 1990.

115 S. Puri and K. Binder, Phys. Rev. A 46, R4487 1992; Phys.

Rev. E 49, 5359 1993.

116 J. F. Marko, Phys. Rev. E 48, 2861 1993.

117 P. Wiltzius and A. Cumming, Phys. Rev. Lett. 66, 3000

1991; C. Harrison, W. Rippard, and A. Cumming, Phys.

Rev. E 52, 723 1995.

118 M. Lo ́pez de Haro, E. G. D. Cohen, and J. M. Kincaid, J.

Chem. Phys. 78, 2746 1983.

NICOS S. MARTYS AND JACK F. DOUGLAS PHYSICAL REVIEW E 63 031205

031205-18