DTIC ADA639827: Mass-Conserved Phase Field Models for Binary Fluids

6.9 MB·English

by  Defense Technical Information Center

Report Documentation Page Form Approved
OMB No. 0704-0188
Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing the collection of information Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington
VA 22202-4302 Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it

does not display a currently valid OMB control number
1. REPORT DATE 3. DATES COVERED
2011 2. REPORT TYPE 00-00-2011 to 00-00-2011
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
Mass-Conserved Phase Field Models for Binary Fluids
5b. GRANT NUMBER
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S) 5d. PROJECT NUMBER
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION
University of South Carolina,Interdisciplinary Mathematics REPORT NUMBER
Institute,Columbia,SC,29208
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)
11. SPONSOR/MONITOR’S REPORT
NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT
Approved for public release; distribution unlimited
13. SUPPLEMENTARY NOTES
14. ABSTRACT
The commonly used incompressible phase field models for non-reactive, binary fluids, in which the
Cahn-Hilliard equation is used for the transport of phase variables, conserve the total volume of each
phase as well as the material volume, but do not conserve the mass of the fluid mixture when the densities
of two components are different. In this paper, we formulate the phase field theory for mixtures of two
incompressible fluids, consistent with the quasi-compressible theory [28], to ensure the conservation of
mass and momentum for the fluid mixture as well as the volume for each fluid phase. In this formulation,
the mass-average velocity is no longer divergence-free (solenoidal) when the densities of two components in
the mixture are not equal, making it a compressible model subject to an internal constraint. An efficient
numerical method is then devised to enforce this compressible internal constraint. Numerical simulations
in confined geometries for both the compressible and the incompressible models are carried out using
spatially high order spectral methods to contrast the model predictions. Numerical comparisons show that
(a) the predictions by the two models agree qualitatively in the situation where the interfacial mixing layer
is thin; and (b) the predictions differ significantly in binary fluid mixtures undergoing mixing with a large
mixing zone. The numerical study delineates the limitation of the commonly used incompressible phase
field model and thereby cautions its predictive value in simulating well-mixed binary fluids.
15. SUBJECT TERMS
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF
ABSTRACT OF PAGES RESPONSIBLE PERSON
a REPORT b ABSTRACT c THIS PAGE Same as 28
unclassified unclassified unclassified Report (SAR)
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
Mass-Conserved Phase Field Models for Binary Fluids
JieShen∗, Xiaofeng Yang†and QiWang‡
Abstract
The commonly used incompressible phase field models for non-reactive, binary fluids, in which
the Cahn-Hilliard equation is used for the transport of phase variables, conserve the total volume of
each phase as well as the material volume, but do not conserve the mass of the fluid mixture when
the densities of two components are different. In this paper, we formulate the phase field theory for
mixtures of two incompressible fluids, consistent with the quasi-compressible theory [28], to ensure
the conservation of mass and momentum for the fluid mixture as well as the volume for each fluid
phase. In this formulation, the mass-average velocity is no longer divergence-free(solenoidal) when
the densitiesof two componentsin the mixtureare notequal, makingit a compressiblemodelsubject
to an internal constraint. An efficient numerical method is then devised to enforce this compressible
internal constraint. Numerical simulations in confined geometries for both the compressible and the
incompressiblemodelsarecarriedoutusingspatiallyhighorderspectralmethodstocontrastthemodel
predictions. Numericalcomparisonsshowthat(a)thepredictionsbythetwomodelsagreequalitatively
inthesituationwheretheinterfacialmixinglayeristhin; and(b)thepredictionsdiffersignificantlyin
binaryfluidmixturesundergoingmixingwithalargemixingzone. Thenumericalstudydelineatesthe
limitation of the commonlyused incompressiblephase field modeland therebycautionsits predictive
valueinsimulatingwell-mixedbinaryfluids.
1 Introduction
Phase field models have been used successfully to study a variety of interfacial phenomena like equi-
librium shapes of vesicle membranes [12, 13, 14, 15, 16, 35], blends of polymeric liquids [36, 37, 38, 17],
multiphase fluid flows [19, 23, 24, 28, 25, 41, 40, 42, 43, 44, 45], dentritic growth in solidification, mi-
crostructure evolution [21, 29, 22], grain growth [8], crack propagation [9], morphological pattern forma-
tion in thin films and on surfaces [26, 30], self-assembly dynamics of two-phase monolayer on an elastic
substrate [27], a wide variety of diffusive and diffusion-less solid-state phase transitions [10, 39], dislo-
cation modeling in microstructure, electro-migration and multiscale modeling [34]. Multiple phase-field
methods can be devised to study multiphase materials [40]. Recently, phase field models are applied to
study liquid crystal drop deformation in another fluid, liquid films, polymer nanocomposites, and biofilms
[19,23,24,28,25,41,40,42,43,44,18,46,5].
Comparing to other mathematical and computational technologies available for studying multi-phase
materials,thephase-fieldapproachexhibitsaclearadvantageinitssimplicityinmodelformulation, easeof
numerical implementation, and the ability to explore essential interfacial physics at the interfacial regions
etc. Computingtheinterfacewithoutexplicitlytrackingtheinterfaceisthemostattractivenumericalfeature
∗
[e-mail],DepartmentofMathematics,PurdueUniversity,WestLafayette,IN46907,USA
†[e-mail], Department of Mathematics and NanoCenter at USC, University of South Carolina, Columbia, SC
29028,USA
‡[e-mail], Department of Mathematics and NanoCenter at USC, University of South Carolina, Columbia, SC
29028, USA;SchoolofMathematics,NankaiUniversity,Tianjin,P.R.China,300071; BeijingComputationalScienceResearch
Center,Beijing,P.R.China,100084
1
of this modeling and computational technology. Since the pioneering work of Cahn and Hilliard in the
50’s of the last century, the Cahn-Hilliard equation has been the foundation for various phase field models
[6, 7]. It arises naturally as a model for multiphase fluid mixtures should the entropic and mixing energy
of the mixture system be known. For immiscible binary fluid mixtures, one commonly uses a labeling or
a phase variable (also known as a volume fraction or an order parameter) φ to distinguish between distinct
fluid phases. For instance φ=1 indicates one fluid phase while φ=0 denotes the other fluid phase in an
immiscible binary mixture. The interfacial region is tracked by 0< φ < 1. Given the historical reason,
most mixing energies are calculated in terms of the volume fraction instead of the mass fraction in the
literature[20,11]. Consequently, thesystemfreeenergyincluding theentropicandmixingcontribution has
been formulated inthevolume fraction aswell[20,11]. Wedenote thesystem free energy forthematerial
system to be modeled by F(φ,∇φ,···). A transport equation for the volume fraction φ along with the
conservation equation ofmomentum and continuity equation constitutes theessential partofthegoverning
system of equations for the binary fluid mixture. The volume fraction serves as an interval variable for the
fluidmixture.
Inthe literature onimmiscible binary mixtures ofincompressible fluids, oneuses theconcept ofchem-
ical potential to formulate the transport equation for the volume fractions of the fluids φ and φ . In this
1 2
formulation, thematerialincompressibility isontheonehandmodeledbythecontinuity equation
∇·v=0, (1.1)
while on the other hand, interpreted as the invariant property of the sum of the volume fractions for the
two fluid components, i.e., φ +φ =1 if we denote φ=φ and φ =1−φ. This assumption is plausible
1 2 1 2
and indeed consistent withthefluidcompressibility (1.1) only ifthe twocomponents areeither completely
separated by phase boundaries when their densities are not equal or possibly mixed when the densities are
identical. Otherwise, there is a potential inconsistency with the usual conservation of mass. This inconsis-
tencyhasbeenidentifiedin[28],butignoredbymanypractitionersusingphasefieldmodelingtechnologies.
Wenotethatthisinconsistency occurs onlyinthemixedregionofthetwoincompressible fluids,wherethe
incompressible condition(1.1)isnolongervalid,indicatingthemixtureisnolongerincompressible despite
that each fluid component participating in mixing is incompressible. This type of fluids is referred to as
quasi-compressible in[28].
Thispaper aimsatdiscussing theinconsistency forbinary mixtures oftwoincompressible fluidsofun-
matched densities andviscosities andproviding aquantitative assessment forthequasi-compressible phase
field model that obeys the conservation of both mass and volume against the incompressible one that only
respects the volume conservation. The paper is organized as follows. First we discuss the mathematical
formulation of the phase field theory for binary viscous fluid mixtures and its various approximations and
their ramifications. Secondly, we develop a new set of numerical algorithms, which enforce the mass con-
servation,tosolvethegoverningsystemoffluidflowequations. Thirdly,weimplementthealgorithmsusing
spatiallyhighorderspectralmethodsanddiscussthediscrepanciesbetweentheadhocincompressiblephase
fieldmodelandthequasi-compressible phasefieldmodelintworepresentative examples.
2 The mathematical model
We revisit the derivation of the governing system of equations for a binary mixture of incompressible
viscous fluids, which includes the transport equation for a phase variable (the volume fraction) and the
conservation equationsformassandlinearmomentum. Theconservation equations forabinarysystemcan
be formulated in two different ways: either as a two fluid model or a one fluid two component model [1,
4,2]. Fornontrivial fluidsimulations, theonefluidmulti-component formulation often yields aconvenient
governingsystemofequationsandeasy-to-implementboundaryconditionsforthemodel’shydro-dynamical
2
variables. The phase field theory falls naturally into the one fluid two component formulation [6]. In the
phase field formulation, chemical reaction between the two distinct components can take place so that one
component can be turned into the other component. However, the overall mass must be conserved. In
this paper, we will not address the phase field formulation with chemical reactions. This topic deserves a
separate discussion ofitsown.
2.1 Governing equations
Inaphasefieldtheory, thetransport equationforthevolumefractionofonefluidphaseisgivenby
φ +∇·(φv)=∇·(λ∇μ), (2.1)
t
where v is an average velocity to be clarified below, λ=λ(φ) is the mobility function, and μ= δF is the
δφ
chemical potential of the material system. The mobility function λ is often taken as a constant λ , but is
0
preferably afunctionofφintheform:
λ=λ φ(1−φ). (2.2)
0
The Cahn-Hilliard equation with the volume fraction dependent mobility is called singular or modified
Cahn-Hilliard equation. Often, it is approximated simply by a constant value λ = λ in studying phase
0
separated, immisciblefluids. Theresultant equation isthewell-knownCahn-Hilliard equation.
Thefreeenergy ofthe mixturesystem isnormally afunction ofthelabeling function ofphase function
anditshigherorderderivatives (onlythefirstorderisincluded hereforbrevity):
F =F(φ,∇φ). (2.3)
In this paper, we consider the mixture of two incompressible fluids with constant mass density ρ and ρ ,
1 2
respectively. Thetotaldensity ofthemixtureisthengivenby
ρ=ρ φ+ρ (1−φ). (2.4)
1 2
Weidentify vas the mass-average velocity for the mixture. Then, the conservation equations for mass and
momentumaregivenby
ρ +∇·(ρv)=0,
t
(2.5)
(ρv) +∇·(ρvv)=∇·(τ)−φ∇μ+F ,
t e
where F is the external force and φ∇μ is the ”elastic force” or the ”surface force” due to the interfacial
e
energy f(φ). The surface force −φ∇μ can be replaced by μ∇φ modulo a surface term which is normally
zero. Inlightofthetransport equation forthevolumefraction, wehave
ρ∇·v=−(ρ −ρ )(φ +v·∇φ)=−(ρ −ρ )(∇·(λ∇μ)−φ∇·v), (2.6)
1 2 t 1 2
wethenderivefromtheaboveand(2.1)that
− −
ρ ρ ρ ρ
∇·v= 2 1[φ +∇·(φv)]= 2 1[∇·(λ∇μ)]. (2.7)
t
ρ ρ
2 2
Itisapparentthatthedivergencefreeconditionforthemass-averagevelocityfieldissatisfiedonlyifρ =ρ
1 2
or ∇·(λ∇μ)=0. Otherwise, the mass conservation equation serves as a constraint for the velocity field,
which determines the undetermined pressure in the constraint hydrodynamic theory for fluidmixtures. We
3
note that ∇·(λ∇μ)is normally not zero for aspatially inhomogeneous system. Hence, as long as ρ (cid:3)=ρ ,
1 2
eq. (2.7)servesasaconstraint.
To close the system of equations, we must come up with aconstitutive equation for the stress tensor τ.
Weconsider themixturemadeupofviscousfluids. Forviscousfluids,thestressconstitutive equation is
τ=τ +2ηD+νtr(D)I, (2.8)
c
where τ is the constraint stress responsible to maintain the constraint eq. (2.7) without any contribution
c
to the entropy production, η is the shear viscosity, ν is the volumetric viscosity, and D is the rate of strain
tensor. The ratio between ν and η depends on the property of the material and is roughly 4.3 for water for
example. The viscosity coefficients for the fluid mixture are interpolated through the volume fraction and
givenby
η=η φ+η (1−φ),ν=ν φ+ν (1−φ), (2.9)
1 2 1 2
whereη1,2,ν1,2 areconstant shearandvolumetric viscosities forfluid1andfluid2,respectively.
Todealwithconstraint (2.7),weaugmentthefreeenergywithatermF¯calledtheconstraint response:
F =F+F¯. (2.10)
Based on the second law of thermodynamics in the form of the Clausius-Duhem inequality, the constraint
response doesnotcontribute toentropy production, i.e.,
δF¯
τ :D− φ˙ =0, (2.11)
c
δφ
whereφ˙ = ∂φ+∇·(vφ)isthematerialderivative. Werewriteeq. (2.7)as
∂t
−
ρ ρ
I:D+ 1 2φ˙ =0. (2.12)
ρ
2
Foreq. (2.11)mustbevalidforallthermodynamic processes thatobeys(2.12),wededucethat
τ =−pI (2.13)
c
where p is the hydrodynamic pressure, and the free energy component (F¯) corresponding to the constraint
response:
−
ρ ρ
F¯= 1 2φp. (2.14)
ρ
2
Ifwechooseμ= δF ,weobtainasetofequationsthatrespecttheconservationofmassandtotalvolume:
δφ
φ +∇·(φv)=∇·(λ∇μ), (2.15a)
t
(ρv) +∇·(ρvv)=ρ[v +v·∇v]=∇·(2ηD+νtr(D)I)−∇p−φ∇μ+F
t t e
=∇·(η∇v)+∇((η+ν)∇·v)−∇p−φ∇μ+F , (2.15b)
e
−
ρ ρ
∇·v= 2 1[∇·(λ∇μ)], (2.15c)
ρ
2
δF δF ρ −ρ
μ= = + 1 2p. (2.15d)
δφ δφ ρ
2
4
Werefer(2.15)asthecompressiblemodel2inthispaper. Ontheotherhand,ifwesetμ= δF inthesystem
δφ
of equations listed above, weobtain another set of equations, which werefer as the compressible model 1.
Withthehelpof(2.15c),thetransport equation forφcanberecastinto
ρ
φ +∇·(φv)= 2 ∇·v (2.16)
t −
ρ ρ
2 1
provided ρ (cid:3)=ρ .
1 2
The above compressible models preserve the mass conservation and are compressible inside the mix-
ing/interfacial region. On the other hand, the incompressible model, in which the mass average velocity
fieldisassumedsolenoidal, consists ofthefollowingequations:
φ +∇·(φv)=∇·(λ∇μ), (2.17a)
t
ρ[v +v·∇v]=∇·(η∇v))−∇p−φ∇μ+F , (2.17b)
t e
∇·v=0, (2.17c)
δF
μ= . (2.17d)
δφ
This model assumes that the flowis incompressible everywhere at the expense of local mass conservation
insidetheinterfacial/mixing region.
Forthecompressible models,wedefinethetotalenergy as
(cid:2)
ρ(cid:4)v(cid:4)2
E(t)= [ +F]dx, (2.18)
2
Ω
wherexistheEuleriancoordinate. Forcompressiblemodel1,therateofchangeinthetotalenergyisgiven
by
(cid:2) (cid:2)
dE ρ −ρ ∂p
=− [λ(cid:4)∇μ(cid:4)2+(τ−τ ):D]dx+ 1 2 φ(v·∇p+ )dx. (2.19)
c
dt ρ ∂t
Ω 2 Ω
Toensurepositivity inthefirstintegral(2ηD+νtr(D)I):D≥0,η≥0,ν+2η ≥0.Unfortunately, wehave
3
nocontrol overthesignforthesecond integral. Forcompressible model2,therateofchange intheenergy
forthissystemofgoverningequations isgivenby
(cid:2) (cid:2)
dE ρ −ρ ∂p dφ
=− [λ(cid:4)∇μ(cid:4)2+(τ−τ ):D]dx+ 1 2 φ( − )dx. (2.20)
c
dt ρ ∂t dt
Ω 2 Ω
Ineithermodels,theenergydissipationcannotberigorouslyestablished. Ontheotherhand,fortheincom-
pressible model,theenergyisdefinedby
(cid:2)
ρ(cid:4)v(cid:4)2
E(t)= [ +F]dx. (2.21)
2
Ω
With the help of the divergence-free condition in the continuity equation, Shen and Yang [32, 33] showed
anenergylawisavailable fortheincompressible system:
(cid:2)
dE
=− [λ(cid:4)∇μ(cid:4)2+(τ−τ ):D]dx≤0. (2.22)
c
dt
Ω
One of the traded-offs in enforcing the mass conservation is losing the energy law. Therefore, both com-
pressiblemodelsalongwiththeincompressiblemixturemodelin[33]arevalidapproximately. Wenotethat
ananalogous phasefieldequation (tomodel1)wasalsoderivedusinganenergyargument in[28].
5
2.2 Choice offree energy
Thefreeenergy F cantakedifferent formdepending ontheapplications. Inthispaper, weconsider the
freeenergy densityinthefollowingform:
1
F(φ,∇φ)=k Tγ[ (cid:4)∇φ(cid:4)2+ f(φ)], (2.23)
B
2
where k is the Boltzmann constant, T is the absolute temperature, and γ is a parameter with the unit of a
B
numberdensityperunitlength. γisinfactproportionaltotheproductofthenumberdensityperunitvolume
andthesquareofthepersistent length.
WefirstlookattheGinzburg-Laudau freeenergy with
1
f(φ)= φ2(1−φ)2 (2.24)
ε2
fortwoimmisciblefluids,whereε>0isasmallparametercharacterizingthehydrophobicpropertybetween
thetwofluids. Therefore,
δF 1
=k Tγ[−∇2φ+ φ(1−φ)(1−2φ)]. (2.25)
δφ B ε2
Wealsoconsider theFlory-Huggins mixingfreeenergy fortwoimmisciblefluidstosimulate thephase
separation dynamics. TheFlory-Huggins mixingfreeenergydensityisgivenby(2.23)with
1 φ 1−φ
f(φ)= [ lnφ+ ln(1−φ)+χφ(1−φ)], (2.26)
ε2 N N
1 2
whereN andN arethepolymerizationindicesforfluid1andfluid2andχisthemixingparameterbetween
1 2
0and2. Ifbothareviscous fluids,weassumeN =N =1. Inthiscase,
1 2
δF 1 lnφ ln(1−φ)
=k Tγ[−∇2φ+ ( + +χ(1−2φ))]+const. (2.27)
δφ B ε2 N N
1 2
2.3 Non-dimensionalization
We denote the characteristic time scale by t and length scale by L . The dimensionless variables are
0 0
definedby
t x vt pt2
t˜= ,x˜ = ,v˜ = 0,p˜= 0 . (2.28)
t0 L0 L0 ρ2L20
Wewill drop the˜on the dimensionless variables inthe following. Wechoose L so that the dimensionless
0
length L =1. We use L = 1 in the following calculations simply for convenience. The dimensionless
y x
modelparameters aredefinedby
Rei,s = ρt02ηLi0, Rei,v = ρt02νLi0, (i=1,2)Λ= λt0LkB40Tγ, ε˜ = Lε0,
(2.29)
Bi=φρ1 +(1−φ), 1 = φ + 1−φ, 1 = φ + 1−φ.
ρ2 Res Re1,s Re2,s Rev Re1,v Re2,v
Here Re and Re denotes the Reynolds number corresponding to the shear and volumetric stress, and Λ is
s v
thedimensionless mobilityparameter. WesetL =1,C= kBTγt0 =1inthisstudyyieldingt = L40ρ2.
y L40ρ2 0 kBTγ
6
Thedimensionless equations forthetwocompressible modelsaregivenby
⎧
⎪⎪⎪ φt+∇·(φv)=∇·(Λ∇μ),
⎪
⎪
⎪
⎪
⎪
⎪⎨ Bi[v +v·∇v]=∇·( 1 ∇v)+∇(( 1 + 1 )∇·v)−∇p−φ∇μ,
t Res Res Rev
(2.30)
⎪
⎪⎪⎪ ∇·v=(1−ρ1)[∇·(Λ∇μ)],
⎪⎪⎪ ρ2
⎪
⎪
⎩
μ=−Δφ+ f˜(φ)(Model1), μ=−Δφ+ f˜(φ)+(ρ1 −1)p(Model2),
ρ2
where f˜(φ)aregivenby(2.24)or(2.26)withεreplacedbyε˜.
The above system is subjected to a set of suitable initial and boundary conditions. For example, if the
mixtureisconfinedinadomainΩ,theboundary conditions are
∂φ ∂μ
| = | =0, v| =0, (2.31)
∂n ∂Ω ∂n ∂Ω ∂Ω
wherenistheoutwardnormal.
3 Numerical schemes
Shen and Yang have developed and studied efficient numerical schemes for the incompressible phase
fieldmodelingreatdetailin[32,33]. Inthissection,wefocusonextendingtheirresultstothecompressible
models given for example by (2.15) for the binary fluid mixture. Since the two compressible models do
not admit a dissipative energy law, it is not possible to construct energy stable schemes for them. Instead,
we shall construct efficient and easy to implement numerical schemes to accommodates the non-vanishing
divergence velocityfieldsoastopreservethemassconservation.
3.1 Discretization intime
Tosimplify thepresentation, weshallpresent onlyfirst-order schemes. Inwhatfollows,thesuperscript
ndenotesthetimelevelandΔt isthetimestepsize.
Schemebasedonamodifiedprojection:
1. Solve(φn+1,μn+1)from:
φn+1−φn +∇·(φnvn)=∇·(Λ∇μn+1), ∂φn+1| =0,
Δt ∂n ∂Ω
(3.1)
μn+1=−Δφn+1+ f˜(φn)+ S (φn+1−φn), ∂μ| =0,
ε˜2 ∂n ∂Ω
whereSisacomputational parameterandεistheparameterinthefreeenergy. Thelasttermisadded
tostabilizetheschemetoallowlargerstepsizes. Itsroleistodampthehighfrequencyorshortwaves
inthenumericalsimulation.
2. Denote
Bin+1=φn+1ρ1 +(1−φn+1), Rens+1=φn+1Re1,s+(1−φn+1)Re2,s,
ρ2 (3.2)
Renν+1=φn+1Re1,ν+(1−φn+1)Re2,ν;
7