Modeling Evaporation and Particle Assembly in ... - ACS Publications

Article

Modeling Evaporation and Particle Assembly in Colloidal Droplets Mingfei Zhao, and Xin Yong Langmuir, Just Accepted Manuscript • Publication Date (Web): 26 May 2017 Downloaded from http://pubs.acs.org on May 29, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Langmuir is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Modeling Evaporation and Particle Assembly in Colloidal Droplets Mingfei Zhao and Xin Yong* Department of Mechanical Engineering, Binghamton University, Binghamton, New York 13902, United States Abstract: Evaporation-induced assembly of nanoparticles in a drying droplet is of great importance in many engineering applications, including printing, coating and thin film processing. The investigation of particle dynamics in evaporating droplets can provide fundamental hydrodynamic insight for revealing the processing-structure relationship in the particle self-organization induced by solvent evaporation. We develop a free-energy-based multiphase lattice Boltzmann method coupled with Brownian dynamics to simulate evaporating colloidal droplets on solid substrates with specified wetting properties. The influence of interface-bound nanoparticles on the surface tension and evaporation of a flat liquid-vapor interface is first quantified. The results indicate that the particles at the interface reduce surface tension and enhance evaporation flux. For evaporating particle-covered droplets on substrates with different wetting properties, we characterize the increase of evaporate rate via measuring droplet volume. We find that droplet evaporation is determined by the number density and circumferential distribution of interfacial particles. We further correlate particle dynamics and assembly to the evaporation-induced convection in the bulk and on the surface of droplet. Finally, we observe distinct final deposits from evaporating colloidal droplets with bulkdispersed and interface-bound particles. In addition, the deposit pattern is also influenced by the equilibrium contact angle of droplet. * E-mail [email protected]

1 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 35

I. Introduction Nanotechnology has created unprecedented functional materials with self-assembled, dissipative one-dimensional or two-dimensional structures characterized by nanometric length scales. These nanomaterials find important applications in biological, optical, and electrical devices,1–3 yet a substantial fraction of them awaits associated manufacturing processes capable of scaling up to mass production. Soft solution processing has been widely recognized as a simple and facile method for synthesizing nanomaterials, utilizing the evaporation of a liquid solution containing nonvolatile solutes.4,5 For example, a droplet with stable nanoparticle (NP) dispersions leaves a colloidal assembly on a solid surface after drying. The self-assembled structures are kinetically trapped in non-equilibrium states, which are sensitively determined by the physicochemical dynamics of the evaporative process.6,7 In particular, the non-uniform deposition of particles within the initial wetted region delineated by the contact line, termed as the coffee-ring effect, has been extensively studied in the past decades.8,9 Inhomogeneous morphologies (e.g., a single “coffee-ring”, stochastically distributed concentric rings, cellular network, etc.)10–12 are often undesirable for surface patterning, coating, and printing applications. The ability to suppress the coffee-ring effect and achieve well-controlled deposit structures is crucial for emerging processing techniques. The recent progress in understanding the mechanism of the coffee-ring formation has elucidated the complex interplay between internal flow fields and particle-substrate interactions.8,12–18 Namely, Deegan et al. attributed the peripheral ring deposition to the outward capillary flow toward the pinned contact line due to faster evaporation at the droplet periphery.8 Hu and Larson discovered that the formation of coffee-ring deposits not only requires pinned contact line but also suppression of Marangoni flow driven by surface tension gradient along the free surface of the droplet. Such Marangoni stresses result from the temperature gradient induced by evaporative cooling.13 Further harnessing of Marangoni effects through binary solvents and additional surface-active surfactants can prevent particle deposition at the droplet perimeter and yield relatively uniform deposits.14,17,18 The interactions between suspended particles and the substrate may also play an important role in the deposition.14,15


Page 3 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Despite the rich inventory of literature on the transport phenomena in drying droplets, the relation between the spatial distribution of particles and the deposition pattern is incomplete. While most studies consider particles dispersed in the bulk of the droplet, few recognized the importance of particle adsorption at the liquid-vapor interface and discussed the deposition of interface-bound particles.19–26 Lin and Jaeger19–21 synthesized highly uniform, long-rangeordered monolayers in colloidal droplet evaporation by leveraging the interface-capturing of surface-active NPs and the ensuing two-dimensional templated assembly at the liquid-vapor interface. Yunker et al.22 suppressed the ring deposition through the jamming of anisotropic particles on the droplet surface. Recent theoretical study by Kang et al.23 further highlights the inevitable particle capturing by the receding liquid-vapor interface and the dominant role of the flow at the free surface in determining the deposition patterns. To our best knowledge, the deposition of colloidal droplets with particles only at the interface has not been explicitly investigated. Here, we examine the evaporation-driven self-assembly and deposition of particlecovered sessile droplets with no particle dispersed in the bulk. To model evaporating colloidal droplets, we develop a coupled lattice Boltzmann-Brownian dynamics model. We investigate the evaporation kinetics in our model and reveal the influence of NPs adsorbed at the liquid-vapor interface. A comprehensive study of the evaporating sessile droplets with NPs on the free surface and the resulting deposit structure is conducted. Below, we first describe the details of the numerical method and then discuss our findings on colloidal droplet evaporation and deposition with only interface-bound NPs. II. Computational Method Traditional computational fluid dynamics methods for multiphase flows generally require specific algorithms for tracking the fluid-fluid interface, such as the level set method for tracking the interface position27 and the volume of fluid method for tracking the volume fraction of each phase or component.28,29 The lattice Boltzmann method (LBM) does not directly solve the governing equations of macroscopic properties derived from conservation laws. Instead, it models a fluid as fictitious particles whose dynamics is described by the discrete Boltzmann equation. Among various LBMs that model multiphase fluids, the free-energy-based formulations30–34 are widely applied due to its thermodynamic consistency and the ability to 3 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 35

explicitly define surface tension and wetting properties of solid boundary.35,36 In this approach, the separation of fluid phases is spontaneous and a diffuse interface evolves naturally according to the transport of fluid particles without the implementation of specific interfacial tracking and reconstruction schemes. Thus, the LBM is computationally efficient in modeling multiphase systems and has been proven to be a powerful tool to investigate the interactions between suspended particles and the surrounding fluid.37–41 To simulate a two-phase fluid filled with NPs, we combine the Zheng-Shu-Chew freeenergy LBM33 (referred to as the ZSC model below) with Brownian dynamics (BD). The hydrodynamics of the NP-laden two-phase fluid is described by the continuity equation, the Navier-Stokes equation, and the Cahn-Hilliard equation, which is the evolution equation for the interface. The three governing equations are given respectively as follows: ∂n r + ∇ ⋅ ( nu ) = 0 ∂t

(1)

r ∂ ( nu ) rr r ur r + ∇ ⋅ ( nuu ) = −∇ ⋅ P + η∇ 2u + G + G drag ∂t

(2)

∂φ r + ∇ ⋅ (φ u ) = θ M ∇ 2 µφ ∂t

(3)

In the ZSC model, n is the fluid density set as the mean value of the two different fluids,33

n = ( ρ H + ρ L ) 2 , where ρ H and ρ L are the densities of high-density phase and low-density r phase, respectively. u is the local velocity, η is the dynamic viscosity, and P is the pressure r r tensor. G and G drag are the forces on the fluid due to the presence of NPs. φ is the order parameter, which takes two distinct values in each of the phases φ * and -φ * , where

φ * = ( ρ H − ρ L ) 2 . θ M is the mobility of the order parameter and µφ is the chemical potential. The corresponding lattice Boltzmann equations33 (LBEs) given in the ZSC model are expressed as follows: r r r 1 r r f i ( x + eiδ t , t + δ t ) − f i ( x , t ) =  f i eq ( x , t ) − f i ( x , t )  + Si τ

(4)

n


Page 5 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

1 r r r r r r r r gi ( x + eiδ t , t + δ t ) =gi ( x , t ) + (1 − q )  gi ( x + eiδ t , t ) − gi ( x , t )  + ( gieq ( x , t ) − gi ( x , t ) ) (5)

τφ

Here, f i and gi are two particle distribution functions, defining macroscopic hydrodynamic variables. In particular, density is defined as n = ∑ i f i and momentum is obtained by nuα = ∑ i f i eiα . The order parameter is calculated by the second particle distribution function as

φ = ∑ gi . Equation 4 is the LBE of the flow field, from which Eqs. 1 and 2 can be recovered by i

the Chapman-Enskog expansion. Equation 5 of the phase field recovers the Cahn-Hilliard equation and thereby captures the motion of the interface. τ n and τ φ are the respective relaxation times for the two LBE, and

is a constant coefficient determined by the value of τ φ ,

q = 1 (τ φ + 0.5 ) . The mobility parameter is given by θ M = q (τ φ q − 0.5 ) Γ , where coefficient Γ is

introduced to control the mobility. f i eq and gieq are the corresponding equilibrium particle distributions. f i eq is defined as f i eq =ωi Ai + ωi n ( 3eiα uα − 3u 2 / 2 + 9uα uβ eiα eiβ / 2 ) , where the coefficients are taken as A1 = 9n / 4 − 15 (φµφ + n / 3) 4 and Ai

i = 2,...,9

= 3 (φµφ + n / 3) . gieq is

r r calculated as g ieq =Ai + Biφ +Ciφ ei ⋅ u , where A1 = −2Γµφ , Ai i ≠1 = Γµφ / 2 , B1 = 1, Bi i ≠1 = 0 , and

Ci = 1/ 2q . The forcing term Si in Eq. 4 is given by42 Si =

eiα ωi ( Fα + µφ ∂αφ ) cs2

(6)

where α stands for the α -direction component, ωi is the weight factor, and cs is the speed of sound. In the two-dimensional nine velocity components (D2Q9) discretization of twodimensional systems,43 we have ω1 =4/9 , ωi

i = 2,...,5

= 1/ 9 , and ωi

i = 6,...,9

= 1/ 36 . In Eq. 6, Fα is the

α -component body force and µφ ∂α φ stands for the α -component of surface force.

The

thermodynamic properties of the system enter the governing equations through the pressure tensor P and the chemical potential , which is derived from a given free energy functional.44 The free energy functional of the NP-laden fluid is given by


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 35

κ n ln n  2  F = ∫ dV ψ (φ ) + ( ∇φ ) + + ψ p dV + ∫ψ s (φs ) dS 2 3  ∫ 

(7)

In the first integral,ψ (φ ) =A(φ 2 − φ *2 ) 2 is the bulk free energy density of the two-phase fluid,33

µφ where A is an amplitude parameter. The second term κ / 2 ( ∇φ ) represents the interfacial 2

free energy density. The thickness of interface W and the surface tension σ are described in terms of A and κ as W = 2κ A φ * and σ = 4 2κ Aφ *3 3 , respectively. The expressions of the pressure

tensor

is

given

by

Pαβ = pδαβ − κ ( ∇φ ) δαβ − ∇α φ∇ β φ    2

p = A ( 3φ 4 − 2φ *2φ 2 -φ *4 ) − κφ∇ 2φ + ( ∇φ ) / 2 + n / 3 2

.

The

chemical

,

potential

where is

µφ = ∂F ∂t = A ( 4φ 3 − 4φ *2φ ) − κ∇ 2φ + ∂ψ p ∂φ . The first-order and second-order derivatives of the order parameter are calculated using the isotropic finite difference scheme developed by Lee and Lin.32 The second integral in Eq. 7 captures the interaction between the NPs and the background fluid, in which we apply a particle-fluid free energy density45,46 r r

r

r

ψ p = ∑ V0 exp ( − r − ri r0 ) φ ( r , t ) − φi ( ri , t ) 

2

(8)

i

r In this equation, the summation is taken over all particles located at corresponding positions ri .

r

r

φ ( r , t ) is the order parameter of the fluid at the position r obtained by linear interpolation. r

φi ( ri , t ) is the order parameter of the i th particle, which defines its wetting behavior. In order to minimize the total free energy of the system, the NPs will move toward the region where the surrounding fluid has the same order parameter as they do. Namely, the NPs with φi = φ * or -φ * is preferentially wetted by the corresponding phase, while φi = 0 results in the NPs that like the interface. Experimentally, the wetting properties of NPs can be varied by surface chemistry modification. V0 is a coefficient that dictates the strength of the NP-fluid interaction. The locality of the NP-fluid interaction is imparted by the exponential function in Eq. (8), which describes the decay of particle influence on the fluid dynamics as the distance between a NP and a fluid element increases. r0 in the exponential function controls the influence range and strength. 6 ACS Paragon Plus Environment

Page 7 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Herein, we consider small NPs with weak coupling between the particle and the fluid. In other words, the presence of NPs should not significantly disturb the surrounding fluid, thus small values of V0 and r0 are chosen.45,46 The characterization of the NP-fluid interaction parameters is detailed in Supporting Information SI1. Furthermore, the interaction between NPs and fluid is treated differently for the liquid r and vapor phases. The order parameter of the fluid at the position of the NP φ ( ri , t ) is used to define the background fluid phase for each NP. When φ ∈ ( −0.99φ * , 0.99φ * ) , the NP is within the interfacial region. The NP is considered to be located in the vapor phase if φ ≤ −0.99φ * and in the liquid phase if φ ≥ 0.99φ * . When particles are in the bulk liquid phase or within the liquidvapor interfacial region, the two-way coupling between NPs and surrounding fluid is imposed through the particle-fluid forces. Having NPs, the external body force Fα in Eq. 6 corresponds to r r r the two forcing terms G and G drag in the Navier-Stokes equation (Eq. 2). The forcing term G accounts for the NP-fluid interaction induced by the particle wetting, and we use the expression45,46 r G = −φ∇ ( ∂ψ p ∂φ )

(9)

r r r By substituting the ψ p expression in Eq. 8 to Eq. 9, we can obtain G = ∑ f i ( r , t ) , where i

r r r r r r f i ( r , t ) = 2φ∇ exp ( − r − ri r0 ) φ ( r , t ) − φi ( ri , t ) 

(10)

r The second forcing term G drag is viscous drag, which can be expressed by45,47 r r r r G drag = −∑ δ ( r − ri )F drag ,i

(11)

i

r Both equations have the summation taken over all the NPs. In Eq. 11, δ ( r ) is the Dirac delta function, reflecting the locality of the reaction acting on the surrounding fluid from the drag r r force on the i th particle F drag ,i . δ ( r ) is approximated in the calculation by a linear extrapolation to the nearest points. The drag force is given as r r r r F drag ,i = −ζ δ r δ t − u ( ri , t ) 

(12)


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 35

r where ζ =6πη RP is the Stoke’s drag coefficient with R p being NP radius, δ r δ t is the velocity r r r of the i th particle, and u ( ri , t ) is the fluid velocity at the position of the i th particle r . The third integral in Eq. 7 accounts for the wetting properties of solid surface on which the sessile droplet sits. ψ s (φs ) is the free energy density of the surface with the form

ψ s (φs ) =-hφs .35 Minimizing the free energy yields an equilibrium boundary condition on the r r surface S , κ n ⋅ ( ∇φ ) s = ∂ψ s ∂φs , where n is the unit normal vector to the surface. Substituting the surface free energy density into this expression leads to the wetting boundary condition,

( ∂φ

∂ n ) s = − h κ . The relationship between the parameter h and the equilibrium contact angle

θeq is given by35 h = 2 (φ * )

2

π   α  α  2κ A sgn  − θ eq  cos   1 − cos    2   3   3 

(13)

where α = arccos ( sin 2 θ eq ) and sgn ( x ) is the sign function. Thus, the first-order derivative in the perpendicular direction of solid surface on the boundary has a fix value of − h κ . For the second-order derivative on the boundary, we use the standard right-handed finite-difference scheme provided by Briant et al.48 By applying the wetting boundary condition, the contact angle is controlled in the simulation while the contact line is free to move. According to the Newton’s second law, the dynamics of NPs in the bulk liquid phase and within the liquid-vapor interfacial region is described by the overdamped Langevin dynamics in the Lagrangian frame of reference as follows,45

r r dt r r r r dri = u ( ri , t ) dt + 2 D p dW ( t ) +  Fi p ( t ) + Fi e ( t ) 

ζ

(14)

Here, we consider the inertia of NPs is negligible compared to the drag force. The positions of NPs are thus updated through Eq. 14 using a second-order Runge-Kutta method.45 The first term r r is the advection of NPs in the surrounding fluid. u ( ri , t ) is the velocity of the surrounding fluid

r at the position of the i th particle ri , which is obtained from the lattice Boltzmann fluid solver. The second term describes the Brownian motion of the NPs, where D p is the diffusion 8 ACS Paragon Plus Environment

Page 9 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

r coefficient of the NPs. W ( t ) is the Wiener process, which is a Gaussian random variable with r r 2 variance W ( t + dt ) − W ( t ) = N ∆t ( N is the dimension of the system) and satisfies the r rr r fluctuation-dissipation theorem.45,46 Fi P ( t ) = ∫ dr f i ( r , t ) is the total force on the i th particle

r r r from the surrounding fluid, where f i ( r , t ) is given by Eq. (11). Fi e captures the excluded volume interactions between the i th particle and other neighboring particles and is defined as

Fi e ( t ) = −

N

∑

(

r r ∂ψ m ri − rj

j =1, j ≠ i

) ∂ ( rr − rr ) i

(15)

j

The summation runs over all the neighboring j particles. To model the particle-particle interaction, we select the repulsive part of the Morse potential

ψ m ( r ) = ε (1 − exp  −λ ( r − re )  ) for r < re 2

(16)

where ε and λ are constants determining the strength and range of the interaction between different particles. re is the specific distance where the force between particles is equal to zero, beyond which the Morse potential transforms from repulsion to attraction. The NPs experience attraction when they move to the vicinity of the wall. The attraction is also expressed by

(

Fi ,ay ( t ) = ∂ψ a ri , y − rwall , y

) ∂( r

i. y

− rwall , y

(

where the similar Morse potential47 ψ a ( r ) = ε 1 − exp  −λ ( r − ra ) 

)

2

)

(17)

for r < ra is used.

We want to point out that the ZSC model uses the average mass density in the LBEs, which is subjected to small variations in the whole flow field. Consequently, the effect of local density and viscosity variations cannot be properly captured in the momentum equation.49,50 We expect the density and viscosity contrast mainly influences Stokes’ drag force and Brownian motion of a particle moving across the liquid-vapor interface. Thus, to capture this behavior, the dynamics and particle-fluid coupling for particles present in the vapor phase are treated differently. In particular, the particle-fluid coupling is removed. In the fluid solver, the two body r r forces G and G drag induced by particles do not exist in the vapor phase. Equation 14 still describes the particle dynamics in the vapor phase, however, without the fluid-particle 9 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 35

r interaction Fi p and the Brownian term. Instead, a long-range attraction representing van der Waals (VDW) attraction between NPs and a phenomenological body force in the negative y direction are introduced. The attractive part of Morse potential with different parameters is used to represent the VDW interaction. The attraction has an effective range re < r < rc , where

rc = re + ( ln 2 ) / λ .47,51 The vertical body force is added due to both physical and computation considerations. Physically, it is common to ignore the gravitational force on nanoparticles suspended in liquid for nanoscopic systems. However, the particles migrating into the vapor phase will certainly experience much stronger gravity, which could significantly influence the dynamics. The vertical force is applied to capture this phenomenon. In the simulation, a particle moving into the vapor phase may only experience van der Waals attraction from neighboring particles, or it will not move. Without the vertical force, the particle will never be absorbed again to the receding interface during evaporation. The physically rooted body force also serves to prevent this simulation artifact of floating particles. In addition, all the NPs will be fixed if their distances between the solid substrate are smaller than d = 0.5 regardless of which phase the NPs belong to, representing permanent deposition to the solid surface. Other depositing particles can lay on the top of the first layer of fixed particles due to the vertical force and particle-particle interactions. We note that the multiphase LBMs incorporating with large density and/or viscosity ratios

34,50,52–54

could capture more precise particle dynamics in the diffusive interfacial region

and the vapor phase. In this work, we expect the density and viscosity contrast would not qualitatively influence the evaporation-driven hydrodynamics in the liquid phase and the concurrent particle self-assembly, under the quasi-steady low-rate diffusion-dominant evaporation in our simulation. However, the density and viscosity contrast will probably play a more significant role in the convection-dominant and unsteady evaporation, which is not considered in this work. In the free-energy LBM, we can simulate the droplet evaporation under isothermal conditions with low evaporation rate (see Supporting Information SI1 for the physical value of evaporate rate) by imposing a Dirichlet boundary condition on the order parameter φ

SH

= φH ,

where S H is an evaporation boundary enclosing the system and the order parameter at the 10 ACS Paragon Plus Environment

Page 11 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

boundary φH < −φ * .55 Here, the implementation of evaporation can be interpreted as imposing a boundary of constant vapor concentration while solving the diffusion equation for the vapor concentration. For an evaporating flat interface studied below, the flat evaporation boundary

φH = −1.5φ * for the outer domain y ≥ 80 is assigned. For droplets, we assign φH = −1.5φ * out a semicircular evaporation boundary with a radius of 98 LB units enclosing the droplet. Our results show that this evaporation model captures the differential evaporation flux along the droplet surface. The simulated velocity fields induced by droplet evaporation using this scheme are also in good agreement with the analytical flow field considering moving contact line,56 which will be further elaborated in the results. This evaporation model also captures the differential evaporation flux along the droplet surface. Our typical simulation domain is 200 ×100 in LB units. The top and bottom boundaries apply the full-way bounce-back boundary condition to recover the no-slip conditions. The horizontal boundaries are periodic boundaries. The binary fluids with different dimensionless density ρ H = 100 and ρ L = 1 are modeled, which results in the equilibrium order parameter

φ * = 49.5 for the liquid phase and -49.5 for the vapor phase and the average density n = 50.5 . We set the thickness of interface W = 5 LB units and the dimensionless surface tension σ =0.1 . The relaxation times are chosen as τ n = 0.875 and τ φ =0.7 . The mobility coefficient is Γ=40 . The viscosity can be calculated by η =n (τ n − 0.5) cs2 , where the speed of sound cs = 3 . The dimensionless radius of NP is R p = 0.032 . In the following simulation, we set the simulation parameters as ε repulsive =1.0 , ε attractive =1.0 × 10 −4 , λ =0.6 , re = 1.2 , and ra = 0.6 . k BT can be converted into LB units as 2.7003 × 10−4 LB unit . According to the Einstein relation D p = k BT

( 6πη R

particle

),

we have the particle diffusivity in LB units as D p = 7.919 ×10−5 .

Considering the much weaker particle absorption energy applied in the simulation than that in physical systems (see Supporting Information SI1 for details), we set the magnitude of the vertical body force to be 3.6 × 10−4 for mimicking the deposition in the vapor phase. The fluid dynamics is solved by the LBM and the particle dynamics is solved by BD. One LB time step contains 20 particle time steps dt . Unless otherwise stated, all the initial droplets have a radius of Rdroplet = 50 and requires 2 × 105 LB time steps to reach the equilibrium state with the 11 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 35

specified contact angle θeq . For static drops, the production runs are carried out for 3.4 × 106 time steps. When evaporation is added to the system, the total time steps for complete dry out of droplets is approximately 8 × 106 . The relation between the dimensionless simulation parameters to physical values is detailed in Supporting Information SI1. An in-house C code of the LBMBD model was developed to carry out the present study and generate the following results. The Laplace’s law of a 2D stationary droplets was verified in our LB model as a typical benchmark for multiphase flows. The 2D Brownian diffusion of the particles in a quiescent LB was verified by the relationship between mean square displacement of the particles and the set particle r diffusivity d 2 = 4 D p t .

III. Results and Discussion A. Evaporating Flat Interfaces with Nanoparticles As discussed in the Introduction section, the primary interest of this work is the dynamics of neutral-wetting NPs ( φ p = 0 ) in evaporating droplets, which preferentially stay on the droplet surface. Based on the comprehensive analysis of the effects of particle-fluid interaction parameters and particle number for static droplets (see Supporting Information SI2), we select the simulation parameters V0 = 5 × 10 −7 and r0 = 0.6 for the evaporating systems. We first examine the effect of interface-bound NPs ( φ p = 0 ) on the evaporation of a flat liquid-vapor interface, where the particle number density (per length) remains constant during evaporation. To model evaporation, we assign Dirichlet boundary condition with φH = −1.5φ * for all the nodes with y ≥ 80 . The evaporation flux from the Cahn-Hilliard equation is defined as J α = −θ M ∂α µ .55 Because the flat interface is along the x direction, we only consider the y component of the evaporation flux normal to the interface. Figure 1 demonstrates the variations of order parameter, evaporation flux, and chemical potential profiles in time for a particle-free flat interface. The order parameter profiles clearly show the reduction of interface height as evaporation continues. As shown in Fig. 1(b), the flux profile shows a peak value near the liquid-vapor interface, followed by a gradual decrease to a smaller value defined as the evaporation boundary. Both the peak and boundary values of the flux reduce as evaporation continue. In addition, the gradient of flux in the vapor region also reduces in time. The behavior of flux is directly related to the 12 ACS Paragon Plus Environment

Page 13 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

chemical potential profiles in Fig.1(c). The chemical potential profile in the vapor phase far from the liquid-vapor interface is monotonic and approximately linear due to the diffusion of phase field induced by the boundary-driven gradient in order parameter shown in Fig. 1(a). The chemical potential exhibits significant variations near the interface, where there exists a maximum gradient in the profile corresponding to the peak flux in Fig. 1(b). Below, our results show that the peak flux is more sensitively influenced by the NPs compared to the boundary flux, so we quantify the evaporation flux using the peak value for comparing interfaces with different

N.


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 35

Figure 1. Profiles of (a) order parameter, (b) y-component evaporation flux and (c) chemical potential at 1× 106 (blue), 2 × 106 (green), and 3.4 × 106 (red) time steps, along the y direction for an evaporating flat interface without NPs. The error bars show the standard deviations along the x direction.


Page 15 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Figure 2. Dependence of (a) averaged y-position of the flat interface and (b) peak evaporation flux at the last time step of the simulations on particle number density ρ NP . The error bars show the standard deviations along the x direction.

Figure 2 plots the average peak value of the evaporation flux at the final time step of the simulations ( 3.4 × 106 time steps) as a function of average particle number density ρ NP (defined as number of particles per units length), where the error bars represent the variations along the x direction. Compared with the particle-free ( ρ NP = 0 ) case, evaporation is noticeably enhanced by adding NPs at the interface, as evidenced by the reduced y position of the interface and the increased peak evaporation flux. The peak flux increases with increasing ρ NP , but the enhancement is not appreciable compared to the substantial difference between the particle-free interface and the interface with the minimum number of particles N = 50 ( ρ NP = 0.25 ). In


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 35

addition, as ρ NP increases, the larger error bars imply that the existence of more NPs at the interface leads to larger perturbation on the phase field and the evaporation flux.

Figure 3. Time evolutions of (a, b) peak evaporation flux and (c, d) interface position. (b) Magnified view of the region in (a) enclosed by the dashed line. (d) Magnified view of the region in (c) enclosed by the dashed line. The error bars represent the standard deviations along the x direction.

To further elucidate the modulation of evaporation by neutral-wetting NPs, a detailed comparison among simulations with different particle number densities is conducted. Time evolutions of the peak flux and the interface positions for the particle-free, N = 50 , and N = 300 ( ρ NP = 1.5) interfaces are presented in Fig. 3. We observe that all peak flux profiles 16 ACS Paragon Plus Environment

Page 17 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

start from the same point, corresponding to the beginning of evaporation. Evaporation is initiated at the evaporation boundary with the Dirichlet condition and gradually propagates toward the liquid-vapor interface, as shown in Fig. S7. The sudden drop of the peak flux in Fig. 3(a) at the early stage is associated to the transient development of evaporation at the liquid-vapor interface. During the propagation, all the simulations produce similar average peak fluxes. Once the evaporation flux reaches the liquid-vapor interface, the peak flux becomes dependent on the condition of the interface. In particular, the evolutions of the peak fluxes for the interfaces with different particle numbers are distinct from each other. Figure 3(b) shows the magnified view of the flux variation at the late stage of the simulation. The evaporation enhancement between the particle-free case and the particle-laden case with ρ NP = 0.25 is apparent, but the difference between ρ NP = 0.25 and ρ NP = 1.5 interfaces is less pronounced. This contrast is consistent with the final flux profile shown in Fig. 2. Corresponding to the propagation process show in Fig. 3(a), the interface positions in Fig. 3(c) keep constant at the early stage. Once evaporation occurs at the liquid-vapor interface, the interface heights start to decrease and the particle-laden interfaces decreases significantly faster than the particle-free interface. Figure 3(d) is the enlarged view of the dotted box in Fig. 3(c). For the late stage, the position of particle-free interface is evidently higher than those of particle-laden interfaces. There is no statistically significant difference between the average interface positions for ρ NP = 1.5 and ρ NP = 0.25 , which may be attributed to the undulations of the particle-laden interfaces. From the analysis of the flat interface, we conclude that the presence of NPs enhances evaporation but the influence of particle number density is weak for homogeneously distributed particles. Note that, the particles in our system are treated as point masses without physical volumes, and thereby do not induce the surface coverage effect, which is the major contribution to evaporation retardation of interface-bound NPs as pointed out by previous MD and MDPD simulations.57,58 However, the coupling between the NPs and the fluid indeed reduces surface tension (as detailed in Supporting Information SI2) and leads to lower energy barrier for phase transition. Overall, evaporation is promoted due to the NP-induced local variation of the fluid phase field. The neutral-wetting NPs in this model can be considered as an analogy to the surfactant monolayers that increase evaporation rate.59

B. Evaporating Particle-Covered Droplets 17 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 35

Figure 4． Time evolutions of the volume changes of particle-laden droplets compared to the particle-free droplet with different equilibrium contact angles of (a) θ eq = 60° , (b) θ eq = 90° , and (c) θ eq = 120° . The dashed lines indicate the early stage of evaporation.

In our simulations, we control the wetting properties of the solid substrate through a given equilibrium contact angle, and therefore the contact line is allowed to move freely. All the droplets with different equilibrium contact angles have the same initial configuration, namely a Rdroplet = 50 droplet with 90° contact angle. In this manner, the volumes of the sessile droplets are


Page 19 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

kept as constant, and the time required for complete evaporation is similar. We assign the evaporation Dirichlet boundary condition with φH = −1.5φ * for the fluid lattices outside a semicircle centered at (100,0 ) with radius Revap = 98 . For all particle-covered droplets, the NPs are randomly introduced in the interfacial region after the droplet reaches the equilibrium state. In the droplet evaporation simulation, the peak evaporation flux cannot be easily defined due to the curved liquid-vapor interface and the influence of solid substrate. Instead, we monitor the volume of droplet to characterize the evaporation. We count the total lattice points that are occupied by the liquid phase to approximate the droplet volume. To quantify the enhancement of evaporate rate due to NPs, a volume difference ∆Ω between the particle-covered droplet and particle-free droplet is then defined as Ω p − Ω 0 for each time step, where Ω p and Ω 0 are the volumes of particle-covered and particle-free droplets, respectively. Figures 4 shows the variation of ∆Ω in time for droplets with different equilibrium contact angles θ eq = 60° , 90° , and 120° . The straight blue lines with ∆Ω = 0 represent the result of the particle-free droplet and serve as the reference. The instantaneous evaporation rate enhancement ( d ∆Ω / dt ) can be characterized by the slope of the volume difference curve. Overall, the results are consistent with our observation in the simulations with evaporating flat interfaces − the interface-bound NPs enhance the evaporation of liquid phase. The more NPs on the surface, the faster droplets evaporate. At the early stage of the droplet evaporation, the profiles of all the three particle-covered sessile droplets exhibit the same behavior as those of the particle-laden flat interface. Namely, at the beginning of droplet evaporation, the volume differences between the particle-free and particle-covered droplets are quite significant, while the difference among particle-covered droplets with increasing particle numbers is not obvious. However, as droplet evaporation continues, the evaporation enhancement starts to show pronounced dependence on the equilibrium contact angle θeq and the particle number N . As shown in Fig. 4(a), the increase of evaporation rate due to NPs for

θ eq = 60° droplets generally diminishes along with time, which is manifested by the reducing slope of ∆Ω curve. However, the droplet with N = 200 NPs exhibits almost constant slope throughout the entire evaporation. Figure 4(b) indicates that the evaporation of θ eq = 90° particlecovered droplets becomes significantly faster after the early stage (slope increases monotonically 19 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 35

after 2 × 106 time steps). Generally, the droplet with the most particles has the steepest slope of the volume difference for the late stage. θ eq = 120° droplets have similar trends among different particle numbers as in θ eq = 60° droplet. As shown in Fig. 4(c), the volume difference profile of the N = 200 droplet has a very large slope at the late stage, while the slopes of the other droplets with fewer particles remain approximately constant throughout evaporation. The evaporation behavior is connected with the particle distribution as elaborated below. Detailed observation of the particle-covered droplet reveals that the interfacial NPs exhibit different distributions on the surface for evaporating droplets having different θeq . Once a quasi-steady evaporation state is established, we demonstrate one representative snapshot for each equilibrium contact angle as shown in Fig. 5. The interfacial particles on the surface of the

θ eq = 60° droplet accumulate around the contact line areas. In contrast, the θ eq = 120° droplet has an apparently high particle concentration at the apex of the droplet. The particle distribution at the surface of the θ eq = 90° droplet is homogenous. Below, we link the advection induced by droplet evaporation to the dynamics of particle assembly to understand the observed distribution patterns. The flow fields inside the droplet and at the liquid-vapor interface are plotted in Fig. 6 for the particle-free and particle-covered droplets. The spurious current in our model is quantified by modeling a stationary droplet on the hydrophilic substrate θ eq = 60° and shows a magnitude on the order of 10−6 . After evaporation starts, the velocity field demonstrated in Figure 6 is on the order of 10−3 . This difference confirms that the spurious current will not influence the evaporation-driven flow. The evaporation of both droplets with θ eq = 60° generates a surface flow in the interfacial region pointing to the contact line region and a clockwise circulation in the bulk, as shown in Figs. 6(a, d). The surface flow will entrain the NPs to the contact line, leading to higher particle concentration at the contact line as shown in Fig. 5(a). For droplets with 90° contact angle in Figs. 6(b, d), the flow fields in the interfacial region are normal to the interface, while the flows in the bulk are much weaker than in θ eq = 60° droplet. The lack of tangential advection results in the uniform distribution of NPs along the interface in Fig. 5(c). The droplet with equilibrium contact angle θ eq = 120° exhibit a different flow pattern compared to θ eq = 60° 20 ACS Paragon Plus Environment

Page 21 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

droplet and θ eq = 90° droplet. Figures 6(c, f) show that the velocity in the interfacial region points to the apex of the droplet and the bulk flow circulation is counterclockwise. Although the existence of NPs induces perturbation on the flow field, the evaporation-drive advection in the particle-covered droplets is consistent with that in the particle-free droplets.

Figure 5. Snapshots of N = 200 evaporating sessile droplets with different equilibrium contact angles of (a) θ eq = 60° , (c) θ eq = 90° , and (d) θ eq = 120° . (b) Magnified view of the contact line region in (a) enclosed by the dashed line. The color maps are constructed based on the value of the fluid order parameter. The green line represents the position of the evaporation boundary. The circular markers represent NPs. The colors of the markers’ edges stand for different particle states: blue for the moving particles and green for the fixed particles. The colors of the marker’s


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 35

faces stand for different regions the particles belong to: red for the liquid-vapor interface and white for the vapor phase. The same markers are used for the following figures with NPs.

The non-uniform distribution of particles on the droplet surface influences the droplet evaporation. When the contact line recedes, the NPs in the interfacial region may deposit to the substrate driven by the short-range substrate attraction, given that the NPs are present in the vicinity of the substrate. In other words, only particles in the contact line region may be deposited. For θ eq = 60° droplets, particles aggregate in the contact line region, which leads to a significant number of NPs deposited to the substrate. The droplets on the hydrophilic substrate have fewer free NPs in the interfacial region comparing to the droplets on other solid substrates, leaving large areas of open region not covered by NPs. Thus, the evaporation rate enhancement on the hydrophilic substrate is weakest. For θ eq = 90° droplet, the increasing degree of evaporation enhancement during evaporation is attributed to the increasing number density associated with the decreasing interfacial area while maintaining a uniform circumferential distribution. We note that the volume of N = 100 droplet decreases slower than that of N = 50 droplet at early stage of the evaporation in Fig. 4(b). The reason is that N = 100 droplet occasionally deposits more NPs on the substrate than N = 50 droplet. For θ eq = 120° , the NPs are guided by the advection to accumulate on the apex of droplet, and the contact line region is depleted. Therefore, even few particles are left on the hydrophobic substrate. Meanwhile, the advection keeps the number of NPs in the interfacial region while the interfacial area decreases, leading to significant increase of particle density and thereby evaporation enhancement. For

N = 50 , 100 , and 150 droplets, the influence of uniformly-distributed particles ( θ eq = 90° ) turns out to be more pronounced than that of apex-concentrated particles ( θ eq = 120° ) on the evaporation enhancement. We believe the degree of evaporation enhancement is determined by combined effects of particle number density and covered interfacial area.


Page 23 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Figure 6. Representative steady-state flow fields inside the (a-c) evaporating particle-free droplet and (d-f) particle-laden droplet with N = 200 . (a, d) droplets having equilibrium contact angle

θ eq = 60° , (b, e) droplets having equilibrium contact angle θ eq = 90° , and (c, f) droplets having equilibrium contact angle θ eq = 120° . The red (green) line illustrates the outer (inner) boundary of the interfacial region. φ = m0.99φ * is used to define the outer (inner) boundary of the interfacial region. The blue arrows are velocity vectors, whose lengths show their actual magnitudes multiplied by 104 .

We claim that the assembly of interface-bound NPs is dominated by the flow field induced by the evaporation of sessile droplet with moving contact line. Masoud and Felske provide analytical solutions of the flow fields for evaporating sessile droplets with moving contact line, considering non-uniform evaporation flux along the droplet height and isothermal conditions (no thermal Marangoni effect).56 Notably, our model captures the same hydrodynamic effect in driving the bulk flow in evaporating droplets under isothermal conditions. Therefore, we attempt to compare our flow fields with the analytical solutions (see Fig. 6 in Ref. 56). In their results, an evaporating θ eq = 60° droplet with free moving contact line generates a large clockwise circulation in the bulk. However, there are several streamlines that generate from the interface and go out of the droplets in the contact line region. In order to verify this detail, the streamlines in our simulation during the early stage of droplet evaporation are plotted in Fig. S8. From Fig. S8(b), we find a few streamlines going out of the droplet near the contact line region, which is consistent with the analytical flow field. On the other hand, the analytical solution of 23 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 35

θ eq = 120° droplet shows two circulations with opposite directions formed in the bulk. The central circulation is counterclockwise as in Fig. 6(c). However, the small circulation near the contact line is not clear in that figure. We plot the streamlines for the early stage of evaporating

θ eq = 120° droplet in Fig. S9. Both circulations can be clearly observed with the directions consistent with the analytical solution. These detailed flow patterns provide strong evidence that our model captures the correct hydrodynamics and its influence on the particle assembly. We note that the small circulation near the contact line in θ eq = 120° droplet is transient in our simulations, which only occur when droplet evaporation begins. The small circulation diminishes once the evaporation is quasi-steady and the flow field is fully developed as shown in Fig. 6(c). This discrepancy may be attributed to the diffuse interface and the evaporation boundary condition applied in our model.

C. Particle Deposition in Colloidal Droplets During evaporation, the NPs are deposited to the substrate by the wall attraction, particleparticle interaction, and the phenomenological vertical force in the vapor phase. The hybrid LBBD simulation allows one to elucidate the detailed structure of the final particle deposits, as well as the relation between structure and evaporation-driven assembly process. Figure 7(a) shows the two-dimensional structure of the final particle deposit produced by an evaporating particlecovered droplet on the hydrophilic ( θ eq = 60° ) substrate. The spatial distribution of NPs indicates that the majority of the particles are deposited to the two contact line regions of the droplet. To quantify the deposit pattern, a particle distribution histogram along the substrate is calculated. Figure 7(c) shows two peaks at the edges of the deposit with a few particles deposited at the center. This edge deposition is induced by the advection-driven aggregation near the contact line as discussed in the previous section.


Page 25 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Figure 7. Particle depositions for droplet with particles in the bulk ( φ p = φ * ) and at the interface (

φ p = 0 ) on the hydrophilic ( θ eq = 60° ) substrate. (a) Two-dimensional structure of the final deposit of the droplet with particles at the interface. (b) Two-dimensional structure of the final deposit of the droplet with particles in the bulk. (c) Particle distribution histograms of the final deposits from the two types of particle-laden droplets. (d) Two-dimensional radial distribution functions obtained from the particle positions in the final deposits. The inset in (d) shows a schematic of the hexagonal ordering.

To highlight the importance of particle location in the drying droplet on the final deposit, we also simulate an evaporating droplet with particles dispersed in the bulk by assigning the particle order parameter as φ p = φ * . The interfacial region of the droplet is free of particles during evaporation. The droplet exhibits constant equilibrium contact angle as set to the very end 25 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 35

of the evaporation stage through continuously receding contact line. The evaporation-driven flow in the bulk is weaker than that at the interface and the circulation (see Figs. 6(a, d)) traps the NPs. In contrast to the particle-covered droplet, the final deposit of the dried particle-suspension droplet is concentrated at the center of the deposit and forms a dome-like structure as shown in Fig. 7(b), and the corresponding histogram has only one large peak. This deposit pattern closely resembles the structure of printed photonic crystal domes in recent experiments.60 In addition, the deposit of particle-suspension droplet is smaller than the one of particle-covered droplet, which can be attributed to weak contact line pinning induced by the interface-bound NPs. The two droplets with different particle distributions produce rather different deposit patterns upon complete evaporation. From the particle radial distribution functions in Fig. 7(d), we find that both deposits exhibit multi-layer structure with two-dimensional hexagonal order. However, better ordering of NPs is formed from the particle-suspension droplet compared to the particlecovered droplet, which is revealed by the lower heights of the peaks corresponding to nearest, second-nearest, and third-nearest neighbors for the deposit of the particle-covered droplet in Fig. 7(d). The position of the nearest neighbor peak characterizes the effective radius of the NPs and is consistent with re = 1.2 . The theoretical ratio of peak heights for 2D hexagonal close-packed structure is 1:1/ 3 :1/ 2 . Notably, the structure of the obtained deposit is marked by a pronounced first peak and an indiscernible second peak, which are attributed to the additional monolayer structure seen in Figs. 7(a, b). As shown in Figs. S10 and S11, both the particle-suspension and particle-covered droplets yield final particle deposits with only one peak in the histogram on the neutral-wetting ( θ eq = 90° ) and hydrophobic ( θ eq = 120° ) substrates. However, there are discernible differences in the heights and positions of peaks, as well as the out-of-plane shape of the deposits. The peak position of the dried particle-suspension droplet is close to the center position of the initial droplet, which indicates that the NPs dispersed in the bulk do not significantly influence the contact line dynamics. Both contact lines recede simultaneously upon evaporation. In contrast, the peak of the dried particle-covered droplet is shifted to one side of the initial droplet due to asymmetric contact line receding, caused by the weak pinning effect of the interface-bound NPs. Compared to the hydrophilic substrate, the particle deposition region reduces as the equilibrium contact angle increases on the neutral-wetting and hydrophobic substrates, while the peak height 26 ACS Paragon Plus Environment

Page 27 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

increases accordingly. This trend implies that the particle influence on the contact line dynamics of the droplets on the hydrophobic substrate is not as significant as it on the hydrophilic and neutral-wetting substrates, which is consistent with our previous discussion of the flow-induced particle aggregation in evaporation droplets. Moreover, the particle-dispersion droplets consistently produce deposits with a dome shape on the neutral-wetting and hydrophobic substrates, similar to the ones obtained on the hydrophilic substrates. In contrast, the deposits of the particle-cover droplets have comparably irregular shapes, which can generally be considered as multilayer films. Our result confirms the direct connection between the hydrodynamics in the evaporating droplet and the deposit pattern of particle assembly.

IV. Conclusions In this work, we present a free-energy-based D2Q9 lattice Boltzmann method coupled with Brownian dynamics to simulate evaporating particle-laden fluidic systems. The Dirichlet boundary condition is assigned to model evaporation. For an evaporating flat liquid-vapor interface, we find that the NPs at the interface enhance evaporation by tracking the changes in the interfacial position and evaporation flux. In particular, we observe that the evaporation flux increases and the final interfacial position significantly decreases after we add a small amount of NPs to the liquid-vapor interface. However, further increase of NPs shows no obvious enhancement. The evaporating particle-covered droplets with moving contact line are investigated on the solid substrates with different wetting properties (e.g., hydrophilic, neutralwetting and hydrophobic). We also observe that large particle number density at the liquid-vapor interface results in significant enhancement in droplet evaporation. Evaporation rate is also found to be influenced by the particle distribution on the surface, which is further attributed to different wetting properties of the substrate. On the hydrophilic substrate, the NPs at the interface are prone to accumulate near the contact line. On the neutral-wetting substrate, the particles keep a relatively homogenous distribution during evaporation. On the hydrophobic substrate, the NPs prefer to aggregate at the apex of the droplet. We reveal that the non-uniform circumferential distribution of particles is induced by the evaporation-driven advection, and the flow field is dependent on the equilibrium contact angle of droplet. Finally, we investigate the final particle deposits of dried particle-covered droplets and compare them with those obtained from droplets with particles dispersed in the bulk. Different solid substrates are considered and the result also shows the particles location has a significant influence on the deposition. The particle-covered 27 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 35

and particle-suspension droplets lead to distinct deposit patterns, which highlight different dynamics of the interface-bound and bulk-dispersed NPs. This work provides fundamental hydrodynamic insight into the evaporation-induced particle self-assembly under isothermal conditions and presents a useful computational tool for predicting the structure of particle deposit.

Acknowledgment This work was supported by the National Science Foundation under grant No. CMMI-1538090. Generous allocations of computing time were provided by the Thomas J. Watson School of Engineering and Applied Science Computational Data Center at Binghamton University. The authors also think Professors Paul Chiarot and Timothy Singler for helpful conversations.

Supporting Information Section 1: details of the mapping between the dimensionless simulation parameters to physical values. Section 2: detailed characterization of the static colloidal droplets and flat interfaces, as well as the influences of particle-fluid interaction parameters and particle number.


Page 29 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

References (1)

Parviz, B. A.; Ryan, D.; Whitesides, G. M. Using Self-Assembly for the Fabrication of Nano-Scale Electronic and Photonic Devices. IEEE Trans. Adv. Packag. 2003, 26, 233– 241.

(2)

Nie, Z.; Petukhova, A.; Kumacheva, E. Properties and Emerging Applications of SelfAssembled Structures Made from Inorganic Nanoparticles. Nat. Nanotechnol. 2010, 5, 15–25.

(3)

Busseron, E.; Ruff, Y.; Moulin, E.; Giuseppone, N. Supramolecular Self-Assemblies as Functional Nanomaterials. Nanoscale 2013, 5, 7098.

(4)

Hench, L. L.; West, J. K. The Sol-Gel Process. Chem. Rev. 1990, 90, 33–72.

(5)

Derby, B. Inkjet Printing of Functional and Structural Materials: Fluid Property Requirements, Feature Stability, and Resolution. Annu. Rev. Mater. Res. 2010, 40, 395– 414.

(6)

Sefiane, K. Patterns from Drying Drops. Adv. Colloid Interface Sci. 2014, 206, 372–381.

(7)

Routh, A. F. Drying of Thin Colloidal Films. Reports Prog. Phys. 2013, 76, 46603.

(8)

Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R.; Witten, T. A. Capillary Flow as the Cause of Ring Stains from Dried Liquid Drops. Nature 1997, 389, 827–829.

(9)

Han, W.; Lin, Z. Learning from “Coffee Rings”: Ordered Structures Enabled by Controlled Evaporative Self-Assembly. Angew. Chemie Int. Ed. 2012, 51, 1534–1546.

(10)

Shmuylovich, L.; Shen, A. Q.; Stone, H. A. Surface Morphology of Drying Latex Films: Multiple Ring Formation. Langmuir 2002, 18, 3441–3445.

(11)

Nguyen, V. X.; Stebe, K. J. Patterning of Small Particles by a Surfactant-Enhanced Marangoni-Bénard Instability. Phys. Rev. Lett. 2002, 88, 164501.

(12)

Deegan, R. D. Pattern Formation in Drying Drops. Phys. Rev. E 2000, 61, 475–485.

(13)

Hu, H.; Larson, R. G. Marangoni Effect Reverses Coffee-Ring Depositions. J. Phys. Chem. B 2006, 110, 7090–7094.


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(14)

Page 30 of 35

Kim, H.; Boulogne, F.; Um, E.; Jacobi, I.; Button, E.; Stone, H. A. Controlled Uniform Coating from the Interplay of Marangoni Flows and Surface-Adsorbed Macromolecules. Phys. Rev. Lett. 2016, 116, 124501.

(15)

Bhardwaj, R.; Fang, X.; Somasundaran, P.; Attinger, D. Self-Assembly of Colloidal Particles from Evaporating Droplets: Role of DLVO Interactions and Proposition of a Phase Diagram. Langmuir 2010, 26, 7833–7842.

(16)

Still, T.; Yunker, P. J.; Yodh, A. G. Surfactant-Induced Marangoni Eddies Alter the Coffee-Rings of Evaporating Colloidal Drops. Langmuir 2012, 28, 4984–4988.

(17)

Sempels, W.; De Dier, R.; Mizuno, H.; Hofkens, J.; Vermant, J. Auto-Production of Biosurfactants Reverses the Coffee Ring Effect in a Bacterial System. Nat. Commun.

2013, 4, 1757. (18)

Marin, A.; Liepelt, R.; Rossi, M.; Kähler, C. J. Surfactant-Driven Flow Transitions in Evaporating Droplets. Soft Matter 2016, 12, 1593–1600.

(19)

Lin, X. M.; Jaeger, H. M.; Sorensen, C. M.; Klabunde, K. J. Formation of Long-RangeOrdered Nanocrystal Superlattices on Silicon Nitride Substrates. J. Phys. Chem. B 2001, 105, 3353–3357.

(20)

Narayanan, S.; Wang, J.; Lin, X.-M. Dynamical Self-Assembly of Nanocrystal Superlattices during Colloidal Droplet Evaporation by in Situ Small Angle X-Ray Scattering. Phys. Rev. Lett. 2004, 93, 135503.

(21)

Bigioni, T. P.; Lin, X.-M.; Nguyen, T. T.; Corwin, E. I.; Witten, T. A.; Jaeger, H. M. Kinetically Driven Self Assembly of Highly Ordered Nanoparticle Monolayers. Nat. Mater. 2006, 5, 265–270.

(22)

Yunker, P. J.; Still, T.; Lohr, M. a; Yodh, a G. Suppression of the Coffee-Ring Effect by Shape-Dependent Capillary Interactions. Nature 2011, 476, 308–311.

(23)

Jafari Kang, S.; Vandadi, V.; Felske, J. D.; Masoud, H. Alternative Mechanism for Coffee-Ring Deposition Based on Active Role of Free Surface. Phys. Rev. E 2016, 94, 63104.

(24)

Yunker, P. J.; Lohr, M. A.; Still, T.; Borodin, A.; Durian, D. J.; Yodh, A. G. Effects of 30 ACS Paragon Plus Environment

Page 31 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Particle Shape on Growth Dynamics at Edges of Evaporating Drops of Colloidal Suspensions. Phys. Rev. Lett. 2013, 110. (25)

Anyfantakis, M.; Geng, Z.; Morel, M.; Rudiuk, S.; Baigl, D. Modulation of the CoffeeRing Effect in Particle/Surfactant Mixtures: The Importance of Particle–Interface Interactions. Langmuir 2015, 31, 4113–4120.

(26)

Li, Y.; Yang, Q.; Li, M.; Song, Y. Rate-Dependent Interface Capture beyond the CoffeeRing Effect. Sci. Rep. 2016, 6, 24628.

(27)

Osher, S.; Sethian, J. A. Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations. J. Comput. Phys. 1988, 79, 12–49.

(28)

Hirt, C. .; Nichols, B. . Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries. J. Comput. Phys. 1981, 39, 201–225.

(29)

Pilliod, J. E.; Puckett, E. G. Second-Order Accurate Volume-of-Fluid Algorithms for Tracking Material Interfaces. J. Comput. Phys. 2004, 199, 465–502.

(30)

Swift, M.; Orlandini, E.; Osborn, W.; Yeomans, J. Lattice Boltzmann Simulations of Liquid-Gas and Binary Fluid Systems. Phys. Rev. E 1996, 54, 5041–5052.

(31)

Inamuro, T.; Ogata, T.; Tajima, S.; Konishi, N. A Lattice Boltzmann Method for Incompressible Two-Phase Flows with Large Density Differences. J. Comput. Phys. 2004, 198, 628–644.

(32)

Lee, T.; Lin, C. L. A Stable Discretization of the Lattice Boltzmann Equation for Simulation of Incompressible Two-Phase Flows at High Density Ratio. J. Comput. Phys.

2005, 206, 16–47. (33)

Zheng, H. W.; Shu, C.; Chew, Y. T. A Lattice Boltzmann Model for Multiphase Flows with Large Density Ratio. J. Comput. Phys. 2006, 218, 353–371.

(34)

Lee, T.; Liu, L. Lattice Boltzmann Simulations of Micron-Scale Drop Impact on Dry Surfaces. J. Comput. Phys. 2010, 229, 8045–8063.

(35)

Briant, A. J.; Wagner, A. J.; Yeomans, J. M. Lattice Boltzmann Simulations of Contact Line Motion. I. Liquid-Gas Systems. Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys.

2004, 69, 1–14. 31 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(36)

Page 32 of 35

Dupuis, A.; Yeomans, J. Lattice Boltzmann Modelling of Droplets on Chemically Heterogeneous Surfaces. Futur. Gener. Comput. Syst. 2004.

(37)

Cates, M. E.; Stratford, K.; Adhikari, R.; Stansell, P.; Desplat, J.-C.; Pagonabarraga, I.; Wagner, a. J. Simulating Colloid Hydrodynamics with Lattice Boltzmann. J. Phys. Condens. matter 2004, 3903, 14.

(38)

Joshi, A. S.; Sun, Y. Multiphase Lattice Boltzmann Method for Particle Suspensions. Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys. 2009, 79, 1–16.

(39)

Joshi, A. S.; Sun, Y. Wetting Dynamics and Particle Deposition for an Evaporating Colloidal Drop: A Lattice Boltzmann Study. Phys. Rev. E 2010, 82, 41401.

(40)

Connington, K. W.; Lee, T.; Morris, J. F. Interaction of Fluid Interfaces with Immersed Solid Particles Using the Lattice Boltzmann Method for Liquid-Gas-Particle Systems. J. Comput. Phys. 2015, 283, 453–477.

(41)

Connington, K. W.; Miskin, M. Z.; Lee, T.; Jaeger, H. M.; Morris, J. F. Lattice Boltzmann Simulations of Particle-Laden Liquid Bridges: Effects of Volume Fraction and Wettability. Int. J. Multiph. Flow 2015, 76, 32–46.

(42)

Huang, H.; Zheng, H.; Lu, X.; Shu, C. An Evaluation of a 3D Free-Energy-Based Lattice Boltzmann Model for Multiphase Flows with Large Density Ratio. Int. J. Numer. Methods Fluids 2009, 63, 1193–1207.

(43)

Succi, S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond; Oxford University Press, 2001.

(44)

Swift, M. R.; Osborn, W. R.; Yeomans, J. M. Lattice Boltzmann Simulation of Nonideal Fluids. Phys. Rev. Lett. 1995, 75, 830–833.

(45)

Ma, Y.; Bhattacharya, A.; Kuksenok, O.; Perchak, D.; Balazs, A. C. Modeling the Transport of Nanoparticle-Filled Binary Fluids through Micropores. Langmuir 2012, 28, 11410–11421.

(46)

Verberg, R.; Yeomans, J. M.; Balazs, A. C. Modeling the Flow of Fluid/particle Mixtures in Microchannels: Encapsulating Nanoparticles within Monodisperse Droplets. J. Chem. Phys. 2005, 123, 224706. 32 ACS Paragon Plus Environment

Page 33 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

(47)

Liu, Y.; Bhattacharya, A.; Kuksenok, O.; He, X.; Aizenberg, M.; Aizenberg, J.; Balazs, A. C. Computational Modeling of Oscillating Fins That “catch and Release” Targeted Nanoparticles in Bilayer Flows. Soft Matter 2016, 12, 1374–1384.

(48)

Briant, A. J.; Papatzacos, P.; Yeomans, J. M. Lattice Boltzmann Simulations of Contact Line Motion in a Liquid-Gas System. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 2002, 360, 485–495.

(49)

Fakhari, A.; Rahimian, M. H. Phase-Field Modeling by the Method of Lattice Boltzmann Equations. Phys. Rev. E 2010, 81.

(50)

Shao, J. Y.; Shu, C.; Huang, H. B.; Chew, Y. T. Free-Energy-Based Lattice Boltzmann Model for the Simulation of Multiphase Flows with Density Contrast. Phys. Rev. E 2014, 89, 33309.

(51)

Bhattacharya, A.; Balazs, A. C. Stiffness-Modulated Motion of Soft Microscopic Particles over Active Adhesive Cilia. Soft Matter 2013, 9, 3945.

(52)

Fakhari, A.; Geier, M.; Lee, T. A Mass-Conserving Lattice Boltzmann Method with Dynamic Grid Refinement for Immiscible Two-Phase Flows. J. Comput. Phys. 2016, 315, 434–457.

(53)

Shao, J. Y.; Shu, C. A Hybrid Phase Field Multiple Relaxation Time Lattice Boltzmann Method for the Incompressible Multiphase Flow with Large Density Contrast. Int. J. Numer. Methods Fluids 2015, 77, 526–543.

(54)

Zu, Y. Q.; He, S. Phase-Field-Based Lattice Boltzmann Model for Incompressible Binary Fluid Systems with Density and Viscosity Contrasts. Phys. Rev. E 2013, 87, 43301.

(55)

Ledesma-Aguilar, R.; Vella, D.; Yeomans, J. M. Lattice-Boltzmann Simulations of Droplet Evaporation. Soft Matter 2014, 10, 8267–8275.

(56)

Masoud, H.; Felske, J. D. Analytical Solution for Stokes Flow inside an Evaporating Sessile Drop: Spherical and Cylindrical Cap Shapes. Phys. Fluids 2009, 21, 42102.

(57)

Yong, X.; Qin, S.; Singler, T. J. Nanoparticle-Mediated Evaporation at Liquid–vapor Interfaces. Extrem. Mech. Lett. 2016, 7, 90–103.

(58)

Cheng, S.; Grest, G. S. Molecular Dynamics Simulations of Evaporation-Induced 33 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 35

Nanoparticle Assembly. J. Chem. Phys. 2013, 138, 64701. (59)

Barnes, G. T. The Effects of Monolayers on the Evaporation of Liquids. Adv. Colloid Interface Sci. 1986, 25, 89–200.

(60)

Kuang, M.; Wang, J.; Bao, B.; Li, F.; Wang, L.; Jiang, L.; Song, Y. Inkjet Printing Patterned Photonic Crystal Domes for Wide Viewing-Angle Displays by Controlling the Sliding Three Phase Contact Line. Adv. Opt. Mater. 2014, 2, 34–38.


Page 35 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

TOC Graphic

ACS Paragon Plus Environment

Modeling Evaporation and Particle Assembly in ... - ACS Publications

Recommend Documents