Influence of Curvature on the Transfer Coefficients for Evaporation and

Mar 18, 2015 - Influence of Curvature on the Transfer Coefficients for Evaporation and Condensation of Lennard-Jones Fluid from Square-Gradient Theory...
0 downloads 8 Views 2MB Size
Article pubs.acs.org/JPCC

Influence of Curvature on the Transfer Coefficients for Evaporation and Condensation of Lennard-Jones Fluid from Square-Gradient Theory and Nonequilibrium Molecular Dynamics Øivind Wilhelmsen,* Thuat T. Trinh, Signe Kjelstrup, and Dick Bedeaux Department of Chemistry, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway ABSTRACT: The influence of curvature on heat and mass transfer across interfaces is often significant in nanoscopic systems. An important example is the initial growth phase of bubbles and droplets in nucleation processes, where they have radii of only a few nanometers. Curvature can be consistently handled within the framework of nonequilibrium thermodynamics by expanding the interface transfer coefficients in the total and Gaussian curvatures. We formulate the coefficients in this expansion analytically and calculate them to second order along the vapor−liquid coexistence line of the Lennard-Jones fluid, which models argon to a good accuracy. To achieve this, squaregradient theory is combined with nonequilibrium molecular dynamics. We discuss how the coefficients depend on the temperature and the truncation value of the interaction potential. The presented coefficients can be used directly to describe heat and mass transfer across interfaces of bubbles/droplets in the Lennard-Jones fluid with cylindrical, spherical, or complex geometries.



INTRODUCTION The field of nanotechnology is rapidly developing and demands knowledge about heat and mass transfer at the nanoscale. Selfassembly of nanoparticles at fluid interfaces (liquid−vapor and liquid−liquid) can for instance enable the preparation of highquality crystals, and nanoparticles at interfaces have been used to modify material stability.1,2 A crucial difference between nano- and macroscopic systems is that interfaces constitute a significant part of the system. Experiments have shown a dependence of the thermal conductivity in nanoparticle suspensions with the size of the particles, where the resistance decreases for smaller particles.3−5 It is very puzzling that nanoporous materials exhibit the opposite behavior with changing pore size.6 It has been hypothesized in the literature that these effects could be explained by the interfacial curvature.7,8 Another example is the irregular shapes of nucleating water clusters, documented for instance by Lengyel et al. in their supersonic expansion experiments.9 Because the nanoclusters deviate much from a spherical geometry, the interface curvature is enormous and will most likely play a key role.8,10 Also for the surface tension, the influence of curvature is poorly understood.11−13 The surface tension of a spherical droplet can be expanded to first order as γ(ξ) = γ0(1 − 2τ/ξ), where ξ is the droplet radius and subscript 0 denotes the flat surface. The sign and magnitude of the first order correction, τ, known as the Tolman length, is still a hot topic of debate.11−14 Interface properties, and how they change under the influence of interface curvature, is of fundamental importance in nanoscience.15−22 We shall in this work investigate the influence of curvature on the interface transfer coefficients of the vapor−liquid interface of the Lennard-Jones (LJ) fluid, which to a good accuracy models argon.23 This gives the © 2015 American Chemical Society

possibility to test the results with both experiments and simulations. Interfaces can be considerable barriers to transfer. For heat transfer, this barrier is quantified by the Kapitza resistance, which has been found in experiments and by nonequilibrium molecular dynamics (NEMD) for a variety of interfaces and substances.16,17,24−26 The interface resistance has also been documented experimentally in evaporation of water, where a an astounding temperature jump of 15.7 K has been measured across the interface.27−29 The transfer of heat and mass across the interface can be described consistently with nonequilibrium thermodynamics:30 1 1 − i = R qqJq′o + R qoμJ To T

(1)

⎛ μo ⎛ 1 μi ⎞ 1⎞ o −⎜ o − i ⎟ + ho⎜ o − i ⎟ = RμoqJq′o + Rμμ J ⎝T T ⎠ T ⎠ ⎝T

(2)

where the superscripts i and o indicate the values of the temperature, T, and the chemical potential, μ, or the enthalpy, h, just inside or outside the interface, respectively. Furthermore, Jq′ is the measurable heat flux, J is the mass flux, and Rij are the interface resistivities, also called interface transfer coefficients. Following Onsager, the matrix of resistivities is symmetric, Rqμ = Rμq. The interface resistivities of the flat interface, Rij,0 are functions of the temperature only. In complex geometries, however, the interface resistivities also depend on the interface curvature. This dependence can be expressed through the two Received: January 21, 2015 Revised: March 16, 2015 Published: March 18, 2015 8160

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C

Figure 1. NEMD simulation cell used in stage 1 of the methodology. The density profile (solid line) is obtained from a NEMD simulation with rc = 2.5, Tcold = 0.85, and Thot = 1.35.

principal curvatures, κ1 and κ2. An equivalent description that is used more commonly is obtained by using the total curvature, H = κ1+κ2, and the Gaussian curvature, K = κ1κ2. In recent work, we proposed the following second-order curvature expansion of the interface resistivities10 R ij(H , K ) = R ij ,0[1 + dijH + νij(dijH )2 + νij̅ dij2K ]

U (r12) = [ULJ(r12) − ULJ(rc)]θ(rc − r12)

(4)

where ULJ(r12) = 4ϵ [(σ/r12)12 − (σ/r12)6] and θ is the Heaviside function. The distance between the particles is r12, the truncation value is rc, the particle diameter is σ, and the minimum of the LJ potential is −ϵ. In the simulations, the variables ϵ, σ, and the mass m of the particles were used to make the variables dimensionless. Stage 1: Local Resistivity of the SGT. In stage 1, we will use NEMD to calibrate the SGT description. To achieve this, we must use a case that can be described simultaneously with NEMD and SGT. We choose the standard procedure used in boundary driven NEMD simulations, where a liquid slab is located in the middle of a rectangular simulation volume with vapor phases on each side, as depicted in Figure 1,35−38 with the dimensions {Lx,Ly,Lz} = {20,20,100} in reduced units. A temperature gradient is imposed by thermostatting the region | z/Lz| < 0.05 to a low temperature, Tcold, and the regions 0.45 < | z/Lz| < 0.5 to a high temperature, Thot. With periodic boundary conditions, the net mass flux is zero. This setup can be handled with most computational programs for NEMD simulations, such as LAMMPS.39 We will now explain how to describe the setup with SGT. The square-gradient model at equilibrium is a minimum of the following functional:

(3)

where i,j = {q,μ}. The length, dij, gives the typical size of a droplet where curvature corrections become important, whereas νij and νij are scalars that are found to be of the order 1. We showed in previous work that the coefficients Rij,0, dij, νij, and νij could be used to describe heat and mass transfer across interfaces of complex geometries, such as an oblate spheroidal droplet (e.g., an M&M candy), a prolate spheroidal bubble (e.g., a rugby ball), or a toroidal bubble (e.g., a doughnut).10 A method was presented that used NEMD to calibrate squaregradient theory (SGT), hence enabling a quantitative mathematical description of the local resistivity across the interface with SGT. The coefficients in eq 3 could then be obtained to a very good accuracy with SGT. We shall in this work complement the method presented previously by formulating the coefficients in the curvature expansion analytically. This allows us to obtain the coefficients to a higher accuracy by solving a set of coupled differential equations. We present also, for the first time, a comprehensive analysis of how the coefficients depend on temperature and on the truncation value of the interaction potential along the vapor−liquid coexistence line of the LJ fluid. The results give the possibility to describe transport of heat and mass across interfaces of nonspherical bubbles/droplets in the LJ fluid/ argon. We expect that the results will give insight into the experiments and simulations dealing with nucleation, where the LJ fluid and argon are popular examples.20,31−34

Ωsgm[T (r), ρ(r)] =



∫ dr ⎢⎣feos (T(r), ρ(r)) ρ(r) − μρ(r) +

⎤ η(T (r), ρ(r)) |∇ρ(r)|2 ⎥ ⎦ 2

(5)

Here, feos is the local Helmholtz energy per particle given by an equation of state (subscript eos), μ is the chemical potential, ρ is the number density, and η is an influence parameter that depends generally on density and temperature. We shall in this work use a temperature- and density-independent influence parameter and describe the local term with the very accurate reference equation of state by Johnson et al.,40 which reproduces results from molecular dynamics (MD) simulations within 2%.40 The equation of state is modified for truncated and shifted potentials according to Smit.41 Using this procedure, we reproduced both the surface tension and the Tolman length from Monte Carlo and MD simulations,13 which shows that SGT correctly accounts for the influence of



THEORY The interface resistivities, Rij, and the coefficients in their curvature expansion are determined in two stages. In the first stage, we calibrate the mathematical formulation from SGT by taking advantage of the atomistic details of NEMD simulations. In the second stage, we use SGT to obtain the coefficients in the curvature expansion up to second order. As example, we consider a fluid of argonlike particles described with the truncated and shifted LJ potential 8161

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C curvature on interfacial properties. We showed that it was crucial to use an accurate equation of state.13 The Euler− Lagrange equations give the following differential equation that defines a stationary state of the square-gradient functional: μ = μeos (T (r), ρ(r)) − η∇2 ρ(r)

where the pressure perpendicular to the interface is represented as

(6)

(19)

Here, μeos is the local contribution from the equation of state. Beyond equilibrium, we assume that thermodynamic functions have the same dependencies on density and its derivatives as at equilibrium. The specific thermodynamic properties (per particle) are42 η f (T (r), ρ(r)) = feos (T (r), ρ(r)) + |∇ρ(r)|2 2ρ(r) (7) u(T (r), ρ(r)) = ueos(T (r), ρ(r)) +

η |∇ρ(r)|2 2ρ(r)

p(T (r), ρ(r)) = peos (T (r), ρ(r)) −

η |∇ρ(r)|2 2

− ηρ(r)∇2 ρ(r) h(T (r), ρ(r)) = heos(T (r), ρ(r)) − η∇2 ρ(r)

On the basis of the Gibbs relation stated in previous work43,44 and the balance equations above, we can derive the entropy production, which according to nonequilibrium thermodynamics gives the following force−flux relation: ⎛ 1 ⎞ ∇⎜ ⎟ = rqq(r)Jq ⎝ T (r) ⎠

(9) (10)

Beyond equilibrium, there is transfer of energy, mass, and momentum through the interface, which can be taken into account by the following balance equations: (11)

∂(ρ v) = −∇·(ρ vv) − ∇p − ∇·γ ∂t

(12)

∂ (ρ(u + 0.5v 2)) = −∇·(ρ v(u + 0.5v 2) + Jq + p v) ∂t (13)

Here, J = ρv is the total mass flux, v is the barycentric velocity, and Jq is the heat flux in the barycentric frame of reference. γ denotes the tension tensor,43 which is represented as

γα , β = η

∂ρ(r) ∂ρ(r) ∂xα ∂xβ

(14)

where α and β denote coordinate components. The pressure tensor is σα , β = psgm δα , β + γα , β

(15)

from which the pressure components perpendicular and parallel to the interface can be obtained. For the flat interface investigated in stage 1, the matrix defining the pressure tensor is diagonal, where two of the diagonal entries are psgm = p∥, and the last entry is the pressure perpendicular to the interface. For the particular case depicted in Figure 1, the mass flux (and mean velocity in z direction) is zero, and eqs 11−13 simplify to 0 = ρvz = J 0=

0=

R ij(ξ) = ℌξ2 ℌ3ξ

∂Jq (17)

∂(p + γzz) ∂z

=

∂p⊥ ∂z



∫−∞ dr

ℌ1 ℌ2 ℌ3

ϕij(r )ext

(21)

where again ξ is the radius. Here, superscript ext means: ϕ(r)ext = ϕ(r) − ϕiΘ(ξ − r) − ϕoΘ(r − ξ). Moreover, ℌ1, ℌ2 , and ℌ3 are Lamé coefficients, and superscript ξ means “evaluated at the dividing surface”. The Lamé coefficients are trivially 1 for a flat surface and functions of the coordinate perpendicular to the interface for curved geometries.8 Equation 21 can also be used to calculate the interface resistivities for curved geometries, with the arguments ϕqq = rqq, ϕqμ = rqq(ho − h), and ϕμμ = rqq(ho −

(16)

∂z

(20)

Both η and peos are available from equilibrium simulations.13 This implies that once an expression for the local resistivity rqq(r) is known it is possible to numerically solve the coupled differential eqs 19 and 20 for a given set of boundary conditions. The solution gives not only the constants, p⊥ and Jq′ but also temperature and density profiles and all properties that depend on these. Hence, rqq(r) represents the missing piece in the nonequilibrium−SGT formulation. In the bulk regions, it is connected to the thermal conductivity λ through rqq(r) = (λT2(r))−1. We will utilize atomistic simulations to obtain the functional form of rqq(r) across the interface in terms of the temperature, density, and density gradients from SGT. In practice, we shall compare the magnitude and position of the temperature jump across the interface in NEMD simulations with those from similar simulations using SGT in order to find a proper mathematical expression for the local resistivity in terms of the temperature, density, and density gradients from SGT. Stage 2: Coefficients in the Curvature Expansion. Both NEMD simulations and SGT give average properties through the interface region that are changing with the coordinate perpendicular to the interface. More practical quantities are therefore the overall interface resistivities, Rij, as defined in eqs 1 and 2. They facilitate a description of heat and mass transfer through the interface with nonequilibrium thermodynamics, where only the value of the temperature and the chemical potentials in the adjacent bulk phases and the temperature and curvature of the interface are needed to calculate the resulting fluxes. The temperature of the interface is close but not equal to the liquid−phase temperature,45 and the curvature is given by the geometry. The link between the local description from SGT and the overall interface resistivities, Rij, is the integral relations. The validity of the integral relations has been verified by NEMD and nonequilibrium SGT.8,36 The integral relations are similar to the Green−Kubo relations in the sense that transport coefficients can be computed using only equilibrium profiles, and for a single-component system, they are represented as46

(8)

∂ρ = −∇·(ρ v) = −∇·J ∂t

∂ 2ρ(z) η ⎛ ∂ρ(z) ⎞ ⎜ ⎟ − ηρ(z) 2 ⎝ ∂z ⎠ ∂z 2 2

p⊥ = peos (T (z), ρ(z)) +

(18) 8162

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C

Figure 2. Two geometries necessary to calculate the coefficients in the curvature expansion of the interface transfer coefficients.

h)2. Hence, all the interface transfer coefficients, Rqq, Rμμ, and Rqμ = Rμq, depend on the local resistivity, rqq. We will describe two approaches to obtain the coefficients in the curvature expansion of the interface resistivities defined in eq 3. Both methods take advantage of the spherical geometry (subscript s, Figure 2a) and the cylindrical geometry (subscript c, Figure 2b), where the curvature expansion becomes ⎡ 2dij dij2(4νij + νij̅ ) ⎤ ⎥ R ij ,s(ξ) = R ij ,0⎢1 + + ⎥⎦ ⎢⎣ ξ ξ2

(22)

⎡ dij dij2νij ⎤ R ij ,c(ξ) = R ij ,0⎢1 + + 2 ⎥ ⎢⎣ ξ ξ ⎥⎦

(23)

procedure was used recently with great success to calculate the coefficients in the curvature expansion of the surface tension;12,47 however, this is the first time the procedure has been formulated for nonequilibrium properties. Even though the quadratic regression procedure has worked excellently for all the cases we have investigated so far, in situations where the curvature dependence of the transport coefficients is weak and the quadratic regression becomes sensitive to numerical inaccuracies, the procedure we will now describe should be preferred. We follow the methodology described by Blokhuis and Bedeaux and expand the density and chemical potential of a spherical droplet in the inverse droplet radius, ξ−1, so that47

The square-gradient model can now be solved for a spherical and a cylindrical geometry from negative curvatures (bubbles) to positive curvatures (droplets) and used together with the integral relations in eq 21 to obtain Rij,s(ξ) and Rij,c(ξ). By performing a quadratic regression of the interface resistivities from negative to positive curvatures for bubbles/droplets of relatively large radii in order to avoid higher order contributions (ξ > 5 in reduced units), we obtain a second-order polynomial in ξ−1 . The regression coefficients in the polynomial correspond to the prefactors in eqs 22 and 23. It turns out that it is sufficient to solve the square-gradient model for a flat interface, a cylindrical geometry, and a spherical geometry to obtain all the relevant coefficients in the curvature expansion, i.e., Rij,0, dij, νij and νij. This approach is most practical and gives the coefficients in the curvature expansion to a good accuracy, as shown previously.10 We emphasize that the interface transfer coefficients of the flat interface and their curvature expansion coefficients are only functions of temperature and are independent of the total volume or number of particles in the simulations, as discussed and showed numerically in previous work.8 An alternative and more elegant procedure that we will pursue in this work is to expand the SGT and the integral relations in ξ−1 and then to define the coefficients analytically as a function of the density expansion functions. A similar

ρs (z) = ρ0 (z) +

μs = μ0 +

1 1 ρ (z) + 2 ρs,2 (z) + ... ξ s,1 ξ

1 1 μ + 2 μs,2 + ... ξ s,1 ξ

(24)

(25)

Here, z = r − ξ, where ξ is the position of the dividing interface. We refer to the functions, ρ0(z), ρs,1(z), etc. as density expansion functions. The above equations can now be used in eq 6, where the Laplace operator is in spherical coordinates. By collecting terms of the same order, one obtains the following set of second-order coupled differential equations:12 ⎛ ∂ 2ρ (z) ⎞ η⎜⎜ 0 2 ⎟⎟ = μeos (ρ0 (z)) − μ0 ⎝ ∂z ⎠

(26)

⎛ ∂ 2ρ (z) ⎞ ⎛ ∂ρ (z) ⎞ ⎛ ∂μ (ρ (z)) ⎞ s,1 ⎟ = −2η⎜ 0 ⎟ + ⎜ eos 0 ⎟ρs,1(z) η⎜⎜ 2 ⎟ ∂ρ ⎠ ⎝ ∂z ⎠ ⎝ ⎝ ∂z ⎠ − μs,1 8163

(27) DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C ⎛ ∂ 2ρ (z) ⎞ ⎛ ∂ρ (z) ⎞ ⎛ ∂ρ (z) ⎞ s,2 ⎟ = −2η⎜⎜ s,1 ⎟⎟ + 2ηz⎜ 0 ⎟ η⎜⎜ 2 ⎟ ⎝ ∂z ⎠ ⎝ ∂z ⎠ ⎝ ∂z ⎠

constraints (eqs 32−34), fully determine the density expansion functions up to second order. These functions can be further used to expand other thermophysical properties in ξ−1, such as the enthalpy and the local resistivity. The last step necessary to define the coefficients in the curvature expansion of the interface resistivities analytically is to expand also the integral relations (eq 21) in ξ−1 for a spherical and a cylindrical geometry. We have carried out this analysis, and the resulting equations defining dij, νij, and νij analytically as functions of the density expansion functions are given in the Appendix.

⎛ ∂μ (ρ (z)) ⎞ ⎟ρs,2 (z) + ⎜ eos 0 ∂ρ ⎠ ⎝ 2 (z) ⎛ ∂ 2μeos (ρ0 (z)) ⎞ ρs,1 ⎜ ⎟⎟ − μs,2 + 2 ⎜⎝ ∂ρ2 ⎠

(28)



where the solution of the first equation gives the density profile across the flat interface. Moreover, μs,1 = 2γ0/Δρ0, where Δρ0 = (ρ0,l − ρ0,g) and ρ0,l and ρ0,g are the liquid and gas coexistence densities, respectively. Furthermore, μs,2 = (−γ0/(Δρ0))(((Δρs,1)/(Δρ0)) + 2τ), where τ is the Tolman length and Δρs,1 = (ρs,1,l − ρs,1,g). We repeat the above procedure for a cylindrical geometry and obtain 2ρc,1(z) = ρs,1(z)

SIMULATION DETAILS We used a collocation method to solve the coupled differential equations for a given set of boundary conditions on a nonuniform grid, defining either the square-gradient model or its curvature expansion (eqs 26−30). The solution was approximated through a third-order polynomial between the grid points. All solutions were obtained to a relative accuracy at every grid point of 10−7 or better. In the NEMD simulations, the whole system was initially held at constant temperature with the Langevin thermostat algorithm48 and equilibrated by performing 2 × 106 steps. The simulation volume ({Lx,Ly,Lz} = {20,20,100}) and the number particles were chosen so that no finite size effects could be observed for the bulk-phase densities or the surface tension. Depending on the temperature, between 10 000 and 12 500 particles were used, where lower temperatures required more particles. Then, the NEMD simulation cell was setup as depicted in Figure 1 and simulated until the system had reached steady state after ∼5 × 106 simulations/steps. After that, another 5 × 106 steps were performed for a production run to collect data at every 50th step. Subsequently, another 106 steps were performed to check convergence, i.e., that temperatures changed less than 1% and the obtained total heat flux changed less than 0.1%. The LAMMPS simulation package was used with a time-step of Δt = 0.002, and periodic boundary conditions were used in all the simulations.39

(29)

⎛ ∂ 2ρ (z) ⎞ ⎛ ∂ρ (z) ⎞ ⎛ ∂ρ (z) ⎞ c,2 ⎟ = −η⎜⎜ c,1 ⎟⎟ + ηz⎜ 0 ⎟ η⎜⎜ 2 ⎟ ⎝ ∂z ⎠ ⎝ ∂z ⎠ ⎝ ∂z ⎠ ⎛ ∂μ (ρ (z)) ⎞ ⎟ρc,2 (z) + ⎜ eos 0 ∂ρ ⎠ ⎝ 2 (z) ⎛ ∂ 2μeos (ρ0 (z)) ⎞ ρc,1 ⎜ ⎟⎟ − μc,2 + 2 ⎜⎝ ∂ρ2 ⎠

(30)

where μc,2 = (−γ0Δρc,1)/2(Δρ0) . The coupled differential equations above determine the density expansion functions for the spherical geometry, ρ0(z), ρs,1(z), ρs,2(z), and for the cylindrical geometry, ρ0(z), ρc,1(z), ρc,2(z), however not uniquely. It can easily be verified that if ρs,1(z) is a solution of eq 27 then so is ρs,1(z) + ζ(∂ρ0(z))/(∂z) where ζ is a constant. Moreover, the number of particles has to be constrained to fix the position of the interface. We have chosen the scaled total amount of particles, Ns = ∫ dz ρ0(z), such that the equimolar surface is located at z = 0. The equimolar dividing surface located at z = 0 gives 2



RESULTS AND DISCUSSION We investigate in this work the interface transfer coefficients of the LJ fluid and the coefficients in their curvature expansion. First, results from stage 1 of the methodology will be discussed, where we obtain the mathematical expression describing the local resistivity of the interface. We then proceed to stage 2 and present the interface transfer coefficients and the coefficients in their curvature expansion, as predicted by SGT. A convenient quantity for describing interface transport, the heat of transfer, will then be discussed, before we exemplify how the curvature expansion can be used in practice to address heat and mass transfer across interfaces of bubbles/droplets with nonspherical shapes. We have investigated the LJ fluid at four truncation values, rc = 2.5, 3, 4, and 6, in reduced units at temperatures ranging from near the triple point to near the critical point. In the figures, however, we report results only for rc = 2.5 and rc = 6 and will refer to results with other truncation values when needed. Moreover, all results will be given in reduced units (LJ units). The results for rc = 6 can to a good approximation be converted to argon by using a well depth of ϵ/kB = 119.8 K, a molecular diameter σ = 0.3405 nm, and a particle mass m = 6.69 × 10−26 kg as parameters,23 where kB is Boltzmann’s constant.10,49 The interface ceases to be a barrier slightly before the critical region, and the characteristic temperature jump



∫−∞ dr [ρ(z)]ext = 0

(31)

with ξ = 0 in the operator indicated by superscript ext, which further leads to the following constraints: ∞

⎛ ∂ρ (z) ⎞ 0 ⎟z 2 ∂z ⎠



∫−∞ dz [ρs,1(z)]ext = ∫−∞ dz ⎜⎝

(32)



∫−∞ dz [ρs,2 (z)]ext = ∞

1 3



⎛ ∂ρ (z) ⎞ 0 ⎟z 3 + ∂z ⎠

∫−∞ dz ⎜⎝





⎛ ∂ρ (z) ⎞ s,1 ⎟⎟z 2 z ∂ ⎝ ⎠

∫−∞ dz ⎜⎜

⎛ ∂ρ (z) ⎞ c,1 ⎟⎟z 2 z ∂ ⎝ ⎠

∫−∞ dz [ρc,2 (z)]ext = 12 ∫−∞ dz ⎜⎜

(33)

(34)

These relations define uniquely the density expansion functions for the equimolar surface. With a different choice of the dividing surface, eqs 32−34 have to be modified. The differential equations (eqs 26−30), together with the 8164

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C

density functional theory descriptions (see Figure 5 in ref 12). This shows that the SGT formulation employed in this work captures the properties of the LJ fluid interface and their dependence on curvature at equilibrium. Figure 4a−d compares density and temperature profiles from SGT and NEMD beyond equilibrium. The figures show that with the proper choice of local resistivity function, rqq (r), SGT gives results that are very close to the time-average results from NEMD. According to kinetic gas theory,51 the origin of the temperature jump across the interface is the interface resistance, located in the so-called “Knudsen layer”. The Knudsen layer has a thickness of approximately one mean free path and is located on the gas side of the interface. In the Knudsen layer, particles having a Maxwellian velocity distribution and the temperature of the gas phase interact with particles coming from the interface, which have a different velocity distribution and temperature.51 It is clear also from NEMD simulations that the interface resistance is on the vapor side of the equimolar surface. An important point elucidated by NEMD simulations, however, is that the interface resistance is only in the very first part of the Knudsen layer. Although the Knudsen layer described in kinetic gas theory has a thickness of approximately one mean free path, the interface resistance from NEMD is positioned in a layer that is only a few molecular diameters thick.51 Contrary to the mean free path in the gas phase that increases at lower temperatures (lower saturation pressures), the layer with the temperature jump in NEMD decreases in size as the temperature decreases (Figure 4c,d). We recommend further investigations to reveal the exact molecular mechanism of the interface resistance of the vapor−liquid interface and to quantify the distance over which the temperature jump occurs. SGT gives continuous profiles through the interface as a function of temperature, density, and derivatives of the density.

across the interface disappears. The critical region is thus not relevant for this discussion and has been omitted from the figures. Stage 1: Local Resistivity of the Interface. The first stage in obtaining the coefficients in the curvature expansion of the interface transfer coefficients was to take advantage of NEMD simulations to calibrate the SGT description. Figure 3

Figure 3. Surface tension of the flat surface from SGT (solid lines), compared to Monte Carlo simulations,50 with rc = 2.5 (circles), rc = 3 (squares), and rc = 6 (diamonds).

shows a comparison of surface tension values for a flat surface obtained from SGT (solid lines) and equilibrium Monte Carlo simulations at different truncation values (circles, squares, and diamonds). As shown in our previous work,13 SGT is very capable of reproducing both the surface tension and the Tolman length from Monte Carlo and MD simulations, in particular at larger truncation values. Moreover, for the coefficients in the curvature expansion of the surface tension, SGT gives nearly identical predictions to more sophisticated

Figure 4. Densities (a and b) and temperature profiles (c and d) from NEMD (solid lines) and the nonequilibrium square-gradient model (red dashed lines) at different truncation values and with different boundary conditions. 8165

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C

Figure 5. Interface resistivities, Rij, and the first-order coefficients in their curvature expansion, dij, with rc = 6 (solid lines) and rc = 2.5 (dashed lines).

region near the critical point at low truncation values (rc = 2.5). There is a small deviation between the density profiles given by the red dashed line and the solid line at the bottom of Figure 4a, near z/Lz = 0.2, because of the tail-corrections of the LJ equation of state used at low truncation values.13 We find that the term α(T)ρ(r)−2|∇ρ(r)|2 describes the resistivity well at low temperatures and that the term β(T)ρ(r)−4|∇ρ(r)|2 dominates at high temperatures. This is captured by strongly temperature-dependent prefactors. For rc 10 = 2.5, we find α = 1.15T−5 r and β = 0.12Tr , with Tr being the reduced temperature. For rc = 3,4, and 6, we find α = 0.29T−5 r and β = 0.20T10 r . The latter expression works very well for all truncation values we tested above rc = 3. We believe that the change in α(Tr) and β(Tr) from rc = 2.5 to rc ≥ 3 is due to the inaccuracies in the equation of state at lower cutoff values. With these simple local resistivity functions as input, the nonequilibrium square-gradient model gives results very close to NEMD. The bottom solid lines in Figure 4c,d show larger fluctuations in the temperature of the gas phase in the NEMD simulations because of the limited number of particles, whereas the nonequilibrium square-gradient model captures the mean behavior well.

To capture the behavior of the interface barrier in the framework of SGT, we search for a function that gives a peak at the interface located closer to the vapor phase. We found that the following expression represents a simple function, which reproduces to a good accuracy the time-average results from NEMD for the LJ fluid: rqq =

⎛ α (T ) β (T ) ⎞ 1 +⎜ 2 + 4 ⎟|∇ρ(r)|2 2 ρ (r) ⎠ λ(T , ρ)T (r) ⎝ ρ (r)

(35)

where λ is the thermal conductivity of the fluid and we have used the model by Bugel and Galliero.52 The second contribution to eq 35 gives a peak in the resistivity at the interface, which is shifted toward the vapor phase compared to the equimolar surface because of the ρ−2 and ρ−4 prefactors. The local resistivity in eq 35 can also be used for curved interfaces. Figure 4 presents the time-average results from the righthand side of the NEMD simulation cell depicted in Figure 1 (solid lines), where the two thermostated regions are also shown. The density profiles from NEMD (solid lines) are reproduced to a good accuracy by the nonequilibrium squaregradient model (red dashed lines), except in the interface 8166

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C

Figure 6. Second-order coefficients in the curvature expansion of the interface resistivities, νij and νij, with rc = 6 (solid lines) and rc = 2.5 (dashed lines).

expansion of the interface resistivities, dij, have the same dimension as the first-order correction of the surface tension, the Tolman length. From recent MD simulations, the Tolman length at similar conditions was −0.1 in reduced units, which is about 45 times smaller than, for instance, dqq with rc = 6.14 It is clear that the interface resistivities depend much more strongly on curvature than the surface tension. The linear coefficients are all within −5 < dij < −1.4 and change little with temperature, in particular for rc = 2.5. By comparing the solid and dashed lines in Figure 5, it is evident that the absolute magnitude of the first-order coefficients increases with the truncation distance. The Tolman length exhibits a similar behavior.12 It is crucial to use an accurate equation of state in the analysis to obtain the interface properties quantitatively, as discussed in previous work.13 Figure 6 reports the second-order coefficients in the curvature expansion of the interface resistivities, νij and νij. They are dimensionless, unlike the second-order curvature corrections of the surface tension, the rigidity constants. In previous work,10 we formulated a novel curvature expansion of the interface transfer coefficients. We assumed that only one length scale was relevant for the expansion to the second order.

Equation 35 is found to work well for vapor−liquid interfaces. However, we expect the second contribution on the right-hand side to depend on the type of interface and to be different for liquid−liquid or liquid−solid interfaces because the interface resistance here is from a different physical mechanism. A similar approach, however, can be carried out for these systems to identify rqq(T,ρ(r)). Stage 2: Coefficients in the Curvature Expansion. In stage 1, we captured the local interface resistance mathematically in the SGT framework. We can now calculate the overall interface resistivities and the coefficients in their curvature expansion. The resulting interface resistivities, Rqq, Rqμ, and Rμμ, for the LJ fluid are plotted in Figure 5 as a function of temperature and for two truncation values. The figures show that all the interface resistivities decrease nearly exponentially with temperature, similar to previous work.49,53 Moreover, all coefficients increase in magnitude by a factor of 2−3 from rc = 2.5 (dashed lines) to rc = 6 (solid lines). This behavior was also observed, at least for the cross-coefficients between heat and mass, in a work by Ge et al. when they increased the truncation value of the LJ potential in NEMD simulations from 1.7 to 2.5 in reduced units.54 The linear coefficients in the curvature 8167

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C

Figure 7. Interface resistivities divided by the interface resistivity of the flat interface, calculated with square-gradient theory and the integral relations (solid lines) and with the curvature expansion (red dashed lines) for spherical and cylindrical bubbles and droplets. The results are for rc = 2.5 and T = 0.7.

Because νij and νij are of order unity for both truncation values and all temperatures investigated, only one length scale is relevant for the curvature expansion to the second order. Equation 3 is therefore the appropriate expression to use for the interface resistivities. In contrast to the first-order coefficients, dij, the second-order coefficients in the curvature expansion, νij and νij, decrease in magnitude with increasing truncation distance. The curvature expansion is only to the second order, and has a limited range of validity. The integral relations defined in eq 21, however, can also be used to calculate the interface transfer coefficients for smaller bubbles/droplets, on the basis of only the equilibrium density profiles through the interface for bubbles/droplets with spherical and cylindrical geometries, as illustrated by the solid lines in Figure 7. The interface transfer coefficients of spherical/cylindrical bubbles and droplets are reproduced well by the curvature expansion at T = 0.7, to a

total curvature (H = κ1 + κ2) of about 0.2 (red dashed lines). Because the second-order expansion gives convex parabolic functions with a minimum at a positive total curvature, it extrapolates better for bubbles than for droplets. The integral relations (solid lines, Figure 7) show that all the interface transfer coefficients decrease monotonically with increasing total curvature. The curvature expansion will only exhibit this behavior until the minimum of the parabola defined by the quadratic curvature expansion is encountered. The range of validity of the curvature expansion varies with temperature, truncation value, and whether the curvature is positive or negative (compare the left- and right-hand sides of Figure 7a-f). Heat of Transfer. Transfer across interfaces can also be formulated using the quantity called the heat of transfer, defined as 8168

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C

Figure 8. Heat of transfer, (superscript g refers to the gas-phase). (a) and the linear term in the corresponding curvature expansion (b), for rc = 2.5 (dashed lines) and rc = 6 (solid lines). The heat of transfer for T = 0.7 and rc = 2.5 for a spherical geometry, as a function of the total curvature (c), from the integral relations (solid line) and from eq 39 (red dashed line).

Figure 9. Local thermal interface resistivity normalized by the thermal resistivity of the flat interface, for a prolate spheroidal droplet (a) and an oblate spheroidal bubble (b), obtained with the curvature expansion and coefficients for rc = 2.5 and T = 0.7.

⎛ J′o ⎞ R qoμ q =− q*o = ⎜⎜ ⎟⎟ R qq ⎝ J ⎠ΔT = 0, J = 0

It is apparent from Figure 8c why the heat of transfer is a convenient quantity. Contrary to the interface resistivities (Figures 7a−f), the heat of transfer changes approximately linearly with the total curvature (H), and higher order coefficients are very small. Consequently, the range of validity of the expansion of the heat of transfer is longer than that of the cross coefficient, Roqμ. The solid line in Figure 8c has been obtained by using the integral relations and can to a very good approximation be reproduced by taking advantage of the following linear curvature expansion of the heat of transfer:

(36)

The heat of transfer is an alternative and convenient way of defining the coupling between heat and mass transfer. The advantage of using the heat of transfer is that it has the same units as the enthalpy and can then easily be compared to, for instance, the vaporization enthalpy. The constitutive relations across the interface can then be reformulated as 1 1 *o ′o o − i = R qq(Jq − q J ) T T

q*o(H , K ) ≈ 1 + (dμq − dqq)H = 1 + dhtH q *o

(37)

0

⎛ μo ⎛ 1 μi ⎞ 1⎞ o −⎜ o − i ⎟ + ho⎜ o − i ⎟ = ( −R qqq*o)Jq′o + Rμμ J ⎝T T ⎠ T ⎠ ⎝T

(39)

The linear coefficient in the curvature expansion of the heat of transfer, dht, is between 1.3 and 2.3 in reduced units, as shown in Figure 8b.

(38) 8169

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C Interface Resistivities of Spheroidal Bubbles and Droplets. The importance of the coefficients presented in Figures 5−7 becomes evident when one wants to describe the growth of bubbles/droplets in the LJ fluid/argon. Because the interface transfer coefficients depend strongly on the curvature, the interface resistivities of nanosized bubbles/droplets are significantly different from the corresponding values at the flat surface (Figure 7). The interface resistivities are higher for bubbles (negative curvature) and lower for droplets (positive curvature). Numerous simulations have shown that the bubble/ droplet formed during a nucleation process generally deviates from a spherical geometry.55−57 However, until now, it has been impossible to deal with the additional complexity of nonspherical shapes, even for the simple LJ fluid. The curvature expansion in eq 3 gives the possibility to deal with heat and mass transfer across interfaces also in complex geometries. As shown in previous work, there are, in principal, no restrictions to which geometries and shapes can be addressed. The LJ fluid is the benchmark fluid for nucleation simulations, and we believe that the presented coefficients have an important role in further understanding of these. We shall demonstrate how the curvature expansion can be used in practice, via the straightforward example of spheroidal bubbles/droplets resembling the shape often observed in nucleation simulations. Figure 9 shows the normalized thermal interface resistivities of droplets and bubbles with oblate and prolate spheroidal geometries. The spheroidal radii are 26 and 15 in reduced units. It is evident that curvature affects droplets and bubbles very differently. Although the thermal interface resistivity of the droplet is considerably smaller than that of a flat interface (Figure 9a), the interface resistance of the bubbles is at some locations almost twice that value (Figure 9b). These differences can be explained by the ratio ℌ1/(ℌ2 ℌ3) in the operator used in the integral relations (eq 21), which increases toward the center of the bubble/droplet. A peak in the local resistivity at the inside of the interface therefore increases the resistance with curvature (bubbles), whereas it decreases with a peak at the outside (droplets). The same behavior should be expected for all interfaces, e.g., in hydrates, crystals, or biological systems. The main consequence of the large curvature dependence of the interface resistivities is that for nanosized bubbles/droplets with nonspherical shapes the resistance to transfer will be highly local. In practice, this means that fluctuations in shape enable the bubble/droplet to enhance growth in some directions. For droplets, the growth will be enhanced at locations with high curvature, which favors a nonspherical shape in the initial growth phase. This agrees with the irregularly shaped water clusters observed in in recent experiments and simulations dealing with supersonic expansions in water.9 The LJ fluid is a much simpler fluid than water, with a well-defined interaction potential. We believe that the coefficients presented in this work, together with simulations with the LJ fluid and experiments with argon, will in the future give further insight into the shape and growth rate of nucleating bubbles and droplets.

obtain these coefficients. They are large compared to similar coefficients for the surface tension, which shows that interface transfer coefficients are strongly influenced by curvature. Their dependence on temperature and truncation distance was analyzed. We used the curvature expansion to calculate the interface resistances of droplets and bubbles having spheroidal shapes. For droplets, transport was enhanced with curvature, whereas the opposite effect was seen for bubbles. The curvature expansion and the corresponding coefficients can be used to give insight into the multitude of experiments and simulations dealing with vapor−liquid nucleation in the LJ fluid/argon.



APPENDIX

Curvature Expansion of the Interface Resistivities

The curvature expansion functions of the enthalpy and local resistivity are hs(z) = h0(z) +

1 1 hs,1(z) + 2 hs,2(z) + ... ξ ξ

rqq ,s(z) = rqq(0)(z) +

(40)

1 1 rqq(s,1)(z) + 2 rqq(s,2)(z) + ... ξ ξ (41)

and equivalently for the cylindrical geometry, where z = r − ξ and ξ is the radius of the droplet. Once the density expansion functions are known, i.e., ρ0(z), ρs,1(z), ρs,2(z), and ρc,1(z), ρc,2(z), for a spherical and cylindrical geometry, respectively, the curvature expansion functions of the enthalpy, h, can be calculated through h0(z) = μ0 + heos(ρ0 (z)) − μeos (ρ0 (z))

(42)

⎛⎛ ∂h (ρ (z)) ⎞ ⎛ ∂μ (ρ (z)) ⎞⎞ eos 0 ⎟⎟⎟ ⎟ − ⎜ eos 0 hs,1(z) = μs,1 + ρs,1(z)⎜⎜⎜ ∂ ρ ∂ρ ⎠⎠ ⎠ ⎝ ⎝⎝ (43)

⎛⎛ ∂h (ρ (z)) ⎞ ⎛ ∂μ (ρ (z)) ⎞⎞ eos 0 ⎟⎟⎟ ⎟ − ⎜ eos 0 hs,2(z) = μs,2 + ρs,2 (z)⎜⎜⎜ ∂ ρ ∂ρ ⎠⎠ ⎠ ⎝ ⎝⎝ ρs,1(z)2 ⎛⎛ ∂ 2heos(ρ0 (z)) ⎞ ⎛ ∂ 2μeos (ρ0 (z)) ⎞⎞ ⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟⎟⎟ + 2 ⎜⎝⎜⎝ ∂ρ2 ∂ρ2 ⎠ ⎝ ⎠⎠ (44)

and the curvature expansion functions of the local resistivity, rqq, are 2 ⎛ α β ⎞⎛ ∂ρ0 (z) ⎞ ⎟ ⎟ ⎜ rqq(0)(z) = rqq , loc(ρ0 (z)) + ⎜⎜ + ⎟ 2 ρ0 (z)4 ⎠⎝ ∂z ⎠ ⎝ ρ0 (z)

(45)

⎛ ∂rqq , loc(ρ (z)) ⎞ 0 ⎟⎟ rqq(s,1)(z) = ρs,1(z)⎜⎜ ∂ ρ ⎝ ⎠



⎛ α β ⎞⎛ ∂ρ0 (z) ⎞⎛ ∂ρs,1(z) ⎞ ⎟⎟⎜ ⎟⎟ ⎟⎜⎜ + 2⎜⎜ + 2 ρ0 (z)4 ⎠⎝ ∂z ⎠⎝ ∂z ⎠ ⎝ ρ0 (z)

CONCLUSIONS We have in this work obtained the interface transfer coefficients of the Lennard-Jones fluid and the corresponding coefficients in a second-order curvature expansion, defined in terms of the total and Gaussian curvatures. Square-gradient theory and nonequilibrium molecular dynamics have been employed to

⎛ α ⎛ ∂ρ (z) ⎞2 2β ⎞ ⎟⎟ρ (z)⎜ 0 ⎟ − 2⎜⎜ + 3 ρ0 (z)5 ⎠ s,1 ⎝ ∂z ⎠ ⎝ ρ0 (z) (46) 8170

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C ρs,1(z)2 ⎛ ∂ 2rqq , loc(ρ0 (z)) ⎞ ⎜ ⎟ + ρ (z ) ⎟ s,2 2 ⎜⎝ ∂ρ2 ⎠

rqq(s,2)(z) =

dq2μR qμ ,0νqμ ∞

=

⎛ ∂rqq , loc(ρ (z)) ⎞ ⎛ ∂ρ (z) ⎞2 ⎛ 0 ⎜⎜ ⎟⎟ + ⎜ 0 ⎟ ⎜⎜ρs,1(z)2 ∂ρ ⎝ ⎠ ⎝ ∂z ⎠ ⎝

∫−∞ dz [rqq(0)Δohc,2 + rqq(c,1)Δohc,1 + rqq(c,2)Δoh0]ext ∞

+−

∫−∞ dz z[rqq(0)Δohc,1 + rqq(c,1)Δoh0] ∞

⎛ 3α 10β ⎞⎟ ⎜⎜ + ⎟ − 2ρs,2 (z) 4 ρ0 (z)6 ⎠ ⎝ ρ0 (z)

+

(53)

and the interface resistivity to mass transfer, Rμμ

⎞ ⎛ α β ⎞ 2β ⎞⎟ ⎛ α ⎜⎜ ⎟ ⎟⎟ + + ⎜⎜ + 3 5 ⎟⎟ 2 ρ0 (z)4 ⎠ ρ0 (z) ⎠⎠ ⎝ ρ0 (z) ⎝ ρ0 (z)

dμμRμμ ,0 =



1 2

∫−∞ dz [2rqq(0)Δohs,1Δoh0 + rqq(s,1)Δoh02]ext



∫−∞ dz z[rqq(0)Δoh02]ext



⎛⎛ ∂ρ (z) ⎞2 ⎞⎛ ∂ρ (z) ⎞⎞ ⎛ ⎜⎜ s,1 ⎟ + 2⎜ ∂ρ0 (z) ⎟⎜ s,2 ⎟⎟ ⎜ ⎟ ⎜⎜⎝ ∂z ⎟⎠ ⎝ ∂z ⎠⎝ ∂z ⎠⎟⎠ ⎝

(54)

2 dμμ Rμμ ,0(4νμμ + νμμ ̅ )

⎛ α 2β ⎞⎛ ∂ρ0 (z) ⎞ ⎟⎟⎜ ⎟ − 4ρs,1(z)⎜⎜ + 3 ρ0 (z)5 ⎠⎝ ∂z ⎠ ⎝ ρ0 (z) ⎛ ∂ρ (z) ⎞ ⎜⎜ s,1 ⎟⎟ ⎝ ∂z ⎠

∫−∞ dz z 2[rqq(0)Δoh0]ext



=

∫−∞ dz [rqq(0)Δohs,12 + 2rqq(0)Δoh0Δohs,2

ext

+ 2Δoh0Δohs,1rqq(s,1) + rqq(s,2)Δoh02] ∞

(47)

−2

∫−∞ dz z[2rqq(0)Δoh0Δohs,1 + rqq(s,1)Δoh02]ext

+3

∫−∞ dz z 2[rqq(0)Δoh02]ext



For the cylindrical geometry, the coefficients follow the same formulas. For ease of notation, we omit (z) in the following formulas to indicate that the variables depend on z. The last step necessary to define the coefficients in the curvature expansion of the interface resistivities analytically is to expand the integral relations (eq 21) in ξ−1 for a spherical and cylindrical geometry. This gives the following expressions for Rqq dqqR qq ,0

1 = 2



∫−∞ dz [rqq(s,1)]

ext





2 dμμ Rμμ ,0νμμ ∞

=

ext



∫−∞ dz z[rqq(0)]

ext

=

∫−∞ dz [rqq(s,2)] ∞

+3



∫−∞ dz z[rqq(s,1)]

ext

ext



∫−∞ dz z 2[rqq(0)]ext

(50)



(51)



ACKNOWLEDGMENTS We would like to thank E. M. Blokhuis for clarifying the role of the dividing surface. Ø.W. would like to thank the Faculty of Natural Sciences and Technology at the Norwegian University of Science and Technology for the Ph.D. grant.



1 2

∫−∞ dz [rqq(0)Δohs,1 + rqq(s,1)Δoh0]ext



∫−∞ dz z[rqq(0)Δoh0]ext



dq2μR qμ ,0(4νqμ ∞ =

AUTHOR INFORMATION

The authors declare no competing financial interest.

the cross coefficient, Rqμ dqμR qμ ,0 =

(56)

o

Notes





+

∫−∞ dz z 2[rqq(0)Δoh02]ext

*E-mail: [email protected].

(49)

∫−∞ dz [rqq(c,2)]ext − ∫−∞ dz z[rqq(c,1)]ext

2 dqq R qq ,0νqq =

+

Corresponding Author

∫−∞ dz z [rqq(0)] 2

∫−∞ dz z[2rqq(0)Δoh0Δohc,1 + rqq(c,1)Δoh02]ext

where Δoh = h − h.



−2





2 dqq R qq ,0(4νqq + νqq̅ ) ext

∫−∞ dz [rqq(0)Δohc,12 + 2rqq(0)Δoh0vhc,2 + 2Δoh0Δohc,1 rqq(c,1) + rqq(c,2)Δoh02]

(48)



(55)

(1) Lin, Y.; Böker, A.; Skaff, H.; Cookson, D.; Dinsmore, A. D.; Emrick, T.; Russel, T. P. Nanoparticle Assembly at Fluid Interfaces: Structure and Dynamics. Langmuir 2005, 21, 191. (2) Bresme, F.; Oettel, M. Nanoparticles at Fluid Interfaces. J. Phys.: Condens. Matter 2007, 19, 413101. (3) Keblinski, P.; Phillpot, S. R.; Choi, S. U. S.; Eastman, J. A. Mechanisms of Heat Flow in Suspensions of Nano-Sized Particles (Nanofluids). Int. J. Heat Mass Transfer 2002, 45, 855. (4) Patel, H. E.; Das, S. K.; Sundararajan, T.; Nair, A. S.; George, B.; Pradeep, T. Thermal Conductivities of Naked and Monolayer Protected Metal Nanoparticle Based Nanofluids: Manifestation of

+ νq̅ μ)

∫−∞ dz [rqq(0)Δohs,2 + rqq(s,1)Δohs,1 + rqq(s,2)Δoh0]ext ∞

−2

∫−∞ dz z[rqq(0)Δohs,1 + rqq(s,1)Δoh0]ext

+3

∫−∞ dz z 2[rqq(0)Δoh0]ext

REFERENCES



(52) 8171

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C Anomalous Enhancement and Chemical Effects. Appl. Phys. Lett. 2003, 83, 2931. (5) Hu, M.; Poulikakos, D.; Grigoropoulos, C. P.; Pan, H. Recrystallization of Picosecond Laser-Melted ZnO Nanoparticles in a Liquid: A Molecular Dynamics Study. J. Chem. Phys. 2010, 132, 164504. (6) Bera, C.; Mingo, N.; Volz, S. Marked Effects on Alloying on the Thermal Conductivity of Nanoporous Materials. Phys. Rev. Lett. 2010, 104, 115502. (7) Lervik, A.; Bresme, F.; Kjelstrup, S. Heat Transfer in Soft Nanoscale Interfaces: the Influence of Interface Curvature. Soft Matter 2009, 5, 2407. (8) Wilhelmsen, Ø.; Bedeaux, D.; Kjelstrup, S. Heat and Mass Transfer through Interfaces of Nanosized Bubbles/Droplets: the Influence of Interface Curvature. Phys. Chem. Chem. Phys. 2014, 16, 10573. (9) Lengyel, J.; Pysanenko, A.; Poterya, V.; Slavíček, P.; Fárník, M.; Kočišek, J.; Fedor, J. Irregular Shapes of Water Clusters Generated in Supersonic Expansions. Phys. Rev. Lett. 2014, 112, 113401. (10) Wilhelmsen, Ø.; Trinh, T.; Kjelstrup, S.; van der Erp, T. S.; Bedeaux, D. Heat and Mass Transfer across Interfaces in Complex Nanogeometries. Phys. Rev. Lett. 2015, 114, 065901. (11) Lei, Y. A.; Bykov, T.; Yoo, S.; Zeng, X. C. The Tolman Length: Is it Positive or Negative? J. Am. Chem. Soc. 2005, 127, 15346. (12) Blokhuis, E. M.; van Giessen, A. E. Density Functional Theory of a Curved Liquid-Vapour Interface: Evaluation of the Rigidity Constants. J. Phys.: Condens. Matter 2013, 25, 225003. (13) Wilhelmsen, Ø.; Bedeaux, D.; Reguera, D. Tolman Length and Rigidity Constants of the Lennard-Jones fluid. J. Chem. Phys. 2015, 142, 064706. (14) van Giessen, A. E.; Blokhuis, E. M. Direct Determination of the Tolman Length from the Bulk Pressures of Liquid Drops via Molecular Dynamics Simulations. J. Chem. Phys. 2009, 131, 164705. (15) Kapitza, P. L. The Study of Heat Transfer in Helium II. J. Phys. (USSR) 1941, 4, 181. (16) Patel, H. A.; Garde, S.; Keblinski, P. Thermal Resistance of Nanoscopic Liquid-Liquid Interfaces: Dependence on Chemistry and Molecular Architecture. Nano Lett. 2005, 5, 2225. (17) Ge, Z.; Cahill, D. G.; Braun, P. V. Thermal Conductance of Hydrophilic and Hydrophobic Interfaces. Phys. Rev. Lett. 2006, 96, 186101. (18) Smith, J. D.; Cappa, C. D.; Drisdell, W. S.; Cohen, R. C.; Saykally, R. J. Raman Thermometry of Free Evaporation from Liquid Water Droplets. J. Am. Chem. Soc. 2006, 128, 12892. (19) Fladerer, A.; Strey, R. Homogeneous Nucleation and Droplet Growth in Supersaturated Argon Vapor: The Cryogenic Nucleation Pulse Chamber. J. Chem. Phys. 2006, 124, 164710. (20) Iland, K.; Wölk, J.; Strey, R.; Kashchiev, D. Argon Nucleation in a Cryogenic Pulse Chamber. J. Chem. Phys. 2007, 127, 154506. (21) Liang, Z.; Sasikumar, K.; Keblinski, P. Thermal Transport across a Substrate-Thin-Film Interface: Effects of Film Thickness and Surface Roughness. Phys. Rev. Lett. 2014, 113, 065901. (22) Davies, G. B.; Krüger, T.; Coveney, P. V.; Harting, J.; Bresme, F. Assembling Ellipsoidal Particles at Fluid Interfaces using Switchable Dipolar Capillary Interactions. Adv. Mater. 2014, 26, 6715. (23) Rowley, L. A.; Nicholson, D.; Parsonage, N. G. Monte Carlo Grand Canonical Ensemble Calculation in a Gas-Liquid Transition Region for 12-6 Argon. J. Chem. Phys. 1975, 17, 401. (24) Eastman, J. A.; Phillpot, S. R.; Choi, S. U. S.; Keblinski, P. Thermal Transport in Nanofluids. Annu. Rev. Mater. Res. 2004, 34, 219. (25) Kim, B. H.; Beskok, A.; Cagin, T. Molecular Dynamics Simulations of Thermal Resistance at the Liquid-Solid Interface. J. Chem. Phys. 2008, 129, 174701. (26) Cheng, C.; Fan, W.; Cao, J.; Ryu, S.-G.; Ji, J.; Grigoropoulos, C. P.; Wu, J. Heat Transfer across the Interface between Nanoscale Solids and Gas. ACS Nano 2011, 5, 10102.

(27) Duan, F.; Badam, V. K.; Durst, F.; Ward, C. A. Thermocapillary Transport of Energy during Water Evaporation. Phys. Rev. E 2005, 72, 056303. (28) Badam, V. K.; Kumar, V.; Durst, F.; Danov, K. Experimental and Theoretical Investigations on Interfacial Temperature Jumps during Evaporation. Exp. Therm. Fluid Sci. 2007, 32, 276. (29) Duan, F.; Ward, C. A.; Badam, V. K.; Durst, F. Role of Molecular Phonons and Interfacial-Temperature Discontinuities in Water Evaporation. Phys. Rev. E 2008, 78, 041130. (30) Kjelstrup, S.; Bedeaux, D. Non-Equilibrium Thermodynamics of Heterogeneous Systems; World Scientific: Hackensack, NJ, 2008. (31) Chen, B.; Siepmann, J. I.; Kwang, J. O.; Klein, M. L. Aggregation-Volume-Bias Monte Carlo Simulations of Vapor-Liquid Nucleation Barriers for Lennard-Jonesium. J. Chem. Phys. 2001, 115, 10903. (32) Wedekind, J.; Wölk, J.; Reguera, D.; Strey, R. Nucleation Rate Isotherms of Argon from Molecular Dynamics Simulations. J. Chem. Phys. 2007, 127, 154515. (33) Hale, B. N.; Thomason, M. Scaled Vapor-to-Liquid Nucleation in a Lennard-Jones System. Phys. Rev. Lett. 2010, 105, 046101. (34) Diemand, J.; Angílil, R.; Tanaka, K. K.; Tanaka, H. Large Scale Molecular Dynamics Simulations of Homogeneous Nucleation. J. Chem. Phys. 2013, 139, 074309. (35) Røsjorde, A.; Kjelstrup, S.; Bedeaux, D.; Hafskjold, B. Nonequilibrium Molecular Dynamics Simulations of Steady-State Heat and Mass Transport in Condensation. II. Transfer Coefficients. J. Colloid Interface Sci. 2001, 240, 355. (36) Simon, J. M.; Bedeaux, D.; Kjelstrup, S.; Xu, J.; Johannessen, E. Interface Film Resistivities for Heat and Mass Transfers-Integral Relations Verified by Non-equilibrium Molecular Dynamics. J. Phys. Chem. B 2006, 110, 18528. (37) Bresme, F.; Armstrong, J. Note: Local Thermal Conductivities from Boundary Driven Non-Equilibrium Molecular Dynamics Simulations. J. Chem. Phys. 2014, 140, 016102. (38) Trinh, T. T.; Vlugt, T. J. H.; Kjelstrup, S. Thermal Conductivity of Carbon Dioxide from Non-Equilibrium Molecular Dynamics: A Systematic Study of Several Common Force Fields. J. Chem. Phys. 2014, 141, 134504. (39) PlimptonS.Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comp. Phys., 1995, 117, 1−19. LAMMPS: Large-Scale Atomic/Molecular Massively Parallel Simulator, Sandia National Lab, 2007 version. http://lammps.sandia.gov (40) Johnson, J. K.; Zollweg, J. A.; Gubbins, K. E. The Lennard-Jones Equation of State Revisited. Mol. Phys. 1993, 78, 591. (41) Smit, B. Phase Diagrams of Lennard-Jones Fluids. J. Chem. Phys. 1992, 96, 8639. (42) Wilhelmsen, Ø.; Bedeaux, D.; Kjelstrup, S.; Reguera, D. Thermodynamic Stability of Nanoscale Bubbles/Droplets: The Square Gradient Theory and the Capillary Approach. J. Chem. Phys. 2014, 140, 024704. (43) Glavatskiy, K. S. Multicomponent Interfacial Transport as Described by the Square Gradient Model; Evaporation and Condensation. Ph.D. thesis, Norwegian University of Science and Technology, 2009. (44) Magnanelli, E.; Wilhelmsen, Ø.; Bedeaux, D.; Kjelstrup, S. Extending the Nonequilibrium Square-Gradient Model with Temperature-Dependent Influence Parameters. Phys. Rev. E 2014, 90, 032402. (45) Johannessen, E.; Bedeaux, D. The Nonequilibrium van der Waals Square Gradient Model. (II). Local Equilibrium of the Gibbs Surface. Phys. A (Amsterdam, Neth.) 2003, 330, 354. (46) Glavatskiy, K. S.; Bedeaux, D. Resistances for Heat and Mass Transfer Through a Liquid-Vapor Interface in a Binary Mixture. J. Chem. Phys. 2010, 133, 234501. (47) Blokhuis, E. M.; Bedeaux, D. Van der Waals Theory of Curved Surfaces. Mol. Phys. 1993, 80, 705. (48) Frenkel, D.; Smit, B. Understanding Molecular Simulation, From Algorithms to Applications; Academic Press: San Diego, CA, 2002. (49) Inzoli, I.; Kjelstrup, S.; Bedeaux, D.; Simon, J. M. Transfer Coefficients for the Liquid-Vapour Interface of a Two-Component Mixture. Chem. Eng. Sci. 2011, 66, 4533. 8172

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173

Article

The Journal of Physical Chemistry C (50) Grosfils, P.; Lutsko, J. F. Dependence of the Liquid-Vapor Surface Tension on the Range of Interaction: a Test of the Law of Corresponding States. J. Chem. Phys. 2009, 130, 054703. (51) Bond, M.; Struchtrup, H. Mean Evaporation and Condensation Coefficients based on Energy Dependent Condensation Probability. Phys. Rev. E 2004, 70, 061605. (52) Bugel, M.; Galliero, G. Thermal Conductivity of the LennardJones Fluid: An Empirical Correlation. Chem. Phys. 2008, 352, 249. (53) Johannessen, E.; Bedeaux, D. The Nonequilibrium van der Waals Square Gradient Model. (III). Heat and Mass Transfer Coefficients. Physica A 2004, 336, 252. (54) Ge, J.; Kjelstrup, S.; Bedeaux, D.; Simon, J. M.; Rousseau, B. Transfer Coefficients for Evaporation of a System with a LennardJones Long-Range Spline Potential. Phys. Rev. E 2007, 75, 061604. (55) Holyst, R.; Litniewski, M. Heat Transfer at the Nanoscale: Evaporation of Nanodroplets. Phys. Rev. Lett. 2008, 100, 055701. (56) Wang, Z.-J.; Valeriani, C.; Frenkel, D. Homogeneous Bubble Nucleation Driven by Local Hot Spots: A Molecular Dynamics Study. J. Phys. Chem. B 2009, 113, 3776. (57) Joswiak, M. N.; Duff, N.; Doherty, M. F.; Peters, B. SizeDependent Surface Free Energy and Tolman-Corrected Droplet Nucleation of TIP4P/2005 Water. J. Phys. Chem. Lett. 2013, 4, 4267.

8173

DOI: 10.1021/acs.jpcc.5b00615 J. Phys. Chem. C 2015, 119, 8160−8173