Diffusion dominated evaporation in multicomponent lattice Boltzmann simulations

0.48 MB·

by  Dennis Hessling

Diffusion dominated evaporation in multicomponent lattice Boltzmann simulations
Dennis Hessling,1,2,∗ Qingguang Xie,2,† and Jens Harting3,2,‡
1Materials innovations institute (M2i), Elektronicaweg 25, 2628 XG Delft,Netherlands
2Department of Applied Physics, Eindhoven University of Technology,
P.O. Box 513, NL-5600MB Eindhoven, The Netherlands
3Helmholtz Institute Erlangen-Nu¨rnberg for Renewable Energy (IEK-11),
Forschungszentrum Ju¨lich, Fu¨rther Straße 248, 90429 Nu¨rnberg, Germany
We present a diffusion dominated evaporation model using the popular pseudopotential multicom-
ponentlatticeBoltzmannmethodintroducedbyShanandChen. Withananalyticalcomputationof
thediffusioncoefficients,wedemonstratethatFick’slawisobeyed. Wethenvalidatetheapplicability
of our model by demonstrating the agreement of the time evolution of the interface position of an
7
evaporating planar film to the analytical prediction. Furthermore, we study the evaporation of a
1
freely floating droplet and confirm that the effect of Laplace pressure is significant for predicting the
0
time evolution of small droplet radii.
2

n PACSnumbers: 47.11.-j,77.84.Nh.
a
J

1 I. INTRODUCTION of individual fluid constituents independently. They can
3 thus help to improve our understanding of evaporation
Evaporatingfluidsareubiquitousinourdailylifeandin driven fluid transport. Simulations of evaporating fluids
n] industrialprocesses,suchasinkjetprinting[1],coating[2] often utilize molecular dynamics (MD) [10–13]. While
y andparticledeposition[3]. Inparticularforsuspensionsor MD offers a very high flexibility in the microscopic de-
d polymersolutions,aswellasfluidsinconfinedgeometries, tails, its computational cost is very high. Therefore, MD
- the evaporation of individual components can induce simulations are limited to very small length and time
u
l fluid flows or a change of relative concentrations leading scales on the nanometre or nanosecond scale [13]. In or-
f to changing rheological and transport properties of the der to reach experimentally relevant scales, a continuum
.
s constituents. For example, the evaporation of a sessile approach is more productive and our method of choice is
c
i colloidal droplet on a substrate leads to a capillary flow the lattice Boltzmann method (LBM) [14–16]. The LBM
s transportingthecolloidalparticlestotheedgeofadroplet, hasgainedpopularityforthesimulationoffluidflowsdue
y
h which finally results in a ring-like deposit [4]. The ring- to its straightforward implementation and parallelization.
p like stains can be a useful tool to deposit particles and Soon after its invention, the LBM was extended to sim-
[ can also be disadvantageous when a uniform pattern is ulate multiple interacting fluid phases and components
desirable. Another example is the evaporation of droplets and today a plethora of multiphase and multicomponent
2
onroughorchemicallypatternedsubstrates. Surrounding methods exists [17–21, 27].
v
7 geometries and the wettability of a substrate have a large The mesoscale nature of the method combined with
3 influence on the lifetime of evaporating droplets [5]. A the possibility to add additional fields, external forces,
6 thorough understanding of this impact of evaporation on suspended objects, thermal noise, or complex boundary
3 the fluid behavior is mandatory to consequently optimize conditionsinaverystraightforwardmannerhasmadethe
0
industrial applications and to improve our fundamental LBMparticularlypopularforapplicationsinmicrofluidics
.
1 understandingofeffectslikefilmformation,dropletdrying, and soft matter physics. Many of the physical systems
0 or droplet spreading. studiedinthesefieldsincludevolatileliquids,wheretheef-
7
There are numerous theoretical [6–8] and experimen- fect of evaporation plays a dominant role. Therefore, it is
1
tal [4, 9] studies of fluid evaporation. While most theo- notsurprisingthatanumberofgroupshassimulatedevap-
:
v retical studies are limited to the macroscopic scale, ex- oratingfluidsusingtheLBMrecently. Ledesma-Aguilaret
Xi periments suffer from difficulties that arise by tuning the al. [22, 23] present a diffusion based evaporation method
individual microscale properties of fluids. The thorough based on the free energy multiphase lattice Boltzmann
r
a understanding of fluid evaporation calls for mesoscopic or method and demonstrate quantitative agreement with
microscopic details and the flexibility to tune the proper- several benchmark cases as well as qualitative agreement
ties of individual fluid constituents independently. This with the experimental data of evaporating droplet arrays.
is possible by means of computer simulations. Computer Jansen et al. [24] study the evaporation of droplets on
simulations allow access to parameters which are not eas- a chemically patterned substrate and qualitatively com-
ily controllable in experiments and to tune the properties parethesimulationresultswithexperimentaldata. Their
method is based on a continuous removal of mass from
the droplet and thus does not allow studying transport
[email protected] processes in the vapor phase. Yan et al. [25] present a
∗
[email protected] thermal model to study the contact line dynamics during
†
‡ [email protected] droplet evaporation where the liquid-vapor phase change
2
is driven by a temperature field and a well defined equa- equilibrium distribution function
tion of state. Joshi and Sun [26] present simulations of
(cid:20) c ·uc (cid:0)uc ·uc (cid:1)
drying colloidal suspensions by means of a modified pseu- feq(ρc,uc )=ω ρc 1+ i eq − eq eq
dopotential multiphase model following Shan and Chen. i eq i c2 2c2
s s
They assign a fixed mass flux to the system boundary (cid:0)c ·uc (cid:1)2(cid:21)
which causes a reduction in vapor concentration and thus + i eq , (2)
2c4
triggersaliquid-vaporphasechangeattheinterface. How- s
ever, their results are purely qualitative and a thorough
analytical understanding of the diffusion in the system is wherecs = √13∆∆xt isthespeedofsoundandωi isaweight
factor defined as ω = 1, ω = 1 and ω = 1 .
missing. The densities are 0defin3ed a1,s...,ρ6c(x,1t8) = ρ 7(cid:80),...,1f8c(x,3t6),
Inthispaperweovercomethislimitationandintroduce 0 i i
where ρ is a reference density, and the velocities are
an alternative evaporation model for the pseudopoten- defined a0s uc(x,t)=(cid:80) fc(x,t)c /ρc(x,t), while the ve-
tial method of Shan and Chen. We focus on the two- i i i
locity in the equilibrium distribution function is uc =
componentversionofthemethod[27],buttheapplication (cid:80) ρcuc/(cid:80) ρc. eq
toanarbitrarynumberofcomponentsandthemultiphase c c
For brevity and numerical efficiency we choose the
pseudopotentialmethodisstraightforward. Generally,the
lattice constant ∆x, the timestep ∆t, the unit mass ρ
pseudopotentialLBMisverypopularduetoitseaseofim- 0
and the relaxation time τc to be unity, which leads to a
plementation and flexibility when combined for example
kinematic viscosity νc = 2τ−1 = 1 in lattice units.
with complex geometries [17, 28], locally varying contact 6 6
The system boundaries are treated as periodic bound-
angles [29], or suspended particles [30]. To trigger evap-
ariesbydefault. Todoso,fluidleavingonesystembound-
oration, we do not impose a mass flux, but instead fix
ary reenters the opposite side and forces are computed
the density of one component at selected boundary sites
acrosstheseperiodicboundaries. Toinhibitflowwallscan
whichinducesadensitygradient. Theevaporationprocess
beconstructedbyinvertingvelocitiesatselectedboundary
is diffusion dominated and can be well described using
sites [14].
Fick’s law with well defined diffusivities. We validate the
applicability of our model by comparing the time depen-
dent simulation results of an evaporating planar film and
B. The pseudopotential multicomponent lattice
a freely floating evaporating droplet with their respective
Boltzmann method
analytical predictions.
This remainder of this paper is organised as follows. For the fluid components introduced above to become
Sec. II introduces the lattice Boltzmann method and our immiscible, Shan and Chen introduced a pseudopotential
extension for evaporating fluids. Our results are shown interaction force
in Sec. III and Sec. IV concludes the article.
(cid:88)(cid:88)
Fc(x,t)=−Ψc(x,t) ω gcc¯Ψc¯(x+e ,t)e (3)
i i i
c¯ i
to achieve separation of the components [27]. This force
II. SIMULATION METHOD is defined as a nearest neighbor interaction between fluid
components c and c¯[27] and scaled through the choice
A. The lattice Boltzmann method of the parameter gcc¯. Here Ψc(x,t) is an effective mass,
defined as
The lattice Boltzmann equation can be obtained from Ψc(x,t)≡Ψ(ρc(x,t))=1−e−ρc(x,t)/ρ0. (4)
spatiallyandtemporallydiscretizingtheBoltzmannequa-
The force is applied to the fluid by adding a shift of
tion. Multiple fluid components c are modeled by fol-
lowing the evolution of the single particle distribution ∆uc(x,t) = τcρFc(cx(,xt,)t) to uceq(x,t) during collision. This
function causes the separation of fluids and the formation of a
diffuseinterfacebetweenthem. Thewidthoftheinterface
separating the regions is typically about 5∆x [31], with a
fc(x+e ∆t,t+∆t)−fc(x,t)= −∆t[fc(x,t)
i i i τc i small dependence on the interaction strength.
−feq(ρc(x,t),uc (x,t))]. (1)
i eq
The single particle distribution functions fc(x,t) at posi- C. The evaporation model
i
tions x alternatively stream, as described by the LHS of
(1), along the i=1,...,19 discretized directions e and When the interaction parameter gcc¯ in the pseudopo-
i
collide,asdescribedbytheRHSof (1)ateverytimestept. tential model is properly chosen [28], a separation of
Throughout this work we utilize 2 components c and c. components takes place. Each component will separate
The collision is achieved by relaxing the probability dis- into a denser majority phase of density ρ and a lighter
ma
tribution functions towards a discretized second-order minority phase of density ρ , respectively.
mi
3
In order to drive the system out of equilibrium, we 10−3
·
impose the density of component c at the boundary sites
x to be of constant value ρc(x ,t)=ρc by setting the 1 Simulation Timesteps
diHstribution function of componHent c toH Theory 100
200
400
fc(x ,t)=feq(ρc ,uc (x ,t)), (5)
i H i H H H 800
cρ 0.5 m e 1600
in which ucH(xH,t) = 0. Setting the velocity to zero at Ti 3200
the system boundary is in agreement with the idea of an
undisturbed large volume outside the system, which does
not cause perturbations in the system. Depending on
the ratio of the minority density ρc and ρc this induces 0
mi H
evaporation or condensation. Furthermore, for simplicity,
0 20 40 60 80 100 120
we ensure total mass conservation within the system by
x
setting the density of component c¯as
FIG. 1: Time evolution of density profiles of component
ρc¯(xH,t)=ρc(xH,t−1)+ρc¯(xH,t−1)−ρcH (6) c. The system is initially empty of c. Diffusion allows
fluid from an infinite reservoir at x=0 to build up
and again ensure an undisturbed flow field by setting gradually changing density profiles in agreement with the
ucH¯ (xH,t)=0, so that the distribution functions of com- analytical solution.
ponent c¯at the evaporation boundary sites x become
H
fc¯(x ,t)=feq(cid:0)ρc¯ ,uc¯ (x ,t)(cid:1). (7)
i H i H H H velocitiesaredefinedasanaverageofthetotalmomentum
before and after each collision [33] as
We note that our model can be easily extended to situa-
tions where mass is not conserved. We ensure the equiv- ρcuc+ρc¯uc¯+ 1(Fc+Fc¯)
alence to the open system by mimicking a zero density U = 2 , (10)
ρc+ρc¯
gradient within the infinite volume outside the system.
1 (cid:20) ρc(ρcuc+ρc¯uc¯)(cid:21)
This way forces from the pseudopotential interactions Uc = ρcuc+Fc+ . (11)
only have an impact in the bulk of the simulation volume 2ρc ρc+ρc¯
and not at the boundary. This can be achieved either
ByperformingaChapman-Enskogexpansion[33],Eq.(9)
by using evaporation boundaries on periodic sites as well
can be rewritten to be identical to Eq. (8), with the
or by using a second layer of evaporation boundary sites.
diffusion coefficients given as [32]
Thereby we enforce the same density and a zero gradient
at the system boundary. (cid:18) 1(cid:19) ρc¯ c2ρcΨc¯gc¯cΨ(cid:48)c
In the case where the set density ρc is not equal to the Dcc = c2 τ − − s ,
H s 2 ρc+ρc¯ ρc+ρc¯
equilibrium minority density ρc , a density gradient in
the vapor phase of component cmiis formed. This gradient (cid:18) 1(cid:19) ρc c2ρc¯Ψcgcc¯Ψ(cid:48)c¯
Dcc¯ = −c2 τ − + s , (12)
causes component c to diffuse towards the minimum. s 2 ρc+ρc¯ ρc+ρc¯
with Ψ(cid:48)c being the spatial derivative of Ψc. We note that
the diffusion is dependent on the symmetric interaction
III. RESULTS AND DISCUSSIONS
strengths gcc¯ and gc¯c, as well as the densities of the two
components ρc and ρc¯.
A. Diffusion
In the limit of small gradients, we can assume that
∇ρc¯(x,t)=−∇ρc(x,t), with which Eq. (8) becomes
In binary fluid mixtures, following Fick’s first law, we
can write the mass flux of component c as jc =−Dc∇ρc, (13)
where
jc =−Dcc∇ρc−Dcc¯∇ρc¯, (8)
(cid:20) 1 c2 (cid:21)
where Dcc, Dcc¯ are the diffusion coefficients. In the Dc = c2s(τ−2)−ρc+sρc¯(ρc¯Ψcgcc¯Ψ(cid:48)c¯+ρcΨc¯gc¯cΨ(cid:48)c) .(14)
Shan-Chen multicomponent method, the mass flux of
component c can be written as [32] To validate the theoretical analysis above, we investi-
gate the diffusion of a component c into a system filled
jc =ρc(Uc−U), (9) withanothercomponentc¯. Weperformasimulationwith
a system size of 125×4×4 and fill the system with fluid
where Uc and U are macroscopic velocities of component c¯of density ρc¯=0.7. A wall of thickness h=2 is placed
c and the binary mixture, respectively. The macroscopic at x=125 and the fluid interaction parameter is set to
4
t = 0
t=0
ρc =0.03
120 H
ρc =0.02
2d H c
0 ρ = 0.03
x
H
121000 H
x
I c
ρ = 0.02
80
H
FIG. 2: (Color online) Schematic representation of the x
60
planar film (front view). We fill the lower half of the
system with fluid c and the upper half with fluid c such
that a fluid-fluid interface forms at xI. The interface 10400
thickness is 2d . To drive the evaporation we impose the
0
boundary condition ρc(x=x )=ρc at the top of the
H H 20
system. A solid surface with thickness h=2 is located at
the bottom.
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
gcc¯=gc¯c =3.6. We utilize an evaporation boundary at
x=0andsetthedensityρc =0.001toensureadiffusive 80 ρc
H
flow that is undisturbed by convection. Then component
FIG. 3: Density profile of fluid c along the x direction
c diffuses into the system. Meanwhile, we numerically
after equilibration (defined as t=0, solid line) and
solve Fick’s second law
density profiles at t=106 timesteps later with boundary
∂ρc densities ρc =0.03 (dashed line) and ρc =0.02 (dotted
=−Dcc∆ρc−Dcc¯∆ρc¯ H H
∂t x line). The magnification depicts of the subtle difference
∂ρc¯ of ρc causing a different density gradient and a different
=−Dc¯c¯∆ρc¯−Dc¯c∆ρc (15) 6H0
∂t time behavior of the moving interface.
todescribethespaceandtimedependentdensityprofileof
componentc. InFig.1wecomparethelatticeBoltzmann
conditions at the bottom, parallel to the interface, while
simulation results (symbols) with the numerical solution
the boundaries normal to the substrate are periodic.
of Eq. (15) (solid lines). From the evaporation boundary
After equilibration the density of fluid c is constant in
with density ρc fluid diffuses into the system. Being a
H boththedenserphase(ρc ≈0.704)andthelighterphase
diffusion process, the rate at which the fluid invades the ma
(ρc4≈00.036), whereas between them a diffuse interface
systemisdependentonthedensitygradient. Thegradient mi
of about 2d =5 lattice units is formed, as shown in the
subsequently decreases and fluid distributes itself further 0
density profile along x direction in Fig. 3 (solid line). We
into the system, aiming to remove the gradient. There
thenapplytheevaporationboundaryconditionbysetting
is a good agreement between simulation results and the
the density at the top boundary ρc(x = 128) to ρc . In
numerical solution, as shown in Fig. 1. H
Fig. 3 we show the density profiles along the x direction
We note that the diffusion equation does not hold at
just after equilibration (solid line) and for evaporation
fluid-fluid interfaces [32]. However, the movement of the
boundarydensitiesρc =0.03(dashedline)andρc =0.02
interface during evaporation is governed by diffusion of (do2tte0d line) after 1H06 subsequent simulation tiHmesteps.
thefluidssurroundingit,whichwedemonstrateasfollows.
Adensitygradientoffluidcisformedinthelighterphase,
resulting in diffusion of fluid c towards the evaporation
boundary. Thus, the interface position decreases with
B. Evaporation of a planar film
time. Itdecreasesfasterforalowerevaporationboundary
density ρc , which indicates that the mass flux increases
We investigate the evaporation of a planar film sitting H
with decreasing the evaporation boundary density.
on a solid substrate, as illustrated in Fig. 2. To do so
If we assume that the fluid densities in the minority
0
we perform simulations with a system size of 128×4×4.
phases vary linearly, Eq. (13) becomes
We fill one half of the system with fluid c and the other
hρcalf=wiρtc¯h fl=ui0d.0c¯4)osfuecqhutahladteansfliutyid(-flρcmuiad=intρercm¯faace=fo0r.m70s, jc =−(Dc(ρmi−ρcH)/(xH −xI −d0))n, (16)
mi mi
at x = 64. We define the position of the interface x where n is the normal vector of the interface. The mass
0 I 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
as the position of ρc−ρc¯=0. The interaction strength flux is approximately proportional to the density differ-
in Eq. (3) is chosen to be gcc¯ = gc¯c = 3.6. We place a ence between the minority density and the evaporation
wall of thickness 2 with simple bounce back boundary boundary density ρ −ρc , which allows to control the
mi H c
ρ
5
0
x
/ 1
I
x
t
h
g
ei
h
R
e I
c R
a H
erf 0.5
t
n
edi ρcH =0.035 2d0
maliz ρρcHc ==00..003200
H
r
o
N
0 0.2 0.4 0.6 0.8 1
FIG. 5: (Color online) Schematic cross-sectional
Timesteps 107 representation of a droplet of radius R surrounded by
· I
another fluid up to the spherical system boundary R ,
FIG. 4: (Color online) Interface position as a function of H
where the evaporation boundary condition is imposed.
time for different evaporation boundary densities
ρc =0.035, ρc =0.03 and ρc =0.02. The theoretical
H H H
prediction Eq. (22) (solid lines) agrees well with the
ρc = 0.035, ρc = 0.03 and ρc = 0.02 along with our
simulation data (symbols) H H H
theoretical analysis Eq. (22) are presented in Fig. 4. We
find excellent quantitative agreement between theory and
simulation.
evaporation rate by varying ρc .
H
Withtheassumptionthatthedensityprofileacrossthe
interface is also linear, the total mass of fluid c in the
C. Evaporation of a freely suspended droplet
system is
Mc(t)=A(cid:2)(x −d )ρc +(x −x −d ) Inthissectionweinvestigatetheevaporationofafreely
I 0 ma H I 0
(ρc +ρc )/2+d (ρc −ρc )(cid:3), (17) floating droplet. A droplet of component c with a radius
mi H 0 ma mi of R is the center of a spherical system of size R and
I H
where A is the area of the cross-section. From Eq. (17), surrounded by component c¯, as shown in Fig. 5. A spher-
we can obtain ical evaporation boundary is applied at RH. Under the
assumption of quasi-static dynamics, the density profile
dx
dMc/dt=A(ρc −ρc /2−ρc /2) I. (18) of component c in the lighter phase satisfies the Laplace
ma mi H dt equation,
Based on the principle of mass conservation, the time
∆ρc(r)=0, (23)
evolution of mass obeys
where the boundary conditions are
(cid:90) t
Mc(t)=Mc(0)−A jc·ndt, (19) ρc(r)| =ρc (24)
0 r=RI+d0 mi
where Mc(0) is the initial mass of fluid c. From Eq. (19), and
we can also get the time derivative of the total mass as ρc(r)| =ρc . (25)
r=RH H
dMc/dt=A|jc|. (20) In spherical coordinates, we obtain the analytical solu-
tion as
By comparing Eq. (18) and Eq. (20), we obtain
R −r R +d
dxI = D(ρcmi−ρcH) . (21) ρc(r)=ρcH−(ρcH−ρcmi)RH −HRI −d0 I r 0, (26)
dt (x −x −d )(ρc −ρc /2−ρc /2)
H I 0 ma mi H where r is the distance from the center of a spherical
WesolveEq.(21)withtheinitialconditionx (t=0)=x coordinate system, originating at the droplet center. In-
I 0
and finally obtain the interface position as a function of serting Eq. (26) into Eq. (13), we obtain the mass flux as
time
x (t)=x +d −(cid:2)(x +d −x )2 jc(r)=−Dc(ρ −ρc ) RH(RI +d0) n , (27)
I H 0 H 0 0 mi H (R −R −d )r2 r
+2 D(ρcmi−ρcH) t(cid:3)1/2. (22) where n is the normal vector Hto theIdrop0let interface.
ρc −ρc /2−ρc /2 r
ma mi H Assuming the density profile across the interface also
Thesimulationresultsofthetimeevolutionoftheinter- satisfies Laplace’s equation, we obtain the total mass of
face position for different evaporation boundary densities fluid c in the system as
6
(cid:90) RI−d0 (cid:90) RI+d0 (cid:18) R +d −rR −d (cid:19)
Mc = 4πr2ρc dr+ 4πr2 ρc −(ρc −ρc ) I 0 I 0 dr
ma mi mi ma 2d r
0 RI−d0 0
(cid:90) RH (cid:18) R −r R +d (cid:19)
+ 4πr2 ρ −(ρ −ρc ) H I 0 dr. (28)
H H mi R −R −d r
RI+d0 H I 0
We simplify Eq. (28) and derive the time derivative of
ρc =0.030
the total mass as 1 ρHc =0.025
H
ρc =0.020
0 H
2π(cid:18) R
dMc/dt= 3 RH2 (−ρH +ρcmi) R/I 0.8
s
+R (−2R ρc +2R ρc −2d ρc +2d ρc ) u
H I H I mi 0 H 0 mi di
(cid:19) a
r
+6R2ρc −6R2ρc −2d2ρc +2d2ρc dR /dt. t 0.6
I ma I mi 0 ma 0 mi I e
pl
(29) ro
d
d
Based on the principle of mass conservation we have ze 0.4
ali
(cid:90) t m
Mc(t)=Mc(0)− 4πr2jc·n dt, (30) or
r N
0 0.2
where Mc(0) is the total initial mass of fluid c in the
system. From Eq. (30) with using Eq. (27), we also get
the time derivative of the total mass as
0 0.2 0.4 0.6 0.8 1
dMc/dt=4πDc(ρmi−ρcH)(R(RI−+Rd0)−RHd ). (31) Timestep ·106
H I 0
FIG. 6: (Color online) Time evolution of the droplet
BycomparingEq.(29)andEq.(31)weobtainthetime radius for evaporation boundary densities ρc =0.03,
H
evolution of the droplet radius as ρc =0.025 and ρc =0.02. Our theoretical analysis
H H
without surface tension Eq. (32) (dashed lines) agrees
dR /dt=4πDc(ρ −ρc )(R +d )R quantitatively with the simulation data (symbols) for
I mi H I 0 H
(cid:20) (cid:18) large droplet radii and deviates for small droplet radii.
2π
(R −R −d ) R2 (−ρ +ρc ) Our theoretical analysis which includes the surface
H I 0 3 H H mi
tension Eq. (40) (solid lines) agrees quantitatively well
+R (−2R ρc +2R ρc −2d ρc +2d ρc ) with the simulation data for both large and small droplet
H I H I mi 0 H 0 mi
(cid:19)(cid:21)−1 radii.
+6R2ρc −6R2ρc −2d2ρc +2d2ρc . (32)
I ma I mi 0 ma 0 mi
We solve Eq. (32) numerically with a 4th-order Runge- evolution of the droplet radius well, and quantitatively
Kutta algorithm. For the simulations we initialize a agreeswiththesimulationdataforthedropletatalarger
droplet with a radius of R = 65 and densities ρc = radius. However, it deviates for small droplet radii. This
0 ma
ρc¯ =0.70,ρc =ρc¯ =0.04,inacomputationaldomain can be explained by the fact that we neglected the effect
ma mi mi
of2563. Thedensitiesoffluidcequilibratetoρc ≈0.712 ofsurfacetensiononthedropletevaporation. Thesurface
ma
inside the droplet and ρc ≈0.037 outside. We initiate tension induces a Laplace pressure, which is larger when
mi
the evaporation by setting the density at the spherical the droplet radius becomes small [6]. We can take into
evaporation boundary r =R to ρc =0.03, ρc =0.025, account this effect as follows:
H H H
and ρc = 0.02. In Fig. 6 we compare the analytical ForasphericaldroplettheYoung-Laplaceequationcan
H
solutionEq.(32)(dashedlines)withthesimulationresults be written as
(symbols). We note that Shan-Chen models are exposed
to spurious vaporisation effects once the droplets become 2γ
P(r >R ,t)=P(r <R ,t)− , (33)
small, i.e. when the diameter is around 5−10 lattice I I R
I
units. To avoid the effect of the spurious vaporisation on
the analysis, we only use the simulation data when the where γ is the surface tension, P(r > R ,t) and P(r <
I
diameter of the droplet is larger than 20 The analytical R ,t) are the pressures outside and inside the droplet at
I
solution captures the qualitative features of the time time t, respectively. We can write the pressure inside the
7
droplet as [27] ByinsertingEq.(35)intoEq.(36),weobtainthemajority
density of fluid c inside the droplet as
c2
P(r <R ,t)=c2(ρc +ρc¯ )+ sg Ψ(ρc )Ψ(ρc¯ ).(34) 2γ (cid:18) 1 1 (cid:19)
I s ma mi 2 cc¯ ma mi ρc (t)=ρc (t=0)− − . (37)
ma ma c2 R R
s 0 I
For simplification, in the case of ρc¯ (cid:28)ρc we can write
mi ma The minority density of fluid c outside the droplet can
the pressure in terms of the leading term as
betreatedasproportionaltothemajoritydensityoffluid
c inside the droplet [6]. Thus we obtain
P(r <R ,t)=c2ρc . (35)
I s ma
ρc (t=0)− 2γ( 1 − 1 )
The pressure outside the droplet can be treated as con- ρc (t)=ρc (t=0) ma c2s R0 RI . (38)
stant during evaporation, so that we get mi mi ρcma(t=0)
2γ Forbrevity,wedenoteρcma(t=0)asρcma,0 andρcmi(t=0)
P(r >RI,t) = P(r <RI,t=0)− R as ρcmi,0. We insert Eq. (37) and Eq. (38) into Eq. (28),
0
and after some manipulations, we finally obtain the time
2γ
= P(r <R ,t)− . (36) derivative of the droplet mass as
I R
I
(cid:18)
2π
dMc/dt = (−R2 −2R R −2R d )ρc
3 H I H H 0 H
(cid:32) ρc − 2γ( 1 − 1 )(cid:33)
+ (R2 +2R R +2R d −6R2+2d2) ρc ma,0 c2s R0 RI
H I H H 0 I mi,0 ρc
ma,0
(cid:32) (cid:33)
−2γρc 1
+ (R R2 +dR2 +R2R +2R R d +d2R −2R3+2R d2) mi,0
I H H I H I H 0 0 H I I c2ρc R2
s ma,0 I
(cid:19)
2γ 1 1 −2γ
+ (6R2−2d2)(ρc − ( − ))+(2R3−2R d2)( ) dR /dt. (39)
I ma,0 c2 R R I I c2R2 I
s 0 I s I
WecompareEq.(39)withEq.(31)andgettheequation
for dR /dt including the effect of surface tension as
I
(cid:20) (cid:18)
2π
dR /dt = 4πDc(ρ −ρc )(R +d )R (R −R −d ) (−R2 −2R R −2R d )ρc
I mi H I 0 H H I 0 3 H I H H 0 H
(cid:32) ρc − 2γ( 1 − 1 )(cid:33)
+ (R2 +2R R +2R d −6R2+2d2) ρc ma,0 c2s R0 RI
H I H H 0 I mi,0 ρc
ma,0
(cid:32) (cid:33)
−2γρc 1
+ (R R2 +dR2 +R2R +2R R d +d2R −2R3+2R d2) mi,0
I H H I H I H 0 0 H I I c2ρc R2
s ma,0 I
2γ 1 1 −2γ (cid:19)(cid:21)−1
+ (6R2−2d2)(ρc − ( − ))+(2R3−2R d2)( ) . (40)
I ma,0 c2 R R I I c2R2
s 0 I s I
WesolveEq.(40)numericallyandcomparethetheoretical be neglected. This result is of particular importance for
prediction with simulation data in Fig. 6. The theoreti- lattice Boltzmann simulations of evaporating droplets
cal analysis including the effect of surface tension (solid since the typical number of lattice nodes available to re-
lines) agrees quantitatively well with the simulation data solve the radius of a single droplet is often limited. This
(symbols) for both large and small droplet radii. Thus, holds in particular for systems involving a large number
we confirm that the effect of surface tension becomes sig- of droplets.
nificant when the droplets become small and must not
8
IV. CONCLUSION other fluid, as an extension of the famous Epstein-Plesset
theory [6]. While the original publication assumes an
infinite system, we extended the model towards a finite
We presented a diffusion dominated evaporation model system size. We demonstrate a good agreement between
usingthepopularpseudopotentialmulticomponentlattice theory and simulation if one takes into account the effect
Boltzmann method introduced by Shan and Chen. The of surface tension causing a high Laplace pressure and
evaporation is induced by imposing the density of one an increased evaporation rate in the case of small droplet
component at the system boundary while ensuring total radii.
mass conservation, which causes diffusion of components As an outlook we note that our method is not only
driven by a density gradient. The diffusion coefficients suitable to simulate evaporating fluids, but that it is
depend on the densities of the fluids as well as the in- straightforwardtoapplyittoinvestigatethecondensation
teraction strength parameters of the Shan-Chen model. of droplets. Therefore, our method can be a powerful
With the analytically determined diffusion coefficients we tool for exploring both evaporation and condensation
confirm that the diffusion obeys Fick’s law. processes in complex fluidic systems.
We derived a theoretical model for the time evolution
of the interface position of an evaporating planar film
ACKNOWLEDGMENTS
under the quasi-static assumption. Our theoretical model
predicts that the evaporation flux is proportional to the
D.HesslingandQ.Xiecontributedequallytothiswork.
density difference between the minority density and the
Financial support from the Materials innovation institute
evaporation boundary density ρ −ρc , while the time
mi H (M2i, www.m2i.nl, project number M61.2.12454b), Oc´e-
evolution of the interface position obeys the expected t0.5
TechnologiesB.V.andNWO/STW(STWprojectnumber
law. We then carried out simulations which are in good
13291) as well as the allocation of computing time at the
quantitative agreement with our analytical model.
High Performance Computing Center Stuttgart and the
Furthermore, we derived analytical models describing Ju¨lich Supercomputing Centre are highly acknowledged.
the evaporation of a floating droplet surrounded by an- We thank D. Lohse and I. Jimidar for fruitful discussions.
[1] L. Wu, Z. Dong, M. Kuang, Y. Li, F. Li, L. Jiang, and A. Narv´aez, B. D. Jones, J. R. Williams, A. J. Valoc-
Y. Song. Adv. Funct. Mater., 25:2237, 2015. chi, and J. Harting. Comput. Geosci., 20:777, 2016.
[2] C. N. Kaplan, N. Wu, S. Mandre, J. Aizenberg, and [18] A. K. Gunstensen, D. H. Rothman, S. Zaleski, and
L. Mahadevan. Phys. Fluids, 27:092105, 2015. G. Zanetti. Phys. Rev. A, 43:4320, 1991.
[3] W. Han and Z. Lin. Angew. Chem. Int. Ed., 51:1534, [19] M. R. Swift, W. R. Osborn, and J. M. Yeomans. Phys.
2011. Rev. Lett., 75:830, 1995.
[4] R.D.Deegan,O.Bakajin,T.F.Dupont,G.Huber,S.R. [20] X. He, Sh. Chen, and R. Zhang. J. Comput. Phys.,
Nagel, and T. A. Witten. Nature, 389:827, 1997. 663:642, 1999.
[5] H. Gelderblom, A. G. Mar´ın, H. Nair, A. van Houselt, [21] G. Falcucci, S. Ubertini, C. Biscarini, S. Di Francesco,
L. Lefferts, J. H. Snoeijer, and D. Lohse. Phys. Rev. E, D. Chiappini, S. Palpacelli, A. De Maio, and S. Succi.
83:026306, 2011. Commun. Comput. Phys., 9:269, 2011.
[6] P.S.EpsteinandM.S.Plesset. J. Chem. Phys.,18:1505, [22] R. Ledesma-Aguilar, D. Vella, and J. M. Yeomans. Soft
1950. Matter, 10:8267, 2014.
[7] Y. O. Popov. Phys. Rev. E, 71:036313, 2005. [23] G. Laghezza, E. Dietrich, J. M. Yeomans, R. Ledesma-
[8] J. M. Stauber, S. K. Wilson, B. R. Duffy, and K. Sefiane. Aguilar, E. S. Kooij, H. J. W. Zandvliet, and D. Lohse.
Langmuir, 31:3653, 2015. Soft Matter, 12:5787, 2016.
[9] A.G.Mar´ın,H.Gelderblom,D.Lohse,andJ.H.Snoeijer. [24] H. P. Jansen, K. Sotthewes, J. van Swigchem, H. J. W.
Phys. Rev. Lett., 107:085502, 2011. Zandvliet,andE.S.Kooij. Phys.Rev.E,88:013008,2013.
[10] L.N.Long,M.M.Micci,andB.C.Wong. Comput.Phys. [25] Q. Li, P. Zhou, and H. J. Yan. Langmuir, 32:9389, 2016.
Commun., 96:167, 1996. [26] A. S. Joshi and Y. Sun. Phys. Rev. E, 82:041401, 2010.
[11] F. Wang and H. Wu. Soft Matter, 9:5703, 2013. [27] X. Shan and H. Chen. Phys. Rev. E, 47:1815, 1993.
[12] W. Chen, J. Koplik, and I. Kretzschmar. Phys. Rev. E, [28] N. Martys and H. Chen. Phys. Rev. E, 53:743, 1996.
87:052404, 2013. [29] J. Harting, C. Kunert, and H. J. Herrmann. Europhys.
[13] D. Lohse and X. Zhang. Rev. Mod. Phys., 87:981, 2015. Lett., 75:328, 2006.
[14] S. Succi. The Lattice Boltzmann Equation: For Fluid [30] F. Jansen and J. Harting. Phys. Rev. E, 83:046707, 2011.
Dynamics and Beyond. Oxford University Press, 2001. [31] S. Frijters, F. Gu¨nther, and J. Harting. Soft Matter,
[15] R.Benzi,S.Succi,andM.Vergassola.Phys.Rep.,222:145, 8:6542, 2012.
1992. [32] X. Shan and G. Doolen. Phys. Rev. E, 54:3614, 1996.
[16] D.Raabe. Modell. Simul. Mater. Sci. Eng.,12:R13,2004. [33] X. Shan and G. Doolen. J. Stat. Phys., 81:379, 1995.
[17] H. Liu, Q. Kang, C. R. Leonardi, S. Schmieschek,