Solid–Liquid Interface Thermal Resistance Affects the Evaporation

Article pubs.acs.org/Langmuir

Solid−Liquid Interface Thermal Resistance Affects the Evaporation Rate of Droplets from a Surface: A Study of Perfluorohexane on Chromium Using Molecular Dynamics and Continuum Theory Haoxue Han,*,† Christiane Schlawitschek,‡ Naman Katyal,¶ Peter Stephan,‡ Tatiana Gambaryan-Roisman,‡ Frédéric Leroy,† and Florian Müller-Plathe†,§ †

Theoretische Physikalische Chemie, Eduard-Zintl-Institut für Anorganische und Physikalische Chemie, Technische Universität Darmstadt, Alarich-Weiss-Straße 8, 64287 Darmstadt, Germany ‡ Institute of Technical Thermodynamics, Technische Universität Darmstadt, Alarich-Weiss-Straße 10, 64287 Darmstadt, Germany ¶ Department of Applied Chemistry, Institute of Technology, BHU, Varanasi 221005, India § Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, United States ABSTRACT: We study the role of solid−liquid interface thermal resistance (Kapitza resistance) on the evaporation rate of droplets on a heated surface by using a multiscale combination of molecular dynamics (MD) simulations and analytical continuum theory. We parametrize the nonbonded interaction potential between perfluorohexane (C6F14) and a face-centered-cubic solid surface to reproduce the experimental wetting behavior of C6F14 on black chromium through the solid−liquid work of adhesion (quantity directly related to the wetting angle). The thermal conductances between C6F14 and (100) and (111) solid substrates are evaluated by a nonequilibrium molecular dynamics approach for a liquid pressure lower than 2 MPa. Finally, we examine the influence of the Kapitza resistance on evaporation of droplets in the vicinity of a three-phase contact line with continuum theory, where the thermal resistance of liquid layer is comparable with the Kapitza resistance. We determine the thermodynamic conditions under which the Kapitza resistance plays an important role in correctly predicting the evaporation heat flux.

I. INTRODUCTION The evaporation of a droplet plays an important role in many industrial applications, such as cooling, printing, washing, and coating, hence the study of the evaporation of droplets has attracted considerable attention. Understanding the heat transfer across solid−liquid interfaces is of great importance when studying the evaporation phenomena. A key quantity that thermally characterizes a solid−liquid interface is the interfacial thermal resistance (Kapitza resistance1) RKap = ΔTsl/q where q is the heat flux and ΔTsl is the temperature difference between the solid and liquid. Over the past decade, considerable efforts have been made to directly probe the Kapitza resistance through both experiments2,3 and molecular simulations.4−15 At the macroscopic scale, the Kapitza resistance for liquids in contact with a solid surface at room temperature can be neglected compared to the resistance of the liquid itself. However, when the characteristic length of the system drops to the nanoscale, the role of the Kapitza resistance becomes much more significant. In the case of an evaporating droplet on a heated surface, the evaporation heat flux is significantly enhanced at the threephase contact line where the liquid film thickness is several tens of nanometers.16 The thermal resistance of liquid layer is very low and therefore comparable to the Kapitza resistance. The region of small spatial dimension between the nonevaporating adsorbed film and the macroscopic meniscus is referred to as the “microregion”.17 The low thermal resistance of the liquid © XXXX American Chemical Society

layer in the microregion leads to a very strong local evaporation rate. As a result, the heat transport and phase change in this region exert a strong influence on overall heat transport and the evaporation rate in a macroscopic system, such as a single groove,18 a growing bubble,16 or a sessile droplet.19 To describe the stationary shape of a thin liquid film in the microregion of an evaporating droplet and predict the evaporation heat flux, the Kapitza resistance should be properly taken into account in the numerical model. Here, we study the role of solid−liquid interface thermal resistance (Kapitza resistance) on the evaporation rate of droplets on a heated surface by using a multiscale combination of molecular dynamics (MD) simulations and analytical continuum theory. At the molecular level, we parametrize the nonbonded interaction potential between perfluorohexane (C6F14) and a face-centered-cubic solid surface to reproduce the experimental wetting behavior of C6F14 on black chromium through the solid−liquid work of adhesion (quantity directly related to the wetting angle). We investigated two different atoms alignments of the solid surface: the (100) and (111) direction. The interfacial thermal conductance between C6F14 and both solid substrates is evaluated by a nonequilibrium molecular dynamics approach for a liquid pressure lower than 2 MPa. The obtained Kapitza resistance is entered into the Received: April 25, 2017 Published: May 11, 2017 A

DOI: 10.1021/acs.langmuir.7b01410 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

Fourier series for torsional energetics, and Coulomb and LennardJones terms for the nonbonded interactions. The parameters are the force constants k, the r0 and θ0 reference values, the Fourier coefficients V, the partial atomic charges q, and the Lennard-Jones radii and well-depths, σ and ϵ.

macroscopic continuum evaporation equations to examine the influence of the Kapitza resistance on evaporation of droplets in the vicinity of a three-phase contact line. We determine the thermodynamic conditions under which the Kapitza resistance plays an important role in correctly predicting the evaporation heat flux. To achieve a reliable determination of the Kapitza resistance, a precise estimation of the interfacial interaction is required. For a weakly polar solid surface, the solid−liquid interaction is dominated by van der Waals (vdW) forces, which can be modeled by using the pairwise 12-6 Lennard-Jones (LJ) potential Uij = 4ϵslij [(σ/rij)12 − (σ/rij)6] where ϵslij is the binding energy between a solid atom i and a liquid atom j with the intermolecular distance rij, and σ is the molecular diameter. The solid−liquid potential well depth ϵslij is often calculated with the Lorentz−Berthelot mixing rule ϵijsl = λ ϵlϵs where ϵl (ϵs) is the liquid−liquid (solid−solid) potential depth. λ is a tuning coefficient that controls the liquid−solid affinity. In previous studies, λ was varied to calculate the Kapitza lengths of model LJ liquids on a solid surface with different wettability.5,22 The quality of λ determines the interface physics, including the reliability of the obtained Kapitza resistance. A commonly adopted method to characterize the solid−liquid interactions ϵsl in molecular simulations is to reproduce the experimental contact angle θ of a liquid droplet on the substrate.6,20,23 Although widely employed, such a method represents strong limitations. Contact angles cannot be used to characterize wetting surfaces where θ is too small to be directly measured. Conducting such a parametrization with molecular dynamics simulation would require a large solid surface with an extended droplet and hence can be very demanding with respect to computational resources. An alternative equilibrium thermodynamic quantity which describes a solid−liquid interface is the work of adhesion Wsl, which is defined as the reversible work per unit area to separate the liquid and solid until the interaction vanishes. Under a low vapor pressure, the work of adhesion can be expressed in terms of the liquid−vapor surface tension and the contact angle as Wsl = γlv(1 + cos θ).24 Wsl was suggested as a target quantity to optimize the nonbonded interfacial interaction for solid−liquid interfaces25 and was successfully used to parametrize watergraphene24 and n-hexane-graphene31 interfaces. Since Wsl is closely related to the interfacial coupling, it has a direct influence on the solid−liquid thermal transport. Shenogina et al. found a linear relation between Kapitza conductance G at self-assembled monolayer (SAM)-water interfaces and the work of adhesion in molecular dynamics simulations,4 while Park and Cahill recently observed a trend of saturation of G for large Wsl by measuring the interfaces between a water−ethanol mixture and gold nanodisks coated with SAM.3 A similar saturation for solid−water interfaces was confirmed by Merabia et al. through analytical calculations based on a generalized acoustic mismatch model.21

E bond =

∑ kb , i(ri − r0, i)2 i

Eangle =

∑ kθ ,i(θi − θ0,i)2 i

Etorsion =

∑ [1/2V1, i(1 + cos ϕi) + 1/2V2, i(1 − cos 2ϕi) i

+ 1/2V3, i(1 + cos 3ϕi) + 1/2V4, i(1 − cos 4ϕi)]

⎡ q q e2

Enonbond =

∑ ∑ ⎢⎢ i

j>i

i j

⎢⎣ rij

⎡⎛ ⎞12 ⎛ ⎞6 ⎤⎤ σij σij ⎥ + 4ϵij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥⎥ r ⎝ rij ⎠ ⎦⎥⎦ ⎣⎝ ij ⎠

Each atom has an associated two-letter atom type that is used to designate the parameters for atom pairs (bond stretching) or atom triplets (angle bending). The atom types used here are CT (aliphatic sp3 carbon) and F (fluorine). The OPLS-AA parameters for perfluoroalkanes are reported in Tables 1−4.

Table 1. OPLS-AA Bond Stretching Parameters for Perfluorohexane bond

kb (kcal mol−1 Å−2)

r0 (Å)

CT−F CT−CT

367.0 268.0

1.332 1.529

Table 2. OPLS-AA Angle Bending Parameters for Perfluorohexane angle

kθ (kcal mol−1 rad−2)

θ0 (deg)

F−CT−F CT−CT−F CT−CT−CT

77.0 50.0 58.35

109.10 109.50 112.70

Table 3. OPLS-AA Nonbonded Parameters for Perfluorohexane atom type

atom or group

q (e−)

σ (Å)

ϵ (kcal mol−1)

F CT CT CT

F CF3 CF2 CF

−0.12 0.36 0.24 0.12

2.95 3.50 3.50 3.50

0.053 0.066 0.066 0.066

Table 4. OPLS-AA Fourier Coefficients (kcal mol−1) for Perfluorohexane

II. MOLECULAR DYNAMICS MODEL AND COMPUTATIONAL DETAILS

dihedral angle

V1

V2

V3

V4

F−CT−CT−F CT−CT−CT−F CT−CT−CT−CT

−2.500 0.300 6.622

0.000 0.000 0.948

0.250 0.400 −1.388

0.000 0.000 −2.118

The FCC (100) surface is modeled by using the 12-6 LJ potential with σ = 3.405 Å and ϵ = 5.29 kcal/mol.27 The (111) surface is modeled by using a gold lattice with the ReaxFF potential.28 We note that the interaction ϵsl between the solid and liquid is of most interest for determining the interfacial thermodynamic and thermal transport properties, rather than the interaction within the solid ϵs. Nosé− Hoover-like thermostat and barostat were used with a coupling time constant of 0.2 and 2 ps for temperature and pressure, respectively. A cutoff distance of 1.2 nm for the LJ potential and a time step of 2 fs

The all-atom optimized potential for liquid simulations (OPLS-AA) of n-perfluorohexane was employed to model the interaction of the C6F14 liquid.26 n-Perfluorohexane is the primary component of an important cooling fluid FC-72 that is often used in fundamental evaporation and boiling studies.16 The molecular dynamics simulations were conducted by using the LAMMPS package. The potential energy function consists of harmonic bond stretching and angle bending terms, a B

DOI: 10.1021/acs.langmuir.7b01410 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir were used in all simulations. The Coulomb interaction was truncated with a cutoff distance of 1.4 nm. Three simulation setups were constructed for different situations including a liquid−vapor, a solid−liquid−vapor, and a solid−liquid− solid system. The liquid−vapor and solid−liquid−vapor systems were used to determine the surface tension of C6F14 in section III. Then the solid−liquid−vapor system was used to optimize the LJ energy parameter between C6F14 and solid in section IV. The solid−liquid− solid system was used to calculate the interfacial thermal conductance in section V. The model system size, simulation time, and simulation procedures are described in the respective sections.

III. LIQUID−VAPOR SURFACE TENSION OF PERFLUOROHEXANE We first validate the potential for the C6F14 liquid by calculating the liquid−vapor surface tension γlv. The surface tension is obtained from the spatial integration of the stress-tensor anisotropy through γlv =

1 2

Figure 1. (a) Liquid−vapor surface tension γlv of C6F14 determined from a free-standing liquid film and a confined liquid film at 290 K. The results from five independent atomic trajectories are presented. λ = 10−6 was chosen for the solid−liquid−vapor system with only repulsive solid−liquid interaction. The error bars correspond to the standard deviations of 5 calculations, each with a duration of 4 ns. The experimental values are taken from ref 30. Stress-tensor distribution across (b) confined liquid film systems and (c) the free-standing liquid film.

L

∫0 x ⟨P⊥(x) − P (x)⟩dx where x is

the direction normal to the surface, Lx is the dimension of the liquid system in the x direction, and P⊥(x) and P∥(x) are the stress tensor components normal and parallel to the liquid interface, respectively.29,31 Two systems are constructed for the calculation of γlv. The first one is a free-standing liquid film that contains 1000 C6F14 molecules (20,000 atoms) in a 25 nm × 5 nm × 5 nm box which is elongated in the x direction. The system was first at equilibrium at 290 K and 1 atm in the fixed atom number N, system pressure P, and temperature T ensemble (NPT) for 1 ns to reach liquid−vapor equilibrium. Then the system is switched to the fixed atom number, system volume, and temperature (NVT) ensemble to sample the liquid pressure tensor for 20 ns. The second system confines the liquid film between two repulsive solid walls. It contains 500 C6F14 molecules and a solid substrate of 3200 atoms arranged in 16 monatomic layers. The molecules are in contact with the (100) surface of the solid. The surface tension was calculated when the energy parameter ϵsl was turned to a very small value such that the solid−liquid interaction is purely repulsive. The results of the surface tension calculated from these two models and the experimental data are shown in Figure 1. The surface tension of ≈9.5 mN/m from our all-atomic molecular dynamics simulations underestimates the experimental value of ≈12 mN/m by ca. 20%, in agreement with a previous study.31 Therefore, the OPLS-AA potential employed in this work yields a relatively satisfactory agreement with experiments.30 The surface tension for liquid−vapor surface tension in the confinement was calculated and compared with the free liquid case to verify that in the limit of weak solid−liquid coupling, the solid−liquid interfacial tension and the liquid surface tension are identical. This equality ensures the validity of the drysurface approach for calculating the solid−liquid work of adhesion. The stress-tensor distributions are calculated for the same systems as used for surface tension calculations. Both systems are divided into 301 bins in the axial direction (perpendicular to the liquid−vapor interface). The stress tensor for each liquid atom in the system is calculated, by following the algorithm documented in the manual of LAMMPS. The per-atom stress tensor is then summed for all the atoms in each bin. The stress-tensor distributions for the free-standing liquid film and confined liquid film systems are hence obtained by averaging the stress tensor of each bin for 10 ns.

IV. PARAMETRIZATION OF INTERFACIAL ENERGY PARAMETER The solid−liquid−vapor system contains 500 C6F14 molecules and a solid substrate of 3200 atoms arranged in 16 monatomic layers with a 10 nm-long vapor phase of the same cross section. We calculate the solid−liquid work of adhesion Wsl as a function of the energy parameter ϵsl between C6F14 and the solid by using a dry-surface approach. Then we were able to locate the value of ϵsl corresponding to the experimental work of adhesion Wexp sl . The n-C6F14 system was brought into contact with the solid surface. The solid−liquid−vapor system was used to calculate Wsl, as shown in the inset of Figure 2. Wsl is calculated by integrating the free energy change per unit area to turn the solid−liquid interaction into a repulsive potential Wsl = −

∫1

2

⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σij σij ∑ ∑ 4 ϵi ⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ r ⎝ rij ⎠ ⎦ i=1 j=1 ⎣⎝ ij ⎠ Nl

Ns

dλ ϵs S

where S is the contact area of the solid−liquid interface, and Nl(s) refers to the number of liquid (solid) atoms. “1” and “2” represent the effective repulsive state and the actual interaction state, respectively. Precisely, “1” represents the purely repulsive state between the solid and liquid with the parameter λ ϵs = 1.0 × 10−6 (kcal/mol)1/2. “2” represents the actual interaction state where 1.0 × 10−6 ≪ λ ϵs < 1 (kcal/mol)1/2. λ ϵlϵs is the interaction energy between the solid and liquid, and ϵl depends on the liquid atom type, i.e. ϵl can be ϵF for fluor atom or ϵC for carbon atom. We denote the quantity between the angle brackets as I. In Figure 2, we plot ⟨I⟩ and the associated work of adhesion versus λ ϵs . The experimental work of adhesion Wexp sl can be determined from the relation exp exp Wexp sl = γlv (1+cos θ) where γlv and θ are the experimental liquid−vapor surface tension of C6F14 and the static contact angle of a C6F14 droplet on the black chromium surface.19 In the wetting experiment under isothermal conditions, θ was −1 30 measured to be ≈3° and γexp lv = 11.91 mN m . Therefore, the experimental work of adhesion was determined as Wsl = 23.8 C

DOI: 10.1021/acs.langmuir.7b01410 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

Figure 3. (a) Cumulative input and output energy Ecum as a function of time t in the hot and cold heat reservoirs. (b) Histogram of the center-of-mass velocities v ̅ of C6F14 molecules at the center of the system. The black solid curve refers to a normalized Gaussian function

(

mv 2

)

f ( v ̅ ) ∝ exp − 2k̅ T̅ . (c) Temperature distribution T(x) along the B

longitudinal direction ix⃗ of the confined system. The error bars correspond to the standard deviations of temperatures of each bin. The red line is a linear fit to the temperature points in the C6F14 liquid. The temperature of the hot and cold reservoirs is kept constant at 310 and 270 K, respectively. The temperature drop ΔTsl at the interface is estimated by extrapolating the linear fits to the temperature profiles in the liquid and calculating the difference at the interface.

Figure 2. Variation of (a) the time-averaged integrand ⟨I⟩ and (b) the solid−liquid work of adhesion versus λ ϵs at 290 K for the (100) and (111) solid surfaces. The dashed line is a guide to the eye. For a work of adhesion of Wsl = 23.8 mJ m−2, the solid−liquid interaction parameter is determined as λ ϵs = 0.650 (kJ mol−1)1/2 for the (100) surface and λ ϵs = 0.725 (kJ mol−1)1/2 for the (111) surface.

Hoover thermostats. For the system to reach the desired pressure, the equations of motion of the atoms are integrated in the NPT ensemble at constant pressure and temperature for enough long simulation time (a few tens of nanoseconds). Then the system is switched to the NVE ensemble for 20 ns. The total energy of the system is monitored to be conserved to conduct the steady-state NEMD simulation in the following stage: the solid atoms in the atomic planes at the two ends of the system are frozen, and two thermostats with different temperatures are imposed on the two solid surfaces to establish the steady-state heat flux in the normal direction to the solid− liquid interface through the system. The heat source and sink in the solid leads were introduced to induce a steady-state heat flux. The system was divided into slab bins in the longitudinal direction with a bin size of 10 Å. The temperature of each bin was calculated and averaged over 20 ns. The cumulative input and output energy Ecum(t) in the hot and cold heat reservoirs show a good system energy conservation in Figure 3(a). It is very important to confirm the local thermal equilibrium where the temperature is defined in this particle system. We calculated the atomic velocity distribution in the central part of the system at steady state. The entire system is divided into onedimensional bins in the normal direction to the solid−liquid interface with a bin width of 10 Å. The center-of-mass (COM) velocities v ̅ (components in the three axes of Cartesian coordinates) of C6F14 molecules inside the bin at the center of the system are recorded and averaged for 500 ps. The histogram of the velocities is shown in Figure 3(b). It can be observed that the COM velocity follows the Maxwell− Boltzmann distribution with the probability density of

mJ m−2. From Figure 2(b), we therefore determine the corresponding solid−liquid interaction parameter as λ ϵs = 0.650 (kJ mol−1)1/2 for the (100) surface and λ ϵs = 0.725 (kJ mol−1)1/2 for the (111) surface. We shall employ these parameters to investigate the interfacial thermal conductance between the C6F14 liquid and solid.

V. INTERFACIAL THERMAL CONDUCTANCE We used a nonequilibrium setup to study heat flow across the solid-C6F14 interface. Our system contains two solid-C6F14 interfaces, created by coupling the C6F14 liquid to the (100) and (111) surfaces, as is shown in Figure 3. For the (100) surface, each solid part contains 2160 solid atoms arranged in an 8.1 nm long cuboid in the x direction with a cross section of 3.24 × 3.24 nm2. For the (111) surface, each part contains 1800 solid atoms with a cross section of 3.04 × 2.98 nm2. Both systems include 500 C6F14 molecules. Periodic boundary conditions are imposed in the direction parallel to the interface plane (transverse to the current flow). Ten monolayers of solid atoms with a length of Lbath = 5.4 nm next to the frozen atoms at the left and right ends are coupled to Langevin heat baths at temperatures Th = T + ΔT/2 and Tc = T − ΔT/2, where T = 290 K and ΔT = 40 K. The bath time constant is chosen as tb = 1 ps, ensuring that the bath-induced mean free path of phonons satisfies Λb = vgtb ≤ 2Lb (vg = 3260 m/s is the speed of sound) so that phonons arriving at thermostats are dissipated before reflecting from the fixed ends. To prepare the systems for thermal transport studies, we first performed NPT runs for all systems at 1 atm and 300 K for 1 ns. Anisotropic pressure coupling allowed independent variation of the dimension perpendicular to the surface. The NPT run was followed by a NVT run for 1 ns, where the solid and C6F14 temperatures were maintained at 300 K using Nosé−

(

mv 2

)

f ( v ̅ ) ∝ exp − 2k̅ T̅ , where m̅ is the mass of a C6F14 molecule, B

kB is the Boltzmann constant, and T is the temperature. The temperature T corresponding to the velocity distribution is D

DOI: 10.1021/acs.langmuir.7b01410 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

in cooling applications,16 the thermal conductance is rather insensitive to the liquid pressure. Despite the interfacial interaction being parametrized with the same work of adhesion, we observe that the (111) solid surface has higher thermal conductances than the (100) surface for all pressures beyond the range of error bars. This could be explained first by a higher surface atomic density on the (111) surface n(111)/n(100) = 2/ √3. In addition, the mass density distributions as plotted in Figure 4(b) and (c) suggest that the first adsorbed layers of nC6F14 have a higher population at the (111) surface compared to the (100) surface. This is in line with the observation of water density on silicon and gold surfaces.7 Therefore, we show that apart from solid−liquid affinity,4 the atomic topology in the vicinity of the interface is also an important factor to define interfacial thermal transport. At a higher temperature of T = 320 K slightly below the saturation point of n-C6F14 at Tsat = 329 K, G remains close to the values at room temperature as shown in Figure 4(a). The three most important factors affecting the thermal conductance are (1) the wall−fluid interaction strength and (2) the influences of surface atomic density,14 as well as (3) the density of first adsorbed liquid layer.15,35 In the pressure range under investigation in this work, we do not observe a considerable increase of the density of the first adsorbed liquid layer, which explains the insensitivity of the thermal conductance with the pressure. In the literature, it was shown that higher liquid pressures can effectively increase the density of the first adsorbed liquid layer, thus changing the thermal conductance.35 Nevertheless, the pressure regime where such considerable difference is observed is around several hundreds of megapascals (for example water in ref 35), which is far beyond the pressure range of droplet evaporation under normal conditions.

determined to be 289 K which agrees relatively well with the local temperature of the central bin. The linear profile demonstrates the establishment of steadystate heat flow in the system. The heat flux through the system is calculated as q = 1/A·∂Ecum/∂t. A typical temperature profile obtained from the nonequilibrium simulation is shown in Figure 3(c). The temperature drop ΔTsl at the interface is estimated by extrapolating the linear fits to the temperature profiles in the liquid and calculating the difference at the interface. The thermal conductance can be obtained as G = q/ ΔTsl. We note that the measured temperature jump is very sensitive with the definition of the location of the solid−liquid interface. A solid−liquid boundary defined at the innermost solid layer in contact with liquid would result in a smaller temperature drop compared to a boundary defined at the first peak of the liquid density, and hence a lower Kapitza resistance would be obtained. In the present work, the temperature gradient in the liquid corresponding to the highest thermal conductance is measured to be ∇T = 0.174 K/Å. The average distance between the innermost solid layer and the first liquid density peak is around 4.7 Å, hence resulting in a temperature difference of 0.818 K. Therefore, the largest error introduced by the definition of the interface position is ±0.409 K. Since the temperature profile in the solid can be considered to be constant, the error of the interface temperature jumps due to the definition of the interface position ϵT = ±0.409 K. The error introduced is 7.8% for the highest thermal conductance and 4.2% for the lowest one, within the error bars reported in Figure 1. We hence conclude that the position of the solid− liquid interface is not crucial for the thermal conductance. Therefore, in our work, we define the solid−liquid boundary at the midpoint of the innermost solid layer in contact with liquid and first peak of the liquid density. Different thermostat temperatures ranging from 30 to 50 K were applied to the systems, and a linear dependence was found between heat flux q and ΔTsl, indicating that the system under investigation was in the linear response regime. The thermal conductances G at the solid-n-C6F14 interface with respect to liquid pressures p below 2 MPa are shown in Figure 4(a) for the (100) and (111) solid surfaces. In this pressure regime, which corresponds to the operation condition

VI. EVAPORATION IN THE VICINITY OF THE THREE-PHASE CONTACT LINE The values of solid−liquid interfacial resistance (inverse of thermal conductance) computed in previous sections are small at a macroscopic scale. For simplicity, let us assume G = 10 MWm−2 K−1 (or the Kapitza resistance RKap = G−1 = 1 × 10−7 Km2 W−1). The equivalent thermal resistance of a liquid layer of C6F14 results in a thickness of around 5 nm (assuming a thermal conductivity value κ = 5 × 10−2 Wm−1 K−1). Therefore, we suggest that the Kapitza resistance is important in processes for film thickness lower than 100 nm. One of the typical processes for which the Kapitza resistance may play an important role is evaporation of an ultrathin liquid film in the vicinity of a three-phase contact line. Figure 5 (upper panel) illustrates schematically the structure of the liquid−gas interface in the vicinity of an apparent three-phase contact line, which separates a liquid drop (Figure 5, lower panel) from an apparently dry part of a solid substrate. The gas phase is a pure vapor at saturation, and the wall is hotter than the saturated vapor. If the liquid is highly wetting (as for example, C6F14 on the most technical substrates, including chromium surface), the apparently dry substrate area is covered by a very thin nonevaporating adsorbed film of constant thickness, which is shown at the left part of the image. The adsorbed film, which has typically a thickness of several nanometers, cannot evaporate due to strong adhesion forces between the substrate and the liquid molecules.17,18 The local thickness of the film increases continuously and smoothly in the direction of the droplet center. The vicinity of the apparent contact line on a microscopic scale is a region of a liquid film, in which the angle

Figure 4. (a) Interfacial thermal conductance G of n-C6F14 versus the liquid pressure p against the (100) and (111) solid surfaces. The error bars correspond to the standard deviations of five independent trajectories. The mass density distribution ρ(x) of n-C6F14 at p = 0.39 MPa and T = 300 K in the direction normal to the interface on the (b) (100) and (c) (111) surfaces. E

DOI: 10.1021/acs.langmuir.7b01410 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

Finally, the liquid−vapor interfacial resistance Rint can be computed from the relation33 R int =

2 − f Tsat 2πRTsat 2f hlv2 ρv

(3)

where ρv is the density of the saturated vapor, R is the gas constant, and f is the evaporation coefficient. The liquid−vapor interfacial resistance is a strong function of the pressure p, since Tsat and ρv increase with increasing pressure, whereas hlv decreases with increasing pressure and tends to zero as the pressure approaches the critical pressure pcr. In Figure 6(a) both the interfacial resistances are plotted for

Figure 5. Scheme of the microregion (upper panel). The microregion is in the vicinity of the apparent three-phase contact line (lower panel). h and ξ denote the liquid film thickness and the coordinate along the heated wall toward the center of the droplet.

of inclination of the liquid−vapor interface changes from zero (in the adsorbed film region) to the apparent contact angle formed by the droplet with the substrate (denoted by θ in Figure 5). This region, typically having a length of an order of micrometer, is referred to as microregion.18 In order to compute the shape of the film in the microregion following the continuous models presented in the literature,17,18 consider a film of thickness h, which can vary along a direction ξ, covering a wall of temperature Tw, which exceeds the saturation temperature Tsat. The film thickness varies in direction ξ, which is parallel to the substrate and points toward the drop center. Assume that the film evaporates into a pure vapor. If the film is flat, so that the absolute value of derivative h′ = dh/dξ is small, the evaporative heat flux qev can be computed from the following equation:16

(

Tw − Tsat 1 + qev =

Δp ρl hlv

Figure 6. (a) Kapitza resistance RKap from MD simulations and liquid−vapor interfacial resistance Rint computed from eq 3 for C6F14 using different values of f. (b) Distribution of the film thickness in the microregion. Empty symbols refer to simulations without RKap, while filled symbols denote simulations including RKap. Computations performed for C6F14 with f = 1 and ΔTw = Tw − Tsat = 3 K. (c) Relative values of the liquid thermal resistance, liquid−vapor interfacial resistance, and Kapitza resistance at different positions of the microregion ( f = 1, ΔTw = 3 K, p/pcr = 0.5). (d) Integrated heat flux per unit length of the contact line for f = 1 and different values of reduced pressure p/pcr. Empty symbols refer to simulations without RKap, while filled symbols denote simulations including RKap.

)

R liq + R int + RKap

(1)

Equation 1 is an extension of an expression for evaporative heat flux from the literature,18 with the solid−liquid interfacial resistance RKap added to the thermal resistance of the liquid layer Rliq = h/κ and the liquid−vapor interfacial resistance Rint. In the above ρl and hlv denote the density of liquid and latent heat of evaporation, respectively, and Rliq = h/κ is the thermal resistance of the liquid layer. The symbol Δp denotes the sum of a capillary pressure difference and disjoining pressure32 h″ A Δp = γlv + 3 2 1.5 (1 + h′ ) h

C6F14 as functions of pressure, whereas the Kapitza resistance has been assumed to stay constant in the relevant range. In these computations we have assumed the value RKap = 10−7 Km2 W−1, which corresponds to the range of G displayed in Figure 4(a). The parameter f may vary between zero and unity, and its value is unknown for most systems. Therefore, in this work we present the results for three different values of f. It is clearly seen that the liquid−vapor interfacial resistance exceeds the Kapitza resistance in the whole range of pressures. It is also clear that Rint decreases with increasing f. The dependence of Rint on pressure is nonmonotonous. The values of Rint and RKap are close to each other for f = 1 in the pressure range 0.3 < p/pcr < 0.7. It can be concluded that in the case of C6F14 the film evaporation rate can never be governed by the Kapitza

(2)

where A denotes the Hamaker constant. For the combination of C6F14 and chromium,34 the Hamaker constant is A = 4.37 × 10−21 J. The disjoining pressure expresses the action of the adhesion forces between the substrate and the liquid molecules which prevent the evaporation of the adsorbed film. The action of the disjoining pressure reduces the volatility of the liquid molecules and thus reduces the evaporation rate. Analysis of eqs 1 and 2 shows that the evaporation rate and evaporation heat flux decrease with decreasing of film thickness and is equal to zero as the film thickness reaches a critical value (it is referred to as adsorbed film thickness). F

DOI: 10.1021/acs.langmuir.7b01410 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

liquid layer becomes even lower than the resistances at both interfaces. However, the evaporation rate in this part of the microregion is very low due to the action of disjoining pressure. One of the most important characteristic values relevant to evaporation in the vicinity of the three-phase contact line is the integrated heat flux per unit length of the contact line, which can be computed from the relation18,19

resistance. However, the Kapitza resistance can affect the evaporation rate, as soon as its value is comparable with the liquid−vapor interfacial resistance (this condition is illustrated in Figure 6(a)) and with the thermal resistance of the liquid film (this condition is fulfilled for films thinner than 10 nm). The distribution of film thickness and, therefore, of the evaporation rate in the microregion (h ≤ 100 nm) between the adsorbed film and the macroscopic meniscus is determined by the balance between the evaporation of liquid in the microregion and the flow of liquid from the macroscopic meniscus into the microregion. In the framework of lubrication approximation the distribution of the film thickness in steady state can be determined from the following equation18 qev hlv

+

1 d ⎛ 3 dΔp ⎞ ⎜h ⎟=0 3νl dξ ⎝ dξ ⎠

Q̇ =

∫0

ξmr

qev dξ

(5)

where ξmr is the length of the microregion. The integrated heat flux as a function of wall overheat ΔTw = Tw − Tsat is represented in Figure 6(d). It is clearly seen that addition of the Kapitza resistance leads to the decrease of the integrated heat flux. The absolute contribution of the Kapitza resistance increases with the increasing wall overheat. However, the influence of the choice of the condensation coefficient f is much stronger than that of the addition of the Kapitza resistance.

(4)

where ξ refers to the coordinate along the wall with ξ = 0 as the contact line, νl denotes the kinematic viscosity of the liquid phase, and qev and Δp are defined in eqs 1 and 2, respectively. The distribution of film thickness in the microregion obtained by a numerical solution of eq 4 for C6F14 with and without the Kapitza resistance (assuming RKap = 0 in eq 1) is depicted in Figure 6(b). It is seen that the Kapitza resistance leads to decreasing of the film thickness at all locations and to decreasing of the apparent contact angle which corresponds to the angle of inclination of the liquid−vapor interface at the edge of the microregion. The decreasing of film thickness and inclination angle can be explained as follows. Additional thermal resistance (the Kapitza resistance) leads to decreasing of the evaporation rate for each value of film thickness. As a result, the flow rate from the macroscopic meniscus to the microregion decreases. This flow is driven by the gradients of capillary pressure and disjoining pressure, both of them determined by the derivatives of film thickness. Decreasing the flow rate means decreasing of the film thickness gradient or the inclination angle of the liquid−vapor interface. In Figure 6(b) the film thickness in the microregion is shown. The adsorbed film exists for ξ < 0, while the macroscopic drop is adjacent to the microregion (ξ > 0.5 μm). For all values of the reduced pressure, p/pcr, the values of the film thickness at all locations within the microregion are reduced as soon as the Kapitza resistance is taken into account. This means that the Kapitza resistance leads to decreasing of the film thickness at all locations. Additionally, it is clearly seen that the inclination angle of the h(ξ) curve(liquid−vapor interface) at the right edge of the microregion (ξ = 0.5), which corresponds to the apparent contact angle, decreases for the simulation cases where the Kapitza resistance is taken into account. This means that the Kapitza resistance leads to decreasing of the apparent contact angle. It is instructive to examine the relative values of the liquid thermal resistance, liquid−vapor interfacial resistance, and Kapitza resistance at different positions of the microregion (see Figure 6(c)). The relation between RKap and Rint stays constant in the whole region. The liquid−vapor interfacial resistance is around 25% higher than the Kapitza resistance. At the right edge of the microregion, where it merges to the macroscopic meniscus, the values of RKap and Rint are much smaller than Rliq. The relative values of RKap and Rint increase monotonously as the film thickness decreases in the direction of adsorbed layer. As the film thickness approaches the minimal value (the thickness of adsorbed layer), the relative resistance of

VII. CONCLUSIONS We investigated the role of the Kapitza resistance on the evaporation rate of droplets on a heated surface by using a multiscale combination of molecular dynamics (MD) simulations and analytical continuum theory. The thermal conductances between C6F14 and (100) and (111) solid substrates are evaluated by a nonequilibrium molecular dynamics approach for a liquid pressure lower than 2 MPa. The obtained Kapitza resistance is entered into the macroscopic continuum evaporation equations to examine the influence of the Kapitza resistance on evaporation of droplets in the vicinity of a three-phase contact line. We determine the thermodynamic conditions under which the Kapitza resistance plays an important role in correctly predicting the evaporation heat flux.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Haoxue Han: 0000-0002-2637-118X Frédéric Leroy: 0000-0002-8470-5811 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Support from the Sonderforschungsbereich Transregio 75 of the Deutsche Forschungsgemeinschaft is gratefully acknowledged. F.M.P. thanks Athanassios Panagiotopoulos and his group for the great hospitality at Princeton University. We thank Vikram Ardham for his help on the technical details of the dry-surface calculation. We thank the Lichtenberg-High Performance Computer at TU Darmstadt for the computing ressources.

■

REFERENCES

(1) Kapitza, P. L. The study of heat transfer in helium II. J. Phys. (USSR), 1941, 4 (1−6), 181−210. (2) Ge, Z.; Cahill, D. G.; Braun, P. V. Thermal conductance of hydrophilic and hydrophobic interfaces. Phys. Rev. Lett. 2006, 96 (18), 186101. (3) Park, J.; Cahill, D. G. Plasmonic Sensing of Heat Transport at SolidLiquid Interfaces. J. Phys. Chem. C 2016, 120, 2814.

G

DOI: 10.1021/acs.langmuir.7b01410 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

Example of Water on Graphene Surfaces. J. Phys. Chem. C 2015, 119, 28470. (25) Leroy, F.; Müller-Plathe, F. Dry-Surface Simulation Method for the Determination of the Work of Adhesion of SolidLiquid Interfaces. Langmuir 2015, 31, 8335. (26) Watkins, E. K.; Jorgensen, W. L. Perfluoroalkanes: Conformational analysis and liquid-state properties from ab initio and Monte Carlo calculations. J. Phys. Chem. A 2001, 105 (16), 4118−4125. (27) Heinz, H.; Vaia, R. A.; Farmer, B. L.; Naik, R. R. Accurate simulation of surfaces and interfaces of face-centered cubic metals using 12−6 and 9−6 Lennard-Jones potentials. J. Phys. Chem. C 2008, 112 (44), 17281−17290. (28) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. A 2008, 112 (5), 1040−1053. (29) Gloor, G. J.; Jackson, G.; Blas, F. J.; de Miguel, E. Test-area simulation method for the direct determination of the interfacial tension of systems with continuous or discontinuous potentials. J. Chem. Phys. 2005, 123 (13), 134703. (30) Dataphysics. http://www.surface-tension.de/ (accessed May 12, 2017). (31) Ardham, V. R.; Deichmann, G.; van der Vegt, N. F.; Leroy, F. Solid-liquid work of adhesion of coarse-grained models of n-hexane on graphene layers derived from the conditional reversible work method. J. Chem. Phys. 2015, 143 (24), 243135. (32) Deryaguin, B. V.; Nerpin, S. V.; Churayev, N. V. Effect of film transfer upon evaporation of liquids from capillaries. Bull. Rilem 1965, 29, 93−97. (33) Schrage, R. W. A theoretical study of interphase mass transfer; Columbia University Press: New York, 1953. (34) Batzdorf, S. Heat transfer and evaporation duringsingle drop impingement onto a superheated wall. Ph.D. Thesis, TU Darmstadt, 2015. (35) Alexeev, D.; Chen, J.; Walther, J. H.; Giapis, K. P.; Angelikopoulos, P.; Koumout-sakos, P. Kapitza resistance between few-layer graphene and water: Liquid layering effects. Nano Lett. 2015, 15 (9), 5744−5749.

(4) Shenogina, N.; Godawat, R.; Keblinski, P.; Garde, S. How wetting and adhesion affect thermal conductance of a range of hydrophobic to hydrophilic aqueous interfaces. Phys. Rev. Lett. 2009, 102, 156101. (5) Barrat, J. L.; Chiaruttini, F. Kapitza resistance at the liquid-solid interface. Mol. Phys. 2003, 101, 1605. (6) Merabia, S.; Shenogin, S.; Joly, L.; Keblinski, P.; Barrat, J. L. Heat transfer from nanoparticles: A corresponding state analysis. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 15113. (7) Pham, A.; Barisik, M.; Kim, B. Pressure dependence of Kapitza resistance at gold/water and silicon/water interfaces. J. Chem. Phys. 2013, 139, 244702. (8) Vo, T. Q.; Kim, B. Transport Phenomena of Water in Molecular Fluidic Channels. Sci. Rep. 2016, 6, 33881. (9) Ho, T. A.; Papavassiliou, D. V.; Lee, L. L.; Striolo, A. Liquid water can slip on a hydrophilic surface. Proc. Natl. Acad. Sci. U. S. A. 2011, 108 (39), 16170−16175. (10) Bin Saleman, A. R.; Chilukoti, H. K.; Kikugawa, G.; Shibahara, M.; Ohara, T. A molecular dynamics study on the thermal transport properties and the structure of the solid liquid interfaces between face centered cubic (FCC) crystal planes of gold in contact with linear alkane liquids. Int. J. Heat Mass Transfer 2017, 105, 168−179. (11) Chilukoti, H. K.; Kikugawa, G.; Shibahara, M.; Ohara, T. Local thermal transport of liquid alkanes in the vicinity of α-quartz solid surfaces and thermal resistance over the interfaces: A molecular dynamics study. Phys. Rev. E 2015, 91 (5), 052404. (12) Giri, A.; Braun, J. L.; Hopkins, P. E. Implications of Interfacial Bond Strength on the Spectral Contributions to Thermal Boundary Conductance across Solid, Liquid, and Gas Interfaces: A Molecular Dynamics Study. J. Phys. Chem. C 2016, 120 (43), 24847−24856. (13) Han, H.; Merabia, S.; Mueller-Plathe, F. Thermal Transport at Solid-Liquid Interfaces: High Pressure Facilitates Heat Flow Through Non-Local Liquid Structuring. J. Phys. Chem. Lett. 2017, 8, 1946− 1951. (14) Ohara, T.; Torii, D. Molecular dynamics study of thermal phenomena in an ultrathin liquid film sheared between solid surfaces: The influence of the crystal plane on energy and momentum transfer at solid-liquid interfaces. J. Chem. Phys. 2005, 122 (21), 214717. (15) Vo, T. Q.; Kim, B. Interface thermal resistance between liquid water and various metallic surfaces. International Journal of Precision Engineering and Manufacturing 2015, 16 (7), 1341−1346. (16) Kunkelmann, C.; Ibrahem, K.; Schweizer, N.; Herbert, S.; Stephan, P.; Gambaryan-Roisman, T. The effect of three-phase contact line speed on local evaporative heat transfer: Experimental and numerical investigations. Int. J. Heat Mass Transfer 2012, 55, 1896. (17) Potash, M.; Wayner, P. C. Evaporation from atwo-dimensional extended meniscus. Int. J. Heat Mass Transfer 1972, 15, 1851−1863. (18) Stephan, P.; Busse, C. A. Analysis of the heat transfer coefficient of grooved heat pipe evaporator walls. Int. J. Heat Mass Transfer 1992, 35, 383−391. (19) Herbert, S.; Fischer, S.; Gambaryan-Roisman, T.; Stephan, P. Local heat transfer and phase change phenomena during single drop impingement on a hot surface. Int. J. Heat Mass Transfer 2013, 61, 605−614. (20) Werder, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Koumoutsakos, P. On the water-carbon interaction for use in molecular dynamics simulations of graphite and carbon nanotubes. J. Phys. Chem. B 2003, 107, 1345. (21) Merabia, S.; Lombard, J.; Alkurdi, A. Importance of viscoelastic and interface bonding effectsin the thermal boundary conductance of solid-water interfaces. Int. J. Heat Mass Transfer 2016, 100, 287. (22) Kim, B. H.; Beskok, A.; Cagin, T. Molecular dynamics simulations of thermal resistance at the liquid-solid interface. J. Chem. Phys. 2008, 129, 174701. (23) Rafiee, J.; Mi, X.; Gullapalli, H.; Thomas, A. V.; Yavari, F.; Shi, Y.; Ajayan, P. M.; Koratkar, N. A. Wetting transparency of graphene. Nat. Mater. 2012, 11, 217. (24) Leroy, F.; Liu, S.; Zhang, J. Parametrizing Nonbonded Interactions from Wetting Experiments via the Work of Adhesion: H

DOI: 10.1021/acs.langmuir.7b01410 Langmuir XXXX, XXX, XXX−XXX