Pinning–Depinning Mechanisms of the Contact Line during

Jun 11, 2018 - E-mail: [email protected]. ... Dynamics of Contact Line Pinning and Depinning of Droplets Evaporating on Microribs. Langmuir. Mazloom...
0 downloads 0 Views 1MB Size
Subscriber access provided by Kaohsiung Medical University

New Concepts at the Interface: Novel Viewpoints and Interpretations, Theory and Computations

Pinning-depinning mechanisms of the contact line during evaporation of micro-droplets on rough surfaces: A lattice Boltzmann simulation Wuzhi Yuan, and Lizhi Zhang Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.8b00857 • Publication Date (Web): 11 Jun 2018 Downloaded from http://pubs.acs.org on June 12, 2018

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 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 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.

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 31 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

Pinning-depinning mechanisms of the contact line during evaporation of micro-droplets on rough surfaces: A lattice Boltzmann simulation Wu-Zhi Yuan1, Li-Zhi Zhang1,2∗ 1. Key Laboratory of Enhanced Heat Transfer and Energy Conservation of Education Ministry, School of Chemistry and Chemical Engineering, South China University of Technology, Guangzhou 510640, China 2. State Key Laboratory of Subtropical Building Science, South China University of Technology, Guangzhou 510640, China ABSTRACT In this study, pinning and depinning of the contact line during droplet evaporation on the rough surfaces with randomly distributed structures is theoretically analyzed and numerically investigated. A Fast Fourier Transformation (FFT) method is used to generate the rough surfaces, whose skewness (Sk), kurtosis (K) and root mean square (Rq) are obtained from real surfaces. A thermal multiphase LB model is proposed to simulate the isothermal pinning and depinning processes. The evaporation processes are recorded with the variations in contact angle, contact radius and drop shape. It is found that the drops sitting on rough surfaces show different behavior from those on smoother surfaces. The former shows a pinned contact line during almost the whole lifetime. By contrast, the latter experiences a stick-slip-jump behavior until the drop disappears. At mesoscopic scale, the pinning of the contact line is actually a slow motion rather than a complete immobilization at the sharp edges. The dynamic equilibrium is achieved by the self-adjustment of the contact line according to each edge. ∗

Author for correspondence. Tel/fax: 86-20-87114264; Email: [email protected]

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

INTRODUCTION Droplet evaporation is not only a phenomenon in nature, but also a process having various practical applications, such as spraying cooling,1, 2 inkjet printing,3-5 surface coatings6, 7

and chip manufacturing.8,

9

Understanding the fundamental mechanisms of droplet

evaporation on solid surfaces is of great importance in these applications. It has been well accepted that there exist three evaporation modes for sessile droplets: the constant contact radius (CCR) mode, the constant contact angle (CCA) mode and the mixed mode.10, 11 In the CCR mode, the contact line is pinned on the solid substrate. Droplets evaporate with constant contact radius and diminishing contact angle. In the CCA mode, the length of the contact line decreases while the contact angle remains constant. In the mixed mode, both the contact line length and contact angle change during evaporation. Numerous studies have been performed to investigate droplet evaporation on various surfaces.12-16 However, up until now, some aspects of droplet evaporation are still not clear, such as the pinning-depinning phenomenon of the contact line during droplet evaporation on heterogeneous surfaces. The complexity of contact line dynamics during droplet evaporation arises from the surface morphology, chemical heterogeneity, and the interfacial wetting states. Chen et al.17 studied the contact line dynamics in terms of force balance between pinning force Fp and depinning forces Fd during droplet evaporation on micro-textured surfaces. They claimed that pinning is caused by Fp>Fd. Orejon et al.18 investigated the motion of the contact line of evaporating droplets on surfaces with varying hydrophobicity. They proposed that the existence of an energy barrier between CCR and CCA state is the origin of the pinning of the

ACS Paragon Plus Environment

Page 2 of 31

Page 3 of 31 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

contact line. However, these two explanations originated from a macro-scale experiment. The motion of the contact line along the micro-structure is hard to observe experimentally with the usual optical methods employed. As an alternative, a meso-scale numerical simulation can capture the dynamics of the contact line and could provide a more detailed information. For nanodroplets composed of hundreds or thousands of liquid molecules, molecular dynamics (MD) simulations have been used to study their contact line dynamics on heterogeneous substrates.19-22 Their simulations reproduce the pinning-depinning phenomenon during droplet evaporation and reveal a molecular origin of the pinning force. However, for microscopic droplets studied here, their size is beyond the scope of MD simulation. In the present work, a hybrid thermal multiphase lattice Boltzmann method (LBM) is employed to study the droplet evaporation on rough surfaces with randomly distributed structures. The LBM can be considered as a mesoscopic method, which provides a viable way to include the microscopic dynamics at the fluid-fluid and fluid-solid interfaces. Therefore, we can reveal the micro-scale effects on droplet evaporation on rough surface.23, 24 In literature, there have been some LB studies on droplet evaporation, mainly on flat and/or chemically heterogeneous surfaces. Li et al.25 studied the mechanisms of droplet evaporation on chemically patterned surfaces and successfully captured the pinning-depinning behavior of the contact line. Further, Rodrigo et al.26 studied the evaporation of planar films and 3D sessile droplets with LBM. These studies are very instructive. However, in all these studies, the surfaces are assumed to be smooth, which is obviously invalid for real surfaces. That’s our purpose. To mimic real surfaces on which the rough structures are randomly

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 31

distributed, we firstly reconstruct the randomly distributed rough surfaces by the Fast Fourier Transformation (FFT).27-29 Then, droplet evaporation is investigated on these rough surfaces. To our best knowledge, this is the first time LBM is used to study the droplet evaporation on the randomly distributed real rough surfaces. The purpose of the present work is to understand the dynamics of the contact line during evaporation of droplets on rough surfaces by LB simulations. A theoretical analysis, with a thermal multiphase Lattice Boltzmann model, is also proposed to interpret the pinning and depinning mechanisms. The rest of paper is organized as follows. In Section 2, a brief introduction about the numerical model and the reconstruction of rough surface are presented. The main results are discussed in Section 3. Finally, a conclusion will be made in section 4. NUMERICAL MODEL AND SURFACE RECONSTRUCTION Lattice Boltzmann Method Based on the simple and famous Bhatnagar-Gross-Krook approximation,30 the LB equation, which governs the evolution of the density distribution function, can be written as

f i (x + e i δt , t + δt ) − f i (x, t ) = −

δt [ f i (x, t ) − f i eq (x, t )] + S i τ

(1)

where fi(x,t) is the density distribution function at the lattice site x and time t, fieq is the equilibrium distribution function, τ is the dimensionless relaxation time and Si represents a general source term. For the D2Q9 model with nine velocity directions at a given point in two-dimensional space, the discrete velocities ei is given by 0 1 0 − 1 0 1 − 1 − 1 1  [e 0 , e1 , e 2 , e 3 , e 4 , e 5 , e 6 , e 7 , e 8 ] =   0 0 1 0 − 1 1 1 − 1 − 1

The equilibrium distribution function fieq for the D2Q9 model are of the form

ACS Paragon Plus Environment

(2)

Page 5 of 31 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

 e ⋅ u eq (e ⋅ u eq ) (u eq ) 2  f i eq (x, t ) = ρwi 1 + i 2 + i 4 −  cs 2cs 2cs2  

(3)

where the weight factors wi are given by w0=4/9, w1-4=1/9 and w5-8=1/36. The fluid density ρ and velocity u can be obtained from the first and second moments of the density distribution functions

ρ = ∑ f i , ρ u = ∑ei fi

(4)

i

i

The viscosity in the LBM is given by

υ = ( τ − 0.5)cs2δt

(5)

where cs = c/ 3 is the lattice sound speed, and c=δx/δt is the lattice speed. δx and δt are the lattice spacing and time step, respectively. In order to calculate the source term in eq 1, the exact difference method (EDM) is adopted. The EDM is directly derived from the Boltzmann equation, which has been applied in a variety of fields with great success.31, 32 The source term is calculated by

S i (x, t ) = f i eq ( ρ , u + ∆ u) − f i eq ( ρ , u)

(6)

where ∆u=Fδt/ρ being the velocity change due to the force term F during time step δt. The force F acting on the fluid consists of fluid-fluid force, fluid-solid force. The fluid-fluid force Fl 8

Fl (x,t ) = −Gψ(x,t )∑ wiψ(x+ e i , t ) e i

(7)

i =1

and the fluid-solid force Fs 8

Fs (x,t ) = −Gadsψ(x,t )∑ wi s(x+ ei , t ) ei

(8)

i =1

where G and Gads are parameters that control the strength of inter-particles and the strength of interaction between the fluid and the solid wall. Various wettability can be achieved by regulating the parameter Gads. s is an indicator function equals 1 or 0 for solid nodes and fluid

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 6 of 31

nodes, respectively. ψ(x, t), the effective mass, is a function of the local density and it depends on the equation of state (EOS). It can be expressed as

ψ (x, t ) =

2( p − ρc s2 ) /( Gc s2 )

(9)

Based on Yuan and Schaefer’s investigation,33 the Peng-Robinson (P-R) EOS is adopted, which provides the maximum density ratio and maintains small spurious currents around the interface. The P-R EOS can be expressed as

p=

ρRT aρ 2 κ(T ) − 1 − bρ 1 + 2bρ − b2 ρ 2

(10)

where

κ (T ) = [1 + (0.37464 + 1.5422ω − 0.26992ω2 )(1 − T / Tc )]2

(11)

2 2 ω is the acentric factor, the attraction parameter a=0.45724R Tc /Pc and the repulsion

parameter b=0.0778TRc/Pc. Tc and Pc are the critical temperature and pressure, respectively. In the present study, we set a=2/49, b=2/21 and R=1.33 Thus ψ can be obtained by ψ =

6 ρRT aρ 2 κ (T ) ρ ( − − ) G 1 − bρ 1 + 2 bρ − b 2 ρ 2 3

(12)

Energy Equation The lattice Boltzmann model is used to simulate the density and velocity fields. In order to simulation the droplet evaporation, a temperature solver is also needed to solve the energy equation. Neglecting the viscous heat dissipation, it is written as34, 35

ρc v ( ∂ tT + u⋅ ∇T ) = ∇ ⋅ ( λ∇T ) − T (

∂p EOS )ρ ∇ ⋅ u ∂T

(13)

where λ is the thermal conductivity and cv is the specific heat at a constant volume. Generally, two schemes are used to solve the energy equation in the thermal LBM: the double distribution function (DDF) approach and the finite-difference scheme. In DDF

ACS Paragon Plus Environment

Page 7 of 31 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

approach, following the same procedure, another distribution function is used to solve the energy equation. Based on the evaluation by Li et al.,36 we adopt the finite-difference scheme in this study, which is simple and accurate compared to the DDF approach. To solve eq 13, here the second-order Runge-Kutta scheme is adopted for the time discretization, i.e.,

T (x, t + δt ) = T (x, t ) +

δt ( h1 + h2 ) 2

(14)

where h1and h2 are given by

h1 = M [T (x, t )]

h2 = M [T (x,t ) +

δt h1 ] 2

(15)

where M(T) represents the right-hand side of eq 13. As for the spatial discretization, the first and second order derivative terms in eq 13 can be calculated by the following second-order isotropic difference schemes, respectively.

∂ i φ(x) ≈

∇ 2φ(x) ≈

1 cδ

∑ w φ(x+ e δ ) e

2 cδ

∑ w [φ(x+ e δ ) − φ(x)]

2 s t

2 2 s t

i

i

t

i

(16)

i

i

i

t

(17)

i

where φ denotes a physical variable. Obviously, using a non-ideal equation of state (EOS), the pseudopotential ψ in eq 12 will be linked to the temperature field. Then the liquid-vapor phase change can be driven by the temperature field through the equation of state (EOS). Boundary Conditions In the present study, for fluid flow, three kinds of boundary conditions are adopted: no-slip boundary conditions at the solid wall, the constant pressure boundary condition at the top boundary and periodic boundary conditions at the left and right boundaries. In the LB

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

framework, for no slip boundary condition, the half way bounce-back scheme is employed;37 for the constant pressure boundary condition, the non-equilibrium extrapolation method proposed by Guo et al.38 is adopted for its good accuracy; and for periodic boundary conditions, the unknown distribution function at one boundary is set to that at the periodic boundary. For phase change, three types of boundary conditions are used: the isothermal wall boundary condition at the fluid-solid interface, the open boundary condition at the top boundary, and periodic boundary conditions at the left and right boundaries. A flow chart of the solution procedure is shown in Figure S1. The detailed solution procedure can be found in Supporting Information. Reconstruction of the Randomly Distributed Rough Surface In order to investigate the effect of rough structures on droplet evaporation, the rough surfaces should be reconstructed before the simulation. There are usually two ways to obtain randomly distributed rough surfaces. One is obtained by experimental measurement, such as AFM, and the other is based on the computer reconstruction. Because the measurement of surface roughness is a time-consuming process and it can only measure some limited sections at a time. The roughness data sets are quite finite. The computer reconstruction is relatively simple and rough surfaces with arbitrary size and roughness can be generated. As a result, based on the limited measured surface data, a statistical method is applied to reconstruct the randomly distributed rough surfaces. Based on previous studies,27-29 a fast and convenient tool in generating rough surfaces, the Fast Fourier Transformation (FFT), is used to generate non-Gaussian randomly distributed rough surfaces with the Sk, K and Rq obtained from the

ACS Paragon Plus Environment

Page 8 of 31

Page 9 of 31 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

real surfaces. The reconstruction of rough surfaces mainly consists of two steps. First, the structural parameters of the rough surfaces are obtained by AFM. Second, based on the measured parameters, the FFT is applied to reconstruct the randomly distributed rough surfaces. The detail of the reconstruction can be found in our previous investigation.39, 40 Four kinds of typical rough surfaces are reconstructed and the rough parameters are listed in the Table 1. Table 1. Parameters for the Rough Surfaces. Surface

K

Sk

Rq (µm)

S1

3.0

-1.2

0.3

S2

3.0

1.2

0.3

S3

3.0

1.2

0.6

S4

5.0

-1.2

0.3

NUMERICAL RESULTS AND DISCUSSION Validation of the code Prior to solving the problem, two benchmarks are adopted to validate the capability of the thermal multiphase LB model for simulating liquid-vapor phase change. Firstly, a stationary non-evaporating droplet is simulated, which is located at the center of a fully periodic domain with mesh size 201× 201. A typical steady-state density contour obtained from simulation is shown in Figure 1(a), in which the reduced temperature is T=0.9Tc and droplet radius is r0=30. To check the coexistence curves, series of simulations are carried out by changing the reduced temperatures T but fixing the r0. The LB simulated coexistence curves for the P-R EOS are in good agreements with the theoretical ones as shown in Figure 1(b), where the discrepancy is due to the thermodynamic inconsistency of the Shan-Chen

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

pseudo-potential model.

Figure 1. (a) The steady-state density contours for T=0.9Tc. (b) Comparison of coexistence curves obtained from LB simulations with the theoretical one predicted by Maxwell equal-area construction.41 Another test is the well-known D2 law, 42, 43 which reads D2(t)/D02=1-kt, where D0 is the initial diameter and D(t) is the droplet diameter at t during evaporation. It is on the basis of the constant temperature of droplet and the constant ambient pressure during the evaporation. In order to keep the vapor pressure constant, a constant pressure boundary condition is implemented. The other conditions are the same as the first validation test. For the temperature, after reached the equilibrium state for the isothermal system T=0.9Tc, the temperature at the boundaries is set to T=0.905Tc. Owing to the temperature difference between the boundaries and droplet, the droplet begins to evaporate. The kinematic viscosity is taken as ν=0.1 in the whole computational domain. The specific heat at constant volume is set to cv=5. Two cases are simulated: Case A with λ=1/3 and Case B with λ=2/3.

ACS Paragon Plus Environment

Page 10 of 31

Page 11 of 31 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. Simulation results of a circular static liquid droplet with evaporation. (a) Temperature distributions and velocity vectors of the evaporating droplet. The black circle represents the liquid-vapor interface. (b) Variation of the square of the non-dimensional droplet diameter with time. The velocity vectors around the droplet and the temperature distributions are shown in Figure 2(a) at t=80000δt for Case A. The radial direction of the velocity vectors normal to the interface are represented as expected. Figure 2(b) presents the time evolution of the normalized squared diameter of the droplet, where D0 is the droplet diameter at t=0. The linear relationship between (D/D0)2 and t can be clearly observed. Droplet Evaporation on the Randomly Distributed Rough Surface Now we turn our attention to the dynamics of evaporating droplets on rough surfaces. Our simulations are carried out in a two-dimensional rectangular domain, as shown in Figure 3. The simulation domain is nx × ny = 1000×1000 (LB units), which are proved to be appropriate after grid dependency tests with three sets of meshes of 800×800, 1000×1000 and 1200×1200. Initially, a droplet with an initial radius r0=150 (LB units) is placed on the rough surfaces re-constructed and is equilibrated for 1.0×105 δt to reach its equilibrium state. After that, the temperature solver is added to model the process of droplet evaporation. In the

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

following simulations, the parameters in lattice units are fixed at: density of liquid ρl=5.9067, density of gas ρv=0.529, which are set as the theoretical densities predicted by the Maxwell equal-area construction41 at the corresponding temperature T=0.9Tc, and the surface tension

γ=0.0237. The intrinsic contact angle is 115° (Gads=-2.0 in eq 8). The temperature of the solid wall is set to Tw=0.905Tc. The kinematic viscosity is taken as ν=0.17 and the thermal diffusivity χ=λ/(ρcv), where λ is the thermal conductivity is set to 0.25 with cv = 2.5.

Figure 3. Schematic diagram of the evaporation of droplet on the randomly distributed rough surface. The red, green and blue colors represent rough surface, liquid droplet and saturating vapor, respectively. The radius of droplet is 150 lattice units. The amplified grid of local grooves is placed in the left corner. The Dynamic Behavior of Evaporating Droplet on Rough Surfaces We firstly investigated the dynamic behaviors of the droplet evaporation on the surface S1. The snapshots of droplet during evaporation are presented in Figure 4. The insets show the motion of the contact line on the left side of the droplet. In the meantime, the variation of the contact angle and the non-dimensional contact radius with respect to time is shown in Figure 5. Due to the heterogeneity, the variation of contact angle is not exactly identical on

ACS Paragon Plus Environment

Page 12 of 31

Page 13 of 31 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

the both sides of the droplet. For convenience, only the contact angle of the droplet on the left side is discussed here. It can be seen from Figure 4(a) that before adding the temperature solver, a droplet with the contact angle of about 126° was formed on the rough surface. After that, the droplet starts to evaporate, resulting from the temperature difference between the droplet and the solid surface. The droplet initially evaporates in the CCR mode. The contact angle decreases gradually while the contact radius remains roughly constant (see Figure 5, I). At a certain time about 1.6×105δt, the contact angle decreases to the equilibrium angle θe =113° (a slightly smaller than that of 115° in isothermal condition due to the Marangoni effect25) and the evaporation shifts to the CCA mode (see Figure 5, II). The contact angle keeps constant while the contact line continues to recede along the surface until it encounters a sharp edge (see Figure 4, b-d), where the contact line is pinned and the droplet evaporates in CCR again (see Figure 5, III). Then, because of further evaporation, the contact angle θ gradually deviates from the equilibrium angle (θθe) gradually approaches

θe,

resulting

in

|cosθH−cosθe|

increases

and

|cosθI−cosθe|

decreases.

Correspondingly, the unbalanced Young’s force FH-y (per unit length) will increase while the one FI will decrease. To achieve a dynamic equilibrium, the length lH gradually decreases, whereas the length lI increases. Eventually, when the contact angle θI reduces to the equilibrium contact angle θe, the length lH reduces to zero and the contact line is totally occupied by the inclined edge BC, marked by the bold red line in Figure 8 (c). The depinning of the contact line occurs and it will slip along the inclined edge from then on.

ACS Paragon Plus Environment

Page 19 of 31 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 8. Droplet evaporation on the surface S1. Local details of the liquid-vapor interfaces at the sharp edge and the density profiles along the solid surface between x=375 and x=405. (a) t=2.5×105δt (b) t=3.0×105δt (c) t=3.2×105δt. The dash line represents the position of sharp edge, at x=391. (III) The origin of jump of the contact line at the junction When the contact angle θI decreases to the equilibrium angle θe at t≈3.2×105δt, the depinning phenomenon occurs. Then, the contact line continuously slip along the inclined edge BC until it reaches the junction C of the two inclined edges, where the contact line quickly crosses the inclined edge CD and arrives at the horizontal edge DE at t≈3.3×105δt, as shown in Figure 7(b). At the junction, the contact line suddenly slides a small distance from the inclined edge BC to the inclined edge CD. However, the contact angle θ changes dramatically from the equilibrium angle θe with the inclined edge BC to a much smaller contact angle (θe-α-β) with the inclined edge CD, as shown in Figure 7(b). Due to the significant deviation from the equilibrium contact angle θe at the junction C, there exists a large unbalanced Young’s force along the inclined edge CD, resulting in the sudden slide of

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

the contact line. Effects of Surface Morphology on the Motion of Contact Line The numerical results of droplet evaporation on the surface S2 is displayed in Figure 9. In the case of the smoother surface S1, the pinning of the contact line is initially observed. Subsequently, the contact line starts to slip along the edge. On the contrary, for drops sitting on the rougher surface S2, the pinning force is so strong that the unbalanced Young’s forces cannot overcome it. Therefore, the droplet evaporates in the CCR mode for almost the whole process of the drop evaporation (see Figure 9(a)). The contact angle decreases gradually while the contact radius remains nearly constant until the drop disappears (see Figure 9(b)). The reason for the different dynamics of contact line on S1 and S2 is the surface morphology, especially the local structure in the vicinity of the contact line. Compared to surface S1, the surface S2 is sharp, featured by a deeper groove and a larger inclined angle. There is a larger gap between the initial contact angle θ and the equilibrium contact angle θe. According to the above analysis, more volume is needed to be evaporated before the contact angle θ decreases to the equilibrium contact angle θe, for the droplet to slip along the inclined edge. Furthermore, due to the deeper groove, the duration of the slip is also prolonged on the surface S2. At last, the size of the droplet is larger compared to the groove and the distance that the droplet slips down the first groove is negligible. All these reasons contribute to the almost immobilization of the evaporating droplet. At t=3.5×105δt, the area of the droplet is reduced by less than 10% upon the initial droplet, but the contact line movement is negligible (see Figure 9(a)).

ACS Paragon Plus Environment

Page 20 of 31

Page 21 of 31 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 9. Droplet evaporation on the rough surface with Sk=1.2, K=3.0 and Rq=0.6. (a) Snapshots of the droplet during evaporation. (b) Time-dependence of contact angle and non-dimensional contact radius. According to above investigations, the surface morphology plays a key role in the dynamics of droplet evaporation. After the initial CCR mode, droplet evaporates in the CCA, CCR and mixed mode randomly, depending on the surface morphology. The duration of the initial CCR mode significantly varies with the local structure in the vicinity of the contact line, such as the height and the slope of the groove. In general, for the randomly distributed rough surface, the surface morphology is not exactly the same even if the rough parameters remain constant. Therefore, in order to gain a better understanding of the effects of surface morphology, twenty surfaces with the same rough parameters are reconstructed for each surface. What’s more, a new index tr is proposed to characterize the droplet evaporation. The tr refers to the fraction of the duration of the first depinning of the contact line in the whole

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

evaporation process. According to tr, the motion of the contact line can be classified into four status: easy to move, 0%~20%; hard to move, 20%~50%; harder to move, 50%~80%; almost immobilization, 80%~100%, which are marked with status 1, 2, 3, 4, respectively.

Figure 10. Distribution of tr for microdroplets on randomly distributed rough surfaces. The tr is divided into four different areas. Figure 10 compares the distribution of tr on each surface in terms of a percentile in the whole lifetime of evaporation. For the surface S1, the surface is much smoother and there are many flat edges on the top of the surface. The tr within 40% accounts for less than 20% of the whole evaporation, indicating that the contact line is easier to move on this surface. When the Sk increases from -1.2 to 1.2, obvious structural changes can be found from the inserted small pictures in Figure 10, S1and S2. With the skewness increasing, there are more sharp edges on

ACS Paragon Plus Environment

Page 22 of 31

Page 23 of 31 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

the top of the rough surfaces, which greatly limit the motion of the contact line. Therefore, on the surface S2, most of the droplets are in Status 4 (see Figure 10, S2). Furthermore, it can be seen from Figure 10 S1 and S3 that with kurtosis increasing, the rough surface has an increasing tendency to form deep, sharp valleys or high peaks. Similar to S2, the increasing sharp edge strongly obstructs the movement of the contact line of the droplet, resulting in decreasing Status 1. At last, the effects of Rq are also investigated. As the Rq increases, the fluctuation of surface roughness becomes greater, thereby forming deeper grooves and larger inclined angles (see Figure 10, S2 and S4). The Status 4 increases from 70% on the surface S2 to 85% on the surface S4. On the rougher surfaces, the droplet prefers to evaporate in the CCR mode, which is attributed to the larger pinning force. From the above analysis, we conclude that the mode of droplet evaporation can be controlled by properly designing the surface roughness, which is of great importance in some applications. CONCLUSION In this work, numerical simulations have been carried out to investigate the droplet evaporation on rough surface with randomly distributed structures. The numerical model is based on the lattice Boltzmann method, which consists of a pseudopotential multiphase lattice Boltzmann model for simulating the density and velocity fields and a finite-difference solver to account for the temperature field and evaporation process. The liquid-vapor phase change is driven by the temperature field via a non-ideal equation of state. Following conclusions are summarized from the study: (1) When a droplet evaporates on rough surfaces, the contact line can be defined as a region with a finite width rather than a strictly one-dimensional object. The pinning of the contact line

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

at the sharp edge is a slow motion rather than being completely immobilized. The dynamic equilibrium is achieved by the self-adjustment of the contact line according to the sharp edges. (2) The numerical results clearly show that after the initial CCR mode, the CCA, the CCR and the mixed modes occur randomly on the rough surfaces depending on the surface morphology. The duration of the initial CCR mode is determined by the local structure in the vicinity of the contact line, such as the height and the slope of the groove. As the height and the slope of the groove increase, the pinning of contact line becomes increasingly strong and thus the duration of the initial CCR mode is prolonged. (3) Even though the modes of droplet evaporation on randomly distributed rough surfaces are random, droplets are likely to evaporate under a specific mode on a properly designed surface. More specifically, droplet on a rough surface with larger Sk, Rq and K are more likely to be in CCR mode during almost the whole process of evaporation. In contrast, on the rough surface with smaller Sk, Rq and K, the contact line is easy to be depinned and the droplet tends to be CCA after a short period of CCR. The obtained new findings are important for understanding the effects of physical roughness on the evaporation of a sessile droplet. It will offer essential guidelines for the design of superhydrophobic surfaces.

ASSOCIATED CONTENT Supporting Information The solution procedure of droplet evaporation. The forces analysis at the sharp edge.

ACS Paragon Plus Environment

Page 24 of 31

Page 25 of 31 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

This information is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION Corresponding Author *E-mail: [email protected] Notes The authors declare no competing financial interest Acknowledgements

The project is supported by the National Key Research and Development Program (No. 2016YFB0901404), and the National Science Fund for Distinguished Young Scholars (No. 51425601). It is also supported by Science and Technology Planning Project of Guangdong Province: Guangdong-Hong Kong Technology Cooperation Funding

Scheme (TCFS), No.2017B050506005.

REFERENCE (1) Cheng, W. L.; Han, F. Y.; Liu, Q. N.; Zhao, R.; Fan, H. L. Experimental and Theoretical Investigation of Surface temperature Non-Uniformity of Spray Cooling. Energy 2011, 36, 249-257. (2) Cheng, W. L.; Zhang, W. W.; Chen, H.; Hu, L. Spray Cooling and Flash Evaporation Cooling: The Current Development and Application. Renew. Sust. Energy Rev. 2016, 55, 614-628. (3) Talbot, E. L.; Yang, L.; Berson, A.; Bain, C. D. Control of the Particle Distribution in

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

Inkjet Printing through an Evaporation-Driven Sol-Gel Transition. ACS appl. mater. inter. 2014, 6, 9572-9583. (4) Singh, M.; Haverinen, H. M.; Dhagat, P.; Jabbour, G. E. Inkjet Printing-Process and Its Applications. Adv. Mater. 2010, 22, 673-685. (5) Lim, T.; Han, S.; Chung, J.; Chung, J. T.; Ko, S.; Grigoropoulos, C. P. Experimental Study on Spreading and Evaporation of Inkjet Printed Pico-Liter Droplet on A Heated Substrate. Int. J. Heat Mass Tran. 2009, 52, 431-441. (6) Zhang, S. J.; Li, Q. W.; Kinloch, I. A.; Windle, A. H. Ordering in A Droplet of an Aqueous Suspension of Single-Wall Carbon Nanotubes on A Solid Substrate. Langmuir 2010, 26, 2107-2112. (7) Marin, A. G.; Gelderblom, H.; Lohse, D.; Snoeijer, J. H. Order-to-Disorder Transition in Ring-Shaped Colloidal Stains. Phys. Rev. Lett. 2011, 107, 085502. (8) Dugas, V.; Broutin, J.; Souteyrand, E. Droplet Evaporation Study Applied to DNA Chip Manufacturing. Langmuir 2005, 21, 9130-9136. (9) Chang, S. T.; Velev, O. D. Evaporation-Induced Particle Microseparations inside Droplets Floating on a Chip. Langmuir 2006, 22, 1459-1468. (10) Picknett, R. G.; Bexon, R. The Evaporation of Sessile or Pendant Drops in Still Air. J. Colloid Interface Sci. 1977, 61, 336-350. (11) McHale, G.; Aqil, S.; Shirtcliffe, N. J.; Newton, M. I.; Erbil, H. Y. Analysis of Droplet Evaporation on a Superhydrophobic Surface. Langmuir 2005, 21, 11053-11060. (12) Xu, W.; Choi, C.H. From Sticky to Slippery Droplets: Dynamics of Contact Line Depinning on Superhydrophobic Surfaces. Phys. Rev. Lett. 2012, 109, 024504.

ACS Paragon Plus Environment

Page 26 of 31

Page 27 of 31 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

(13) Tsai, P.; Lammertink, R. G. H.; Wessling, M.; Lohse, D. Evaporation-Triggered Wetting Transition for Water Droplets upon Hydrophobic Microstructures. Phys. Rev. Lett. 2010, 104, 116102. (14) Sobac, B, Brutin, D. Triple-Line Behavior and Wettability Controlled by Nanocoated Substrates: Influence on Sessile Drop Evaporation. Langmuir 2011, 27, 14999-15007. (15) Kulinich, S.A.; Farzaneh, M. Effect of Contact Angle Hysteresis on Water Droplet Evaporation from Super-hydrophobic Surfaces. Appl. Surf. Sci. 2009, 255, 4056-4060. (16) Nguyen, T. A. H.; Nguyen, A. V.; Hampton, M. A.; Xu, Z. P.; Huang, L. B.; Rudolph, V. Theoretical and Experimental Analysis of Droplet Evaporation on Solid Surfaces. Chem. Eng. Sci. 2012, 69, 522-529. (17) Chen, X. M.; Ma, R. Y.; Li, J. T; Hao, C. L.; Guo, W.; Luk, B. L.; Li, S. C.; Yao, S. H.; Wang, Z. K. Evaporation of Droplets on Superhydrophobic Surfaces: surface Roughness and Small Droplet Size Effects. Phys. Rev. Lett. 2012, 109, 116101. (18) Orejon, D.; Sefiane, K.; Shanahan, M. E. R. Stick-Slip of Evaporating Droplets: Substrate Hydrophobicity and Nanoparticle Concentration. Langmuir 2011, 27, 12834-12843. (19) Wang, F. C.; Wu, H. A. Pinning and Depinning Mechanism of the Contact Line during Evaporation of Nano-Droplets Sessile on Textured Surfaces. Soft Matter 2013, 9, 5703-5709. (20) Zhang, J. G.; Muller-Plathe, F.; Leroy, F. Pinning of the Contact Line during Evaporation on Heterogeneous Surfaces: Slowdown or Temporary Immobilization? Insights from a Nanoscale Study. Langmuir 2015, 31, 7544-7552. (21) Wang, F. C.; Wu, H. A. Molecular origin of contact line stick-slip motion during droplet evaporation. Sci. rep. 2015, 5, 17521.

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

(22) Li, Y. Q.; Wu, H. A; Wang, F. C. Effect of a Single Nanoparticle on the Contact Line Motion. Langmuir 2016, 32,12676-12685. (23) Yu, H.; Chen, J.; Zhu, Y. B.; Wang, F. C.; Wu, H. A. Multiscale transport mechanism of shale gas in micro/nano-pores. Int. J. Heat Mass Tran. 2017, 111,1172-1180. (24) Chen, L.; Kang, Q. J.; Mu, Y. T.; He, Y. L.; Tao, W. Q. A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications. Int. J. Heat Mass Tran. 2014, 76, 210-236. (25) Li, Q.; Zhou, P.; Yan, H. J. Pinning-Depinning Mechanism of the Contact Line during Evaporation on Chemically Patterned Surfaces: A Lattice Boltzmann Study. Langmuir 2016, 32, 9389-9396. (26) Ledesma-Aguilar, R.; Vella, D.; Yeomans, J. M. Lattice Boltzmann Simulations of Droplet Evaporation. Soft Matter 2014, 10, 8267-8275. (27) Wu, J. J. Simulation of rough surfaces with FFT. Tribol. Int. 2000, 33, 47-58. (28) Wu, J. J. Simulation of non-Gaussian surfaces with FFT. Tribol. Int. 2004, 37, 339-346. (29) Manesh, K.; Ramamoorthy, B.; Singaperumal, M. Numerical generation of anisotropic 3D non-Gaussian engineering surfaces with specified 3D surface roughness parameters. Wear 2010, 268, 1371-9. (30) Qian, Y. H.; D'Humières, D.; Lallemand, P. Lattice BGK Models for Navier-Stokes

Equation. Europhys. Lett. 1992, 17, 479-484. (31) Huang, H. B.; Krafczyk, M.; Lu, X. Y. Forcing term in single-phase and Shan-Chen-type multiphase lattice Boltzmann models. Phys. Rev. E 2011, 84, 046710. (32) Li, Q.; Luo, K. H.; Li, X. J. Forcing Scheme in Pseudopotential Lattice Boltzmann

ACS Paragon Plus Environment

Page 28 of 31

Page 29 of 31 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

Model for Multiphase Flows. Phys. Rev. E 2012, 86, 016709. (33) Yuan, P.; Schaefer, L. Equations of State in A Lattice Boltzmann Model. Phys. Fluids 2006, 18, 042101. (34) Li, Q.; Kang, Q. J.; Francois, M. M.; He, Y. L.; Luo, K. H. Lattice Boltzmann Modeling of Boiling Heat Transfer: The Boiling Curve and the Effects of Wettability. Int. J. Heat Mass Tran. 2015, 85, 787-796. (35) Anderson, D. M.; McFadden, G. B.; Wheeler, A. A. Diffuse-interface Methods in Fluid Mechanics. Annu. Rev. Fluid Mech.1998, 30, 139-165. (36) Li, Q.; Zhou, P.; Yan, H. J. Improved Thermal Lattice Boltzmann Model for Simulation of Liquid-Vapor Phase Change. Phys. Rev. E 2017, 96, 063303. (37) He, X. Y.; Zou, Q.; Luo, L. S.; Dembo, M. Analytic Solutions of Simple Flows and Analysis of Nonslip Boundary Conditions for the Lattice Boltzmann BGK model. J. Stat. Phys. 1997, 87,115-136. (38) Guo, Z. L.; Zheng, C. G.; Shi, B. C. Non-equilibrium Extrapolation Method for Velocity and Pressure Boundary Conditions in the Lattice Boltzmann Method. Chinese Phys. 2002, 11, 366-374. (39) Yuan, W. Z.; Zhang, L. Z. Lattice Boltzmann Simulation of Droplets Impacting on Superhydrophobic Surfaces with Randomly Distributed Rough Structures. Langmuir 2017, 33, 820-829. (40) Zhang, L. Z.; Yuan, W. Z. A Lattice Boltzmann Simulation of Coalescence-Induced Droplet Jumping on Superhydrophobic Surfaces with Randomly Distributed Structures. Appl. Surf. Sci. 2018, 436, 172-182.

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

(41) Callen, H. B. Thermodynamics and an Introduction to Thermostatistics. Second ed., Wiley, New York, 1985. (42) Law, C. K. Recent Advances in Droplet Vaporization and Combustion. Prog. Energy Combust. Sci. 1982, 8, 171-201. (43) Li, Q.; Kang, Q. J.; Francois, M. M.; Hu, A. J. Lattice Boltzmann Modeling of Self-propelled Leidenfrost Droplets on Ratchet Surfaces. Soft Matter. 2016, 12, 302-312. (44) Sui, Y.; Ding, H.; Spelt, PDM. Numerical Simulations of Flows with Moving Contact Lines. Annu. Rev. Fluid Mech. 2014, 46, 97-119. (45) Butt, H. J.; Gao, N.; Papadopoulos, P.; Steffen, W.; Kappl, M.; Berger, R. Energy Dissipation of Moving Drops on Superhydrophobic and Superoleophobic Surfaces. Langmuir 2017, 33, 107-116. (46) Marchand, A.; Weijs, J. H.; Snoeijer, J. H.; Andreotti, B. Why is surface tension a force parallel to the interface? Am. J. Phys. 2011, 79, 999. (47) Li, H. X.; Yang, W. L.; Aili, A.; Zhang, T. J. Insights into the Impact of Surface Hydrophobicity on Droplet Coalescence and Jumping Dynamics. Langmuir 2017, 33, 8574-8581. (48) Li, G. Q.; Alhosani, M. H.; Yuan, S. J.; Liu, H. R.; Ghaferi, A. A.; Zhang, T. J. Microscopic droplet formation and energy transport analysis of condensation on scalable superhydrophobic nanostructured copper oxide surfaces. Langmuir 2014, 30, 14498-14511.

ACS Paragon Plus Environment

Page 30 of 31

Page 31 of 31 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

ACS Paragon Plus Environment