Predictions of homogeneous nucleation rates for $n$-alkanes accounting for the diffuse phase interface and capillary waves

TL;DR

This study improves the prediction of homogeneous nucleation rates for n-alkanes by incorporating capillary wave effects and a density gradient theory, addressing discrepancies with experimental data.

Contribution

It introduces a correction for capillary wave effects on surface tension within the density gradient theory for n-alkanes, enhancing nucleation rate predictions.

Findings

Capillary waves decrease surface tension and nucleation rates.

Critical clusters are nearly free of capillary wave effects.

Theoretical corrections have minor impact on temperature dependence.

Abstract

Homogeneous droplet nucleation has been studied for almost a century but has not yet been fully understood. In this work, we used the density gradient theory (DGT) and considered the influence of capillary waves (CW) on the predicted size-dependent surface tensions and nucleation rates for selected $n$-alkanes. The DGT model was completed by an equation of state (EoS) based on the perturbed-chain statistical associating fluid theory (PC-SAFT) and compared to the classical nucleation theory and the Peng--Robinson EoS. It was found that the critical clusters are practically free of CW because they are so small that even the smallest CW wavelengths do not fit into their finite dimensions. The CW contribute to the entropy of the system and thus decrease the surface tension. A correction for the effect of CW on the surface tension is presented. The effect of the different EoSs is relatively small because by a fortuitous coincidence their predictions are similar in the relevant range of critical cluster sizes. The difference of the DGT predictions to the classical nucleation theory computations is important but not decisive. Of the effects investigated, the most pronounced is the suppression of the CW which causes a sizable decrease of the predicted nucleation rates. The major difference between experimental nucleation rate data and theoretical predictions remains in the temperature dependence. For normal alkanes, this discrepancy is much stronger than observed, e.g., for water. The theoretical corrections developed here have a minor influence on the temperature dependency. We provide empirical equations correcting the predicted nucleation rates to values comparable with experiments.

Tables

Table 1. Table 1: Coefficients of Eq. ( 36 ).

DGT with PC-SAFT:	A	B	C
$n$-heptane	-110.7738	62.1559	-6.1818
$n$-octane	-94.01211	53.1095	-6.1818
$n$-nonane	-91.89308	52.5203	-6.1818
$n$-decane	-100.6096	58.3962	-6.1818
CNT with PC-SAFT:	A	B	C
$n$-heptane	-97.9379	53.8308	-3.7779
$n$-octane	-82.0400	45.0167	-3.7779
$n$-nonane	-79.5862	44.0546	-3.7779
$n$-decane	-88.9715	50.0613	-3.7779
DGT with PR:	A	B	C
$n$-heptane	-134.7050	80.0453	-12.0051
$n$-octane	-119.4497	73.1246	-12.0051
$n$-nonane	-119.5625	75.3177	-12.0051
$n$-decane	-130.1146	83.2105	-12.0051
CNT with PR:	A	B	C
$n$-heptane	-124.5077	72.4331	-8.7391
$n$-octane	-109.2277	65.3503	-8.7391
$n$-nonane	-109.5674	67.5294	-8.7391
$n$-decane	-119.1090	74.7644	-8.7391

Table 2. Table 2: Parameters of PR EoS

	$T_c\left(K\right)$	$p_c\left(\mathrm{bar}\right)$	$\omega (-)$
$n$-heptane	540.13	27.36	0.3365
$n$-octane	569.32	24.97	0.3765
$n$-nonane	594.55	22.81	0.4140
$n$-decane	617.70	21.03	0.4495

Table 3. Table 3: Table of parameters for the fit of the experimental surface tension

	$A$ (mN/m)	$B$ (mN/m)	$C$ (mN/m)
$n$-heptane	54.1778	-0.7586	3.9897
$n$-octane	56.5399	-10.4928	8.4723
$n$-nonane	53.9000	-	-
$n$-decane	53.6000	-	-

Equations

$$S=\exp \left(\frac{{\mu }_{\mathrm{V}}-{\mu }_{\infty }}{k_{\mathrm{B}}T}\right),$$

$${\mu }_{\mathrm{L}}(n_{\mathrm{c}})={\mu }_{\mathrm{V}}.$$

$$J={\left[\sum _{n=1}^{\infty }\frac{1}{f_n{C}_n^{\mathrm{eq}}}\right]}^{-1},$$

$$f_n=\nu A_n,$$

$$\nu =C_1\sqrt{\frac{k_{\mathrm{B}}T}{2\pi M}},$$

$$C_{n}^{\mathrm{eq}}=\frac{1}{\vartheta}\exp\Big{(}-\frac{\Delta\Omega_{n}}{k_{\mathrm{B}}T}\Big{)},$$

$$J=\frac{1}{\vartheta }A(n_{\mathrm{c}})\nu \mathcal{Z}\exp \left(-\frac{\Delta \Omega (n_{\mathrm{c}})}{k_{\mathrm{B}}T}\right),$$

$$\mathcal{Z}=\sqrt{-\frac{\Delta {\Omega }^{\prime \prime }(n_{\mathrm{c}})}{2\pi k_{\mathrm{B}}T}}\simeq \frac{2}{{\rho }_{\mathrm{L},\infty }}\sqrt{\frac{{\sigma }_{\infty }}{k_{\mathrm{B}}T}}.$$

$$J=\frac{{\rho }_{\mathrm{V}}{\rho }_{V\infty }}{{\rho }_{L\infty }}\sqrt{\frac{2{\sigma }_{\infty }}{\pi M}}\exp \left(-\frac{\Delta \Omega }{k_{\mathrm{B}}T}\right).$$

$$\Delta \Omega =\frac{1}{3}{A}_{\mathrm{s}}\sigma ,$$

$$\Delta p=\frac{2\sigma }{r_{\mathrm{s}}},$$

$$\Delta \Omega =\frac{16\pi }{3}\frac{{\sigma }^3}{\Delta p^2}.$$

$$f(\rho ,\nabla \rho ,{\nabla }^2\rho )=f^0(\rho )+c_1(\rho ,T){\nabla }^2\rho +c_2(\rho ,T)(\nabla \rho {)}^2+\dots ,$$

$$f(\rho )=f^0(\rho )+c(\rho ,T)\cdot (\nabla \rho {)}^2+\dots .$$

$$\Delta \Omega [\rho (r)]={\int }_0^{\infty }\left[\Delta {\omega }^0(\rho )+\frac{1}{2}c{\left(\frac{\mathrm{d}\rho }{\mathrm{d}r}\right)}^2\right]4\pi r^2\mathrm{d}r,$$

$$\Delta {\omega }^0(\rho )=f^0(\rho )-\rho {\mu }_{\mathrm{V}}+p_{\mathrm{V}}.$$

$$\frac{{\mathrm{d}}^2\rho }{\mathrm{d}{r}^2}+\frac{2}{r}\frac{\mathrm{d}\rho }{\mathrm{d}r}=\frac{1}{c}\Delta \mu (\rho ).$$

$$\rho (r\to \infty )={\rho }_{\mathrm{V}},\frac{\mathrm{d}\rho }{\mathrm{d}r}(0)=0.$$

$$\Delta N={\int }_0^{\infty }\left[\rho (r)-{\rho }_{\mathrm{V}}\right]4\pi r^2\mathrm{d}r.$$

$${\sigma }_{\infty }={\int }_{-\infty }^{\infty }\left[\Delta {\omega }^0(\rho )+\frac{1}{2}c{\left(\frac{\mathrm{d}\rho }{\mathrm{d}z}\right)}^2\right]\mathrm{d}z,$$

$${\sigma }_{\exp }={\sigma }_{\mathrm{bare}}-\frac{3}{8\pi }{k}_{\mathrm{B}}T{q}_{\max }^2,$$

$$q_{\max }=\frac{1}{2.64{\xi }^{+}},$$

$$q_{\max }=\frac{1}{5.17{\xi }^{-}}.$$

$${\xi }^{\pm }={\xi }_0^{\pm }{\left|1-T/T_{\mathrm{c}}\right|}^{-\nu },$$

$$\frac{{\xi }_0^{+}}{{\xi }_0^{-}}\simeq 1.96,$$

$$\frac{{\sigma }_0({\xi }_0^{+}{)}^2}{k_{\mathrm{B}}{T}_{\mathrm{c}}}\simeq 0.39,$$

$$\sigma ={\sigma }_0(1-T/T_{\mathrm{c}}{)}^{2\nu }$$

$${\sigma }_{\mathrm{bare}}={\sigma }_{\exp }\left(1+\frac{3}{8\pi }\frac{T}{T_{\mathrm{c}}}\frac{1}{2.72}\right).$$

$${\lambda }_{\max }=\pi r_{\mathrm{s}},$$

$$q_{\min }=\frac{2\pi }{{\lambda }_{\max }}=\frac{2}{r_{\mathrm{s}}}.$$

Full text

I Introduction

Understanding the homogeneous droplet nucleation is important in many natural and industrial processes such as formation of secondary aerosols in the atmosphere, formation of water droplets in the steam turbines, or nucleation during the gas-cleaning procedures, e.g., in processing of natural gas or in technologies for carbon capture and storage.

Despite many attempts that resolved partial problems of nucleation, there is no complete theory which would give quantitatively correct predictions. The basic approach to nucleation is through the classical nucleation theory (CNT), developed by Becker and DöringBecker and Döring (1935) and extended by Zeldovich.Zeldovich (1942, 1943) The concepts of CNT have been described in several textbooks.Vehkamäki (2006); Kaschiev (2000); Kalikmanov (2013) The CNT models a cluster as a macroscopic liquid droplet with a step-wise phase interface. However, the real thickness of the interface is comparable with the size of the droplets during the nucleation (in the range of nanometers). The goal of this work is to overcome this limitation by considering a diffuse phase interface (a smooth density profile) and by taking the capillary waves into account. The progress of nucleation theory is backed by experimental results. The state of the art in the experimental nucleation science has been recently reviewed by Wyslouzil and Wölk.Wyslouzil and Wölk (2016) In this work, experimental results are used for comparisons without considering experimental circumstances.

Clearly, the step-wise density profile is a crude approximation. A smooth (or diffuse) change of the density between the two phases describes reality better. This approach is considered in the classical density functional theory (DFT) and its approximation, the density gradient theory (DGT). The DGT was first developed in pioneering work of van der Waals,van der Waals (1893); Rowlinson (1979) and later thoroughly elaborated by Cahn and Hilliard.Cahn and Hilliard (1958, 1959) The classical DFT was proposed by EbnerEbner, Saam, and Stroud (1976) and developed by EvansEvans (1979), based on the quantum DFT. Laaksonen and OxtobyLaaksonen and Oxtoby (1995) studied nucleation of binary mixtures using Lennard–Jones fluids. They used DFT for ideal and non-ideal mixtures not obeying the Berthelot mixing rules. They found large discrepancies between DFT and CNT, especially for the non-ideal mixtures. Wilemski and LiWilemski and Li (2004) investigated DGT and mean field DFT of droplets and bubbles near the spinodal. They compared their results to an approximation of DGT near spinodal developed by Cahn and Hilliard Cahn and Hilliard (1959) and found a qualitative agreement. They address limitations of the mean field treatment of real fluids near spinodal since this region also coincides with region of large fluctuations which are not properly represented in the mean-field theory. LutskoLutsko (2008) examined nucleation using the minimal energy free path method for the Lennard–Jones fluid. He incorporated the nudged-elastic band technique which adds fictitious elastic forces to have collection of density profiles evenly spaced. The advantage of this technique is its relatively simple implementation and robustness. Lutsko got good agreement with simulations of the structure of the cluster, such as the excess number of molecules of the critical cluster or density profiles. More recently,Lutsko (2011) Lutsko incorporated DGT with density profiles approximated as piecewise linear functions for the Lennard–Jones fluid. Unlike the present work, he computed the needed influence parameter using the direct correlation function. Simplified density profiles allowed direct comparison with the CNT. Density profiles, excess number of molecules in the cluster and excess free energy computed by Lutsko agreed well with the simulation data. McGrath et al.McGrath et al. (2010) examined the nucleation of argon using DFT and Monte Carlo simulations and compared them to CNT. They found rather large differences for four various potentials, as reported by other studies. DFT results agreed more with Monte Carlo simulations than with CNT. They concluded that one of the reasons why CNT fails to describe nucleation correctly is in improper molecular potentials. Obeidat et al.Obeidat et al. (2010) combined DGT and SAFT (statistical associating fluid theory)Chapman et al. (1989) EoS (equation of state) for methanol and compared it to CNT. To get the results of DGT, they used a finite difference scheme, unlike in this work. As will be described in Sec. VI, we found this method impractical for the systems considered here. McGrath et al. reported that the DGT gives better nucleation rates than the CNT. While the smooth-interface models are clearly closer to reality than CNT and may provide a significant improvement for specific systems, we found that the considering solely the smooth-interface effect is not sufficient to provide quantitatively correct predictions for $n$-alkanes.

The diffuse-interface approaches of DFT and DGT generate an “intrinsic” profile of density through the phase interface. This profile corresponds to averaging over short space of time scales. Besides that, the interface is also disturbed by thermally driven mesoscopic undulations, the capillary waves (CW). The physics of CW was described by Buff, Lovett and Stillinger.Buff, Lovett, and F. H. Stillinger (1965) CW were later investigated by many researchers, e.g., by KayserKayser (1986) who calculated an increase in interfacial tension caused by cutting off CW longer than a certain wavelength, Beysens and RobertBeysens and Robert (1987) and MeunierMeunier and Langevin (1982); Meunier (1987) who combined smoothly changing density approach with CW. Meunier’s approach, developed for planar phase interfaces, is used in this work to estimate the effect of CW on the surface tension for droplets (clusters) relevant to nucleation. Mecke, Dietrich and NapiórkowskiNapiórkowski and Dietrich (1993); Mecke and Dietrich (1999, 2005) derived an effective Hamiltonian combining the effect of DFT and CW and surface tension depending on the wavenumber $q$. Their complex CW model differs significantly from the traditionalBuff, Lovett, and F. H. Stillinger (1965) approach. Conclusions of their theory led to a rather wide discussion.Doerr et al. (1999); Fradin et al. (2000); Lin et al. (2005); Tarazona, Checa, and Chacón (2007); Hiester (2008) Chacón, Tarazona and AlejandreChacón and Tarazona (2003); Tarazona and Chacón (2004); Chacón, Tarazona, and Alejandre (2006) developed a method to determine a definition of the instantaneous dividing surface in molecular simulations applied to atomic fluids and water. They bridged the approaches of CW and interpretation of x-ray reflectivity experiments making it an alternative model of CW. Boltachev and BaidakovBoltachev and Baidakov (2006) provided a new empirical model of CW fitted to ellipsometry data for argon. Their approach provides an alternative for determining the effect of CW on surface tension and thickness of the phase interface.

In this work we address the problem of homogeneous droplet nucleation with DGT using PC (perturbed chain) SAFTGross and Sadowski (2000, 2001) EoS and considering effect of CW. For comparison, we add theoretical results also for CNT and cubic Peng–Robinson (PR)Peng and Robinson (1976) EoS.

This article is organized as follows. In Sec. II we consider some basic relations for nucleation computations and mention their weak points. We connect the nucleation rate to the work needed to form the droplet in supersaturated vapor. In Sec. III we provide mathematical results of DGT to determine the density profiles that correspond to the work of formation of a cluster. In Sec. IV we introduce the PC-SAFT EoS. In Sec. V we characterize the impact of CW on surface tension and show how to eliminate CW from the input of DGT. In Sec. VI we explain mathematical methods for the numerical solution defined in Sec. III. In Sec. VII we provide the resulting nucleation rates of four $n$-alkanes; namely $n$-heptane, $n$-octane, $n$-nonane and $n$-decane, and compare them to the experimental data. We reevaluate the data to determine the temperature dependence of the supersaturation and number of molecules in the critical cluster. We examine impact of CW on the surface tension and the nucleation rate. In the last section, Sec. VIII, we draw conclusions about the nucleation rates, temperature dependence, and influence of CW.

II Nucleation rates

Nucleation of droplets in a homogeneous vapor can occur if the vapor density $\rho_{\mathrm{V}}$ is greater than the density of saturated vapor $\rho_{\mathrm{V\infty}}$ at a given temperature $T$. In further text, subscript $\mathrm{\infty}$ refers to the infinite curvature radius of the planar phase interface under the saturation conditions. The departure of the system from saturation is expressed in terms of supersaturation $S$ defined as

[TABLE]

where $\mu_{\mathrm{V}}$ is the chemical potential of the vapor phase, $\mu_{\infty}$ is the chemical potential at saturation, which is equal for liquid and vapor phases, and $k_{\mathrm{B}}$ is the Boltzmann constant. If $S>1$, the vapor becomes supersaturated. Supersaturated vapor is in a metastable state, i.e., it is stable with respect to small fluctuations but unstable with respect to large fluctuations, particularly in the form of clusters of the newly forming phase.Rowlinson and Widom (1982); Kalikmanov (2013) Homogeneous droplet nucleation is the process of formation of liquid droplets in the supersaturated vapor phase free of aerosol particles, impurities, or solid walls. These droplets can be viewed as clusters of $n$ molecules, referred to as $n$-mers. Microscopically, these clusters grow or shrink if a monomer joins or leaves the cluster. Condensation and evaporation of dimers is less probable except for strongly associating vapors. Of all the $n$-mers, the most important for nucleation is the so-called critical cluster with $n_{\mathrm{c}}$ molecules which has the same probability to grow as to shrink.Becker and Döring (1935) In terms of thermodynamics, this condition can be interpreted as an equality of chemical potential in the liquid core of the critical cluster and the chemical potential of the supersaturated vapor,

[TABLE]

The number of droplets formed in a unit of volume per a unit of time is called the nucleation rate, $J$. Considering the chain of individual growth-shrink processes, it is possible to derive an expression for the nucleation rate in the form of a converging infinite sum,Vehkamäki (2006); Kaschiev (2000); Kalikmanov (2013)

[TABLE]

where $f_{n}$ is the condensation rate, i.e., the number of monomers impacting on the surface of $n$-mer per unit of time; the sticking probability is taken as unity. $C_{n}^{\mathrm{eq}}$ is the equilibrium concentration of $n$-mers at given $T$ and $\mu_{\mathrm{V}}$ in a system which is constrained in the way that the growth over certain size $n_{\mathrm{max}}$ is prohibited. If $n_{\mathrm{max}}$ is substantially larger than the critical size $n_{\mathrm{c}}$, its magnitude does not play any role. When no interactions between large clusters are considered, the limit of maximal size of $n$-mer can be placed to infinity,Vehkamäki (2006) as adopted in Eq. (3). Constrained equilibrium concentrations enter the nucleation theory via the condition of detailed balanceVehkamäki (2006); Kalikmanov (2013) stating that in the equilibrium all microscopic processes must have their time-mirrored counterparts which are equally probable. This condition allows for expression of the evaporation rate in terms of the condensation rate. The latter is given as

[TABLE]

where $A_{n}$ is the surface area of $n$-mer and $\nu$ is the impingement rate of monomers per surface area. This is given by the kinetic theory of gases asVehkamäki (2006); Kalikmanov (2013); Landau and Lifshitz (1980)

[TABLE]

where $C_{1}$ is the concentration of monomers and $M$ is the molar mass. The concentration $C_{n}^{\mathrm{eq}}$ can be expressed in terms of the work of formation $\Delta\Omega_{n}$ of $n$-mer asVehkamäki (2006); Kalikmanov (2013)

[TABLE]

where $\vartheta$ is the volume available for the translation of the cluster and attributes to the translational degrees of freedom of the cluster. The Boltzmann factor $\exp(-\Delta\Omega_{n}/k_{\mathrm{B}}T)$ can be interpreted as a probability of observing at least one $n$-mer in a region of volume $\vartheta$ in a (constrained) equilibrium system. The original choice by Becker and DöringBecker and Döring (1935) was $1/\vartheta=\rho_{\mathrm{V}}$, which, however, has been shown to affect the dependence of cluster concentrations on supersaturation in a theoretically unacceptable manner.Barrett (1997); Girshick and Chiu (1990); Hrubý (2004); Reguera and Reiss (2004) In this work, we adopt $1/\vartheta=\rho_{\mathrm{V}\infty}$, corresponding to the correction by Courtney.Courtney (1961) The CNT modification by Girshick and ChiuGirshick and Chiu (1990) corresponds to $1/\vartheta=\rho_{\mathrm{V}\infty}e^{\theta}$, where $\theta$ is the so-called dimensionless surface tension (variable consisting of the surface tension, liquid density, and thermal energy $k_{\mathrm{B}}T$). The reasoning behind the Girshick-Chiu modification is that the expression for equilibrium concentration of clusters, Eq. (6), should hold down to monomers – provided that the dimensionless surface tension is size-independent. Here, however, we consider the size dependency and, therefore, this correction was not adopted.

The sum in Eq. (3) can be approximated by an integral. The integrand is a function of $n$ showing a sharp peak at $n_{\mathrm{c}}$, which is due to the exponential part (the work of formation) in Eq. (6). Compared to that, the condensation rate varies modestly with $n$ and, therefore, can be approximated as $f_{n}\simeq f_{n_{\mathrm{c}}}$.Vehkamäki (2006); Kaschiev (2000); Kalikmanov (2013) In a more refined treatment, $1/f_{n}$ can be approximated as a linear function of $n$ in the neighborhood of $n_{\mathrm{c}}$. The proportional part cancels out upon the integration. The peaked function $1/C_{n}^{\mathrm{eq}}$ is then replaced by a Gaussian function and, integrated, yielding toVehkamäki (2006); Kaschiev (2000); Kalikmanov (2013)

[TABLE]

where $\mathcal{Z}$ is the Zeldovich factor defined as

[TABLE]

In Eq. (8), $\Delta\Omega^{\prime\prime}$ is the second derivative of the work of formation with respect to the number of molecules in the cluster, evaluated for the critical cluster, $\Delta\Omega^{\prime\prime}(n_{\mathrm{c}})=\partial^{2}\Delta\Omega/\partial n^{2}(n_{\mathrm{c}})$, $\rho_{\mathrm{L,\infty}}$ is the liquid density, and $\sigma_{\infty}$ is the surface tension for the saturated state. The rightmost expression of the Zeldovich factor is an approximation based on CNT. Although different expressions will be used for $\Delta\Omega$, the Zeldovich factor can be kept in the classical form because it is rather insensitive to the modifications considered. With Eqs. (5) and (8), and using $\vartheta$ corresponding to Courtney’s correction,Courtney (1961) nucleation rate given by Eq. (7) becomes

[TABLE]

The work of formation $\Delta\Omega$ of the critical cluster can be expressed in terms of the Gibbsian thermodynamics of phase interfaces as

[TABLE]

where $A_{\mathrm{s}}=4\pi r_{\mathrm{s}}^{2}$ is the area of the surface of tension of the cluster and $r_{\mathrm{s}}$ is its radius. The surface of tension is defined by postulating the Young–Laplace equation in its standard form,Rowlinson and Widom (1982); Vehkamäki (2006)

[TABLE]

where $\Delta p=p_{\mathrm{L}}-p_{\mathrm{V}}$ is the Laplace pressure. The work of formation is connected to the surface tension via

[TABLE]

We note that the Laplace pressure $\Delta p$ for the critical cluster is related to supersaturation given by Eq. (1) through the equality defined in Eq. (2) and EoS employed for calculation of chemical potential and density.

The work of formation for the CNT is obtained by substituting the flat-interface surface tension $\sigma_{\infty}$ into Eq. (12). Inversely, if $\Delta\Omega$ is known, this equation can be used to compute the (size-dependent) surface tension for the critical cluster. Typically, further simplifications are employed in CNT. In particular, the compressibility of the liquid is neglected and the gas phase is assumed as an ideal gas. These simplifications were not employed in this work; the CNT work of formation was computed from appropriate EoS (PR EoS and PC-SAFT EoS).

III Density gradient theory (DGT)

The DGT operates with a density $\rho$ smoothly changing throughout the phase interface. The local Helmholtz energy density is assumedCahn and Hilliard (1958) as a function of density and invariants of the spatial derivatives of density – the gradient squared $(\nabla\rho)^{2}$, and the Laplacian $\nabla^{2}\rho$. The Helmholtz energy is expanded in terms of these invariants around the Helmholtz energy of the homogeneous fluid $f^{0}(\rho)$ which can be computed from EoS in the region between saturated vapor and saturated liquid. We note that this region involves metastable vapor, metastable liquid, and the region between spinodals, which is unstable in the homogeneous case, but inside the phase interface it is stabilized by the density gradient. Retaining only linear terms in a Taylor expansion with respect to variables $(\nabla\rho)^{2}$ and $\nabla^{2}\rho$, the Helmholtz energy density readsCahn and Hilliard (1958, 1959); Baidakov et al. (2002)

[TABLE]

where $c_{1}$ and $c_{2}$ are Taylor coefficients. Using the Green (divergence) theorem, the second derivative $\nabla^{2}$ in Eq. (13) can be substituted to the first derivative which leads toCahn and Hilliard (1958)

[TABLE]

In Eq. (14), the second term accounts for the inhomogeneity induced by the interface. Coefficient $c=-\partial c_{1}/\partial\rho+c_{2}$ is called the influence parameterCahn and Hilliard (1958) and can be related to the direct correlation function.Davis (1996) In this work, it will be evaluated from the experimental data (as discussed below) and assumed as density-independent.

The work of formation of a droplet according to DGT can be expressed as a volume integral consisting of homogeneous and gradient parts,Cahn and Hilliard (1958); Davis (1996)

[TABLE]

where $r$ is the radial coordinate (the distance from the center of the droplet). Term $\Delta\omega^{0}$ is the difference of the grand-potential density of the inhomogeneous system with respect to the homogeneous one at the same chemical potential (equal to the negative value of the vapor pressure $p_{\mathrm{V}}$),

[TABLE]

We note that the work of formation, Eq. (15), is not a mere function of the density, but a functional since its value depends on the whole density profile $\rho(r)$ and not just on the local value. The work of formation has a saddle point for the density profile corresponding to the so-called critical cluster; the work of formation of the critical cluster is minimal with respect to all the properties in the functional space except for density profile variations changing the size of the cluster (shifting the phase interface closer to or further from the center), for which it is maximal.Rowlinson and Widom (1982); Barrett (1997) Setting the functional derivative of Eq. (15) equal to zero as a condition for the saddle point leads to an Euler–Lagrange equationCahn and Hilliard (1959)

[TABLE]

In Eq. (17), $\Delta\mu=\mu^{0}(\rho)-\mu_{\mathrm{V}}$ is the difference of chemical potential of a hypothetical homogeneous fluid at a given local density $\rho$ and the actual chemical potential $\mu_{\mathrm{V}}$. In the critical cluster, the chemical potential is everywhere equal to chemical potential of the homogeneous vapor phase. To solve the Euler–Lagrange equation, two boundary conditions are needed,Cahn and Hilliard (1959)

[TABLE]

Other important outputs of the DGT calculations are the excess number of molecules $\Delta N$ and the surface tension of the saturated state $\sigma_{\infty}$ (corresponding to the planar phase interface). The excess number of molecules is a difference of the number of molecules in a system containing a cluster to the number of molecules in a system of the same volume containing homogeneous vapor only. It can be computed as

[TABLE]

The surface tension for a planar phase interface separating saturated liquid and saturated vapor can be obtained in a relation similar to Eq. (15),

[TABLE]

where $z$ is the coordinate perpendicular to the phase interface. Surface tension can also be obtained from experimental data; therefore, Eq. (20) can be used for computation of the (density-independent) influence parameter $c$.Hrubý, Labetski, and van Dongen (2007)

IV PC-SAFT equation of state

In this work, DGT was combined with the PC-SAFTGross and Sadowski (2000, 2001) EoS, which belongs to the family of modern SAFT-type equations of state being developed since 1990’s.Chapman et al. (1989); Müller and Gubbins (2001); Tan, Adidharma, and Radosz (2008) The SAFT-type equations produce a qualitatively correct isotherm in the two-phase region, the so-called van der Waals loop. On the other hand, they provide remarkably good results for the liquid phase compared, e.g., to the cubic equations. As a consequence, the SAFT-type equations represent a convenient tool for modeling of vapor-liquid phase interfaces and for modeling nucleation of droplets and bubbles using DFT and DGT.

The SAFT-type EoSs are defined in the form of the Helmholtz energy, which is given as a sum of the ideal gas part $F_{\mathrm{id}}$ evaluated from the isobaric heat capacity of the ideal gas and the residual part $F_{\mathrm{res}}$ defined by the SAFT terms. An important advantage of the SAFT-type EoSs is that the residual part of the Helmholtz energy is expressed as a sum of individual contributions accounting for various types of intermolecular interactions,Tan, Adidharma, and Radosz (2008) e.g., the van der Waals attractions, Coulombic forces, or hydrogen bonds.Müller and Gubbins (2001); Tan, Adidharma, and Radosz (2008) In PC-SAFT, the residual part consists of the hard chain contribution $F_{\mathrm{hc}}$ representing the reference fluid and the perturbation contribution $F_{\mathrm{disp}}$. The hard chain contribution and the perturbation contribution are defined by a set of three molecular parameters: the segment number $m$, the segment diameter $\sigma$, and the energy parameter $\epsilon/k_{\mathrm{B}}$. Values of these parameters are usually correlated to the saturation pressure and the liquid density.Gross and Sadowski (2001) Due to their reasonable physical basis and possibility of considering various contributions in $F_{\mathrm{res}}$, the SAFT-type equations can be applied to a large variety of substances ranging from relatively simple components such as gases, $n$-alkanes or halogenated hydrocarbonsVinš and Vacek (2009); Vinš, Hrubý, and Planková (2013) to polar substances,Vinš, Hrubý, and Planková (2012); Schäfer, Sadowski, and Enders (2014) polymers,Gross et al. (2003) and associating fluids.Tan, Adidharma, and Radosz (2008); Grenner et al. (2007) Consequently, these equations of state are becoming popular both in scientific and engineering applications.Müller and Gubbins (2001); Tan, Adidharma, and Radosz (2008) Molecular parameters $m$, $\sigma$ and $\epsilon/k_{\mathrm{B}}$ for selected $n$-alkanes used in this study were taken from the original work of Gross and Sadowski.Gross and Sadowski (2001)

V Capillary waves

Thermal motion of molecules causes that the interface is not a smooth surface but it is rather disturbed by the capillary waves.Buff, Lovett, and F. H. Stillinger (1965) We will use some results of the so-called mode-coupling theory developed by MeunierMeunier (1987) which considers effect of the broadening of the phase interface due to the CW on the surface tension for a planar phase interface.

We define a “bare” surface tension $\sigma_{\mathrm{bare}}$ as the surface tension of a phase interface which is forced to be perfectly planar, cleared of CW. The bare surface tension still accounts for some thermal motion of molecules, which, however, is only correlated over a lengthscale comparable to the bulk correlation length. The surface tension determined in macroscopic experiments $\sigma_{\mathrm{exp}}$ includes the effect of CW. Allowing CW adds some disorder to the system which leads to an increase of the surface excess entropy. Consequently, multiplied by temperature, the CW reduce the surface tension asMeunier (1987)

[TABLE]

where $q_{\mathrm{max}}$ is the upper cutoff of the wavenumbers considered. $q_{\mathrm{max}}$ corresponds to CW with the smallest wavelength that are not already accounted for in DGT.

MeunierMeunier (1987) estimated $q_{\mathrm{max}}$ based on the mode-coupling theory as

[TABLE]

where we corrected the original value Meunier (1987) 2.55 to the new value $2.64$ according to the recent results of the critical scaling theory.Anisimov and Bertrand (2010) $\xi^{+}$ is the correlation length of the supercritical fluid. Since we are dealing with subcritical temperatures, it is more consistent (but numerically equivalent) to reformulate Eq. (22) in terms of $\xi^{-}$, the correlation length of the subcritical fluid,

[TABLE]

Both correlation lengths obey a scaling law

[TABLE]

where $\nu=0.63$ is a critical exponent. Critical amplitudes $\xi^{\pm}_{0}$ in Eq. (24) are related as

[TABLE]

which explains the mathematical equivalence of Eqs. (22) and (23). Critical amplitude $\xi^{+}_{0}$ can be obtained from another critical scaling ratio

[TABLE]

where $T_{\mathrm{c}}$ is the critical temperature. The critical amplitude $\sigma_{0}$ for the surface tension has been obtained by fitting scaling relation

[TABLE]

to the experimental surface tension data (for the list of data see Appendix B).

Using Eqs. (23) to (27), the bare surface tension in Eq. (21) can be related to the experimental value as

[TABLE]

This is equivalent to the result of Gross,Gross (2009) aside from the updated numerical value in the denominator. Effect of the CW is strongest at the critical point, where the bare surface tension is higher by $4.42\,\%$ than the experimental value, and it remains significant down to the triple point temperature. We note that Eq. (28) is universal, i.e., substance independent.

Eq. (21) assumes an infinitely large system. In a finite system, the maximal wavelength of the CW is limited by the system’s linear dimension. Correspondingly, a lower cutoff of the wavenumber $q_{\mathrm{min}}$ should be considered. This is also the case of the nanodroplets relevant to nucleation. We estimate the maximal wavelength as a half of the circumference of the droplet,

[TABLE]

where radius of the droplet is represented by the radius of the surface of tension, $r_{\mathrm{s}}$. This estimate corresponds to the oblate-prolate ellipsoid oscillations, which is the lowest mode of CW on a sphere. The lower cutoff of the wavenumber is then

[TABLE]

The critical cluster radius depends on the temperature $T$ and the supersaturation $S$. We evaluated both cutoffs for the experimental nucleation data considered here for the comparison with theoretical predictions (more in Sec. VII). This was done in order to determine the range of the CW wavenumbers relevant for the critical clusters.

The upper cutoff $q_{\mathrm{max}}$ was directly computed for a given temperature $T$ rounded to a certain value for each series (every individual experiment differed in the temperature somewhat). To evaluate the lower cutoff $q_{\mathrm{min}}$ using Eq. (30), the radius $r_{\mathrm{s}}$ was obtained by interpolating the data computed by DGT for the experimental supersaturation $S_{\mathrm{exp}}$. For every rounded temperature we got several cutoffs $q_{\mathrm{min}}$ depending on the supersaturation of the data point; for clarity we considered only the minimal $q_{\mathrm{min}}$ for given rounded temperature so that the range $(q_{\mathrm{min}},q_{\mathrm{max}})$ was the widest possible. Fig. 1 shows both cutoffs for the experimental data computed using Eqs. (23) and (30). The seemingly outlaying $q_{\mathrm{min}}$ point in the upper right corner of Fig. 1 corresponds to data by Wagner and StreyWagner and Strey (1984) at higher supersaturations and thus smaller critical clusters than the other data. It can be seen that $q_{\mathrm{min}}>q_{\mathrm{max}}$ in all cases. Therefore, the range $(q_{\mathrm{min}},q_{\mathrm{max}})$ is an empty interval; this means that the critical clusters are too small to accommodate any CW. All thermal fluctuations are bulk-like density fluctuations included in the continuous density profile model. Consequently, the bare spherical phase interface obtained by DFT or DGT is an appropriate representation for the critical clusters.

VI Numerical method

Eq. (17) with boundary conditions given by Eq. (18) form a boundary value problem. This simply looking problem has several difficulties. First, the density profile near the gaseous phase has a sharp, corner-like, shape; at the gas-phase side, its slope changes abruptly. This “stiff” behavior is caused by the strongly non-linear right-hand side of the differential equation, Eq. (17). The second difficulty is that for large droplets the density profile in the interior of the droplet is almost perfectly flat. The difference of the density in the droplet center $\rho(0)$ to the density $\rho_{\mathrm{L}}$, following from the equality of chemical potentials, Eq. (2), approaches zero rapidly as the droplet becomes large compared to the thickness of the phase interface. This situation occurs for droplets much larger than the typical critical clusters, so that the resulting uncertainty does not influence the nucleation rate computations presented here. However, it represents an obstacle in studying the asymptotic behavior of the size-dependent surface tension when the droplet radius approaches infinity.

In a previous work,Hrubý, Labetski, and van Dongen (2007) the boundary value problem was solved using a matlabMathworks.com (2010) procedure bvp4. This procedure uses a finite difference scheme: the original interval is divided into subintervals, where the right-hand side of differential equation (17) is integrated. An initial approximation of the solution is required. In many general cases, the initial guess can be very crude. However, in case of Eq. (17) with right-hand side computed from a realistic EoS, this initial approximation had to be very close to the final solution to ensure a convergence of the iteration procedure. Due to these problems, this integration procedure was considered as impractical for routine calculations.

To overcome the above-mentioned difficulties, a new algorithm was developed based on the simple but robust shooting method. The boundary value problem was converted into an initial value problem by estimating the density in the center of the droplet $\rho_{0}$, so the initial conditions were

[TABLE]

The initial value problem was then solved using an implementation of the variable-step Runge–Kutta method (matlab function ode45). The original boundary condition in Eq. (18), in theoretical formulation placed to infinity, can be modified to $\rho(R)=\rho_{\mathrm{V}}$ where $R$ is finite but sufficiently large, such that the theoretical density profile at $R$ is sufficiently close to the vapor density. In the shooting method, the center density $\rho(0)$ is searched for in an iteration procedure such that the right-hand boundary condition is matched.

As shown in Fig. 2, numerical solutions of Eq. (17) with conditions defined by Eq. (31) indicate an oscillatory behavior. We note that these oscillations are not a result of a particular numerical procedure but rather they are an intrinsic feature of Euler-Lagrage equation, Eq. (17). In this example, a decrease from density in the center $\rho_{0}=713\,\mathrm{kg\,m^{-3}}$ down to the vapor density $\rho_{\mathrm{V}}=0.36\,\mathrm{kg\,m^{-3}}$ occurs within the first 3 nm. Then the curve detaches from the vapor density and oscillates around density $\rho\simeq 250\,\mathrm{kg\,m^{-3}}$, corresponding to a physically unstable (but numerically stable) solution of the equality of chemical potentials, Eq. (2). This oscillatory behavior was observed for estimates of $\rho_{0}$ lying between the correct value and the saturated liquid density $\rho_{\mathrm{L}}$. Solutions of Eq. (17) for guesses of $\rho_{0}$ slightly below the correct value approached the vapor density and then detached downwards, quickly crossing the zero density. Because of this behavior, enforcing the boundary condition at certain fixed coordinate $R$ was found impractical. Therefore, the second condition of Eq. (31) was changed to

[TABLE]

In addition, estimates undershooting the vapor density have been rejected. As a result, the proper value $\rho_{0}$ was approached from below (since lower $\rho_{0}$ corresponded to higher $\rho_{\mathrm{V}}$). To achieve a fast convergence, we attempted to compute a next approximation of $\rho_{0}$ based on the linear extrapolation of the previous two guesses and the corresponding difference of the left and right sides of Eq. (32). Unfortunately, this modified secant method appeared to be unreliable because the numerical solution of Eq. (17) exhibits a quasi-random error when considered as a function of parameter $\rho_{0}$. The magnitude of this numerical error is small enough so that it does not affect the computations of the work of formation. However, the non-smooth dependence on $\rho_{0}$ caused convergence failures for $\rho_{0}$ guesses close to the correct boundary value solution. This problem was remedied by changing the secant method to a slower but more robust bisection method. Additional improvementsPlanková (2011) were developed to make the algorithm more robust and precise. In particular, the inner and outer parts of the density profile in intervals, where the compressibility only negligibly differs from compressibility at $\rho_{0}$ or at $\rho_{V}$, were substituted by an analytical solution of a linearized Eq. (17).

VII Results and discussion

We performed DGT and CNT computations for four alkanes; namely $n$-heptane, $n$-octane, $n$-nonane, and $n$-decane, for temperatures corresponding to nucleation rate experimental data by Rudek et al.,Rudek et al. (1996) Hung et al.,Hung, Krasnopoler, and Katz (1989) Luijten,Luijten (1998) Viisaanen et al.Viisanen, Wagner, and Strey (1998), and Wagner and Strey,Wagner and Strey (1984)

• •

$n$-heptane (C7): 249 K, 259 K, 268 K, 276 K,

• •

$n$-octane (C8): 241 K, 248 K, 258 K, 267 K, 287 K, 298 K, 302 K,

• •

$n$-nonane (C9): 220 K, 230 K, 240 K, 247.6 K, 257 K, 267 K, 272 K, 284 K, 298 K, 313 K,

• •

$n$-decane (C10): 268 K, 276 K, 281 K, 294 K, 303 K, 316 K, 330 K.

A comprehensive experimental study of nucleation in $n$-alkanes from $n$-heptane to $n$-decane in nitrogen or argon has recently been provided by Ghosh et al.Ghosh et al. (2010) We did not include it in the present comparisons because the very high nucleation rates (of the order of $10^{22}$ to $10^{23}$ $\mathrm{m^{-3}s^{-1}}$) make the comparison with the other data difficult. However, we note that their comparison with CNT prediction leads to similar conclusions, particularly concerning the strong temperature dependence of the deviation (cf. their Fig. 9).

Computations were done using PC-SAFT EoS and PR EoS. Influence parameters were computed using Eq. (20) based on the “bare” surface tension of the planar phase interface that was obtained from the experimental surface tension using Eq. (28). (See Appendix B for more information about the experimental surface tensions.) The computations of the density profiles for critical clusters were performed using the algorithm described in Sec. VI. The nucleation rates were evaluated using Eq. (9), DGT works of formation were obtained using Eq. (15) for the computed density profiles, the works of formation for CNT were found using Eq. (10) with $\sigma=\sigma_{\infty}$.

Figs. 3 to 6 show nucleation rates depending on the supersaturations, Eq. (1), computed using the various combinations of DGT and CNT with PC-SAFT and PR EoSs compared with the experimental data in all figures. One color and one symbol always correspond to one temperature. Lines depict the theoretical data: solid lines with light-colored symbols are the DGT and PC-SAFT computations, dashed lines with light-colored symbols correspond to CNT and PC-SAFT, solid lines with white symbols to DGT and PR, dashed lines with white symbols to CNT and PR.

The darkest symbols depict the experimental data. The experiments were not conducted at the same nominal temperatures $T_{0}$ which correspond to the temperatures in the figures as the experimental temperatures $T_{\mathrm{exp}}$ slightly varied. To enable a direct comparison, the experimental supersaturations $S_{\mathrm{exp}}$ were corrected to match the nominal temperature. A Taylor expansion was done for $\ln S_{\mathrm{exp}}$ around the experimental temperature $T_{\mathrm{exp}}$ to match the desired temperature $T_{0}$,

[TABLE]

where we defined $S_{0}=S(T_{0},J_{\mathrm{exp}})$. However, the derivative ${\partial\ln S}/{\partial T}$ would be difficult to obtain, hence it was substituted by

[TABLE]

These derivatives were evaluated by differentiating Eq. (9).

The slopes of experimental data with respect to supersaturation relatively well correspond to the theoretical slopes. However, all the lines show a systematic temperature-dependent shift.

Another way how to study the quality of theoretical predictions is via ratios of experimental to predicted nucleation rates. Figs. 7 and 8 show ratios of theoretical and experimental nucleation rates corresponding to the experimental temperatures and supersaturations as functions of inversed reduced temperatures. Numerical data were twice interpolated using cubic splines; first, $\ln J$ was interpolated with respect to $\ln S$ on two theoretical isotherms neighboring $T_{\mathrm{exp}}$ such that their $\ln S$ corresponded to the experimental value $\ln S_{\mathrm{exp}}$. Then these two nucleation rates were again interpolated with respect to temperature to match $T_{\mathrm{exp}}$.

Fig. 7 gives the ratios of the experimental nucleation rates to the nucleation rates computed using DGT and PC-SAFT EoS. In the figure, one color and one symbol correspond to one experimental data set. Results of DGT with PC-SAFT EoS and CNT with PC-SAFT EoS are compared in Fig. 8. The DGT results show slightly smaller temperature trend of the deviation from experimental data.

Clearly, the experimental data show a strong temperature trend which is not predicted by any of the theories. The ratios of experimental to theoretical nucleation rates for the various alkanes appear to follow almost a single linear function of the reciprocal reduced temperature,

[TABLE]

where $A$ and $B$ are empirical parameters. Empirical correction, Eq. (35), is of the same form as proposed by Wölk et al.Wölk et al. (2002) for correcting the CNT predictions of nucleation rates for water. As discussed by McGraw and LaaksonenMcGraw and Laaksonen (1996), the dependence of the ratio of experimental to theoretical nucleation rates on supersaturation is rather weak; in other words, the dependence of the nucleation rate on the supersaturation is predicted relatively well by the theory. However, some discrepancies can be observed for the data of $n$-nonane, for which the experiments were done in a wide range of supersaturations. A failure of CNT to properly describe the supersaturation dependence accurately has been discussed by Girshick.Girshick (2014) We propose an additional term describing the dependence on the supersaturation,

[TABLE]

We used the least-squares method to determine the parameters in Eq. (36) for $n$-nonane for both nucleation theories and both EoSs. Except for $n$-nonane, the nucleation rate data for selected $n$-alkanes are available in a rather narrow supersaturation range. This makes it impossible to obtain a reliable value for parameter $C$. Therefore, we used the numerical value of coefficient $C$ for $n$-nonane as an approximation of the supersaturation dependence for other alkanes. Figs. 9 and 10 show the ratios of experimental nucleation rates $J_{\mathrm{exp}}$ to predictions of Eq. (36). These figures show that the correlation successfully captures the variability of the experimental data. Values of coefficients $A,B$, and $C$ are listed in Tab. 1 for DGT and CNT combined with PC-SAFT and PR EoSs.

Typically, the experiments provide nucleation rates in a relatively narrow range determined by the experimental method. Therefore, it is practical to represent experimental data as a temperature dependence of supersaturation $S_{\mathrm{ref}}$ corresponding to a chosen reference nucleation rate $J_{\mathrm{ref}}$. Figs. 11 and 12 show a temperature dependence of the supersaturation depending on the inverse reduced temperature $T_{\mathrm{c}}/T$ for the selected four substances. Numerical calculations were performed in a wider range of temperatures compared to experimental data in order to provide a better picture of the temperature trend. Experimental data by Rudek et al.Rudek et al. (1996) and Hung et al.Hung, Krasnopoler, and Katz (1989) were measured in a substantially lower range of nucleation rates than by Luijten,Luijten (1998) Viisaanen et al.,Viisanen, Wagner, and Strey (1998) and Wagner and Strey Wagner and Strey (1984). For this reason, we used two reference nucleation rates $J_{\mathrm{ref}}=7.57\times 10^{4}\mathrm{m}^{-3}\mathrm{s}^{-1}$ and $J_{\mathrm{ref}}=3.22\times 10^{13}\mathrm{m}^{-3}\mathrm{s}^{-1}$, which were obtained as geometrical averages of the nucleation rates included in each data subset. We note that Figs. 11 and 12 have have the same ranges of both $x$ and $y$ axes.

In order to reduce the experimental data obtained for various nucleation rates to a given reference nucleation rate $J_{\mathrm{ref}}$, the theoretical data were interpolated using a cubic spline. For experimental data, we assumed a linear dependence,

[TABLE]

where coefficients $b_{1}$ and $b_{2}$ were fitted using the least squares method for each isothermal experimental data set. One symbol in Figs. 11 and 12 (and later also in Figs. 13 and 14) corresponds to one isothermal data set of the measurements, typically including five measurements. Supersaturations $S$ range from 5.39 to 43.18 for the case of $J_{\mathrm{ref}}=7.57\times 10^{4}\mathrm{m}^{-3}\mathrm{s}^{-1}$ (Fig. 11) and from 74.94 to 181.525 for the case with $J_{\mathrm{ref}}=3.22\times 10^{13}\mathrm{m}^{-3}\mathrm{s}^{-1}$ (Fig. 12).

Figs. 13 and 14 show a temperature dependence of the number of molecules in the critical cluster $N^{*}$ depending on the inverse reduced temperature $T_{\mathrm{c}}/T$. The data was fitted to match one reference nucleation rate as discussed above. For the DGT predictions, the number of molecules $N^{*}$ was evaluated as

[TABLE]

where $\Delta N_{\mathrm{N^{*}}}$ is the excess number of molecules given by Eq. (19) and $V_{\mathrm{N^{*}}}$ is the volume of the droplet consisting of $N^{*}$ molecules. A rational estimate of volume is the volume of a sphere enclosed by the surface of tension, $V_{\mathrm{N^{*}}}=4/3\pi r_{\mathrm{s,N^{*}}}^{3}$.

Regarding the experimental data, the number of molecules of the critical cluster was evaluated using the nucleation theorem,Oxtoby and Kashchiev (1994)

[TABLE]

As the derivative can be fitted by Eq. (37), $N^{*}=b_{2}-1$.

Again, we split the experimental data into two groups and used reference nucleation rates $J_{\mathrm{ref}}=7.57\times 10^{4}\mathrm{m}^{-3}\mathrm{s}^{-1}$ and $J_{\mathrm{ref}}=3.22\times 10^{13}\mathrm{m}^{-3}\mathrm{s}^{-1}$. The number of molecules in critical cluster $N^{*}$ in Figs. 13 and 14 range from 34.2 to 68.9 for the case of $J_{\mathrm{ref}}=7.57\times 10^{4}\mathrm{m}^{-3}\mathrm{s}^{-1}$ and from 10.1 to 16.4 for the case with $J_{\mathrm{ref}}=3.22\times 10^{13}\mathrm{m}^{-3}\mathrm{s}^{-1}$.

As can be seen in Figs. 11, 12, 13, and 14, the theories provide qualitatively correct predictions, despite the large disagreement in terms of the nucleation rate. However, a noticeable disagreement can be seen in the slopes. Predicted derivatives of the supersaturation $S$ and of the number of molecules in the critical cluster $N^{*}$ with respect to temperature at a constant nucleation rate differ from the experimental slopes, especially at low temperatures.

As discussed in Sec. V, CW broaden the phase interface and somewhat decrease the experimental surface tension with respect to the bare surface tension corresponding to a phase interface constrained to be free of CW. Fig. 15 shows the surface tension depending on the Laplace pressure $\Delta p$ computed using DGT combined with the two EoSs (PR and PC-SAFT), and with two different definitions of the influence parameters. The figure demonstrates an example for $n$-heptane at temperature of $276$ K. The Laplace pressure is connected with the surface-of-tension radius of the droplet $r_{s}$ via the Young–Laplace equation, Eq. (11). Therefore, $\Delta p=0$ Pa at left corresponds to an infinite radius, i.e., the planar phase interface with surface tension $\sigma_{\infty}$. The size-dependent surface tension was obtained from Eq. (12) based on the work of formation, Eq. (15). At the right end, the surface tension according to DGT vanishes at the the spinodal. Lines $c_{\mathrm{bare}}$ connect results for influence parameters computed using Eq. (20), where $\sigma_{\infty}$ was cleared of CW using Eq. (28). This is a prediction considered as correct, because, as discussed in Sec. V, the surface of the critical clusters is practically free of CW and, at the same time, the bare surface tension is the appropriate value to compute the influence parameter as DGT does not include CW. Lines $c_{\mathrm{exp}}$ correspond to the influence parameters computed via Eq. (20) using the experimental surface tension. This is the standard approach; however, it is considered as less accurate. From Fig. 15 it can be seen that $\sigma_{\mathrm{bare}}>\sigma_{\mathrm{exp}}$. The surface tension curves for both EoSs start at either $\sigma_{\mathrm{bare}}$ or at $\sigma_{\mathrm{exp}}$ at $\Delta p=0$ Pa, but then tends to a spinodal value given by the particular EoS. The inset figure shows a detail of the trends for low Laplace pressures. The range of Laplace pressures for the critical clusters in the experiments is on the descending branches of the curves close to the maximum. In the case shown in Fig. 15, the Laplace pressures range from 27.3 to 29.9 MPa, corresponding to critical radii of 1.64 to 1.50 nm. It can be seen that CW have a significant effect on the surface tension. This effect is bigger or at least comparable (in other cases) with the effect of replacing the cubic PR EoS by the more accurate PC-SAFT EoS. A reason can be seen in a crossing of the lines predicted by both EoSs just in the relevant range.

Fig. 16 shows nucleation rates depending on the supersaturation of $n$-heptane at temperature of $276$ K. Again, the two EoSs and two different values for the influence parameters were used, similarly as in Fig. 15. The stars depict the experimental data. Constraining the phase interface to be free of CW increases the surface tension and, consequently, enhances the work of formation and lowers the nucleation rates. The effect is quite large, which proves the importance of such a treatment. We note that the quantitative agreement of the experimental data with the theoretical prediction shown in Fig. 16 is rather fortuitous: due to the strong temperature effect, the differences are larger at other temperatures.

VIII Conclusions

We computed nucleation rates of four alkanes: $n$-heptane, $n$-octane, $n$-nonane, and $n$-decane, using the density gradient theory complemented by an original treatment of capillary waves, and the classical nucleation theory. For both approaches we used a simple cubic Peng–Robinson EoS, which essentially is an empirical modification of the van der Waals EoS, and the PC-SAFT EoS, which is a modern EoS based on molecular arguments. The computed predictions were compared to the experimental nucleation rate data.

An original algorithm was developed to solve an Euler–Lagrange equation of DGT which together with the boundary conditions forms a boundary value problem. Several enhancements of the numerical method were developed in order to ensure a robust convergence to the physically-proper density profiles.

For computing the nucleation rates, we used an expression based on the classical nucleation kinetics with corrected supersaturation dependence, Eq. (9). The work of formation was obtained using DGT combined with our CW treatment or, in case of CNT, using the Gibbs formula based on the surface tension for the planar phase interface.

In our investigation of the effect of capillary waves on the surface tension and nucleation rates we found that the critical clusters are so small that the shortest wavelengths of CW are already too large to fit to the surface of the nanoscopic droplet representing the critical cluster. This means that the thermal motion of the molecules and the corresponding fluctuations of density and other local quantities are completely accounted for in the DGT description. Therefore, for a consistent computation of the works of formation using either DGT or CNT, it is appropriate to remove the effect of CW from the planar surface tension used in the computations. This was achieved using Eq. (21) based on the theory of Meunier.Meunier (1987) In DGT, the influence parameter was fitted to reproduce the bare surface tension corresponding to the planar phase interface free of CW.

Three contributions to the surface tension of a droplet can be recognized in Fig. 15. First is the Tolman’s linear effect which accounts for the mild increase of surface tension (for low Laplace pressure $\Delta p$). This effect is related to symmetry properties of the density profile, expressed by the Tolman length. Second, quadratic effect, appears to be decisive in bending the curvature such that the surface tension reaches a maximum and then decreases towards the spinodal. The quadratic effect can be attributed to the finite thickness of the phase interface and to changes in its shape as the critical cluster becomes smaller. Third is the CW effect which decreases the surface tension of nano-droplets as a part of the spectrum of the CW that do not fit to the finite dimensions of the droplet. Suppression of CW leads to a considerable enhancement of the predicted surface tensions and, in turn, to a decrease in the nucleation rates.

A pre-requisite for realistic modeling using DGT is an accurate EoS. The PC-SAFT EoS provides a more accurate representation of thermodynamic properties in the stable liquid and gas regions covered with experimental data. In the Supplementary Material 2, we provide some comparisons showing that PC-SAFT EoS is clearly superior to the PR EoS, which is generally considered as a reasonable standard for modeling hydrocarbon systems. Besides, the molecular foundations of PC-SAFT EoS allow to consider the predictions for the metastable and unstable regions as plausible. Fig. 2 in the Supplementary Material 2 shows that the van der Waals loops predicted by both EoSs differ substantially. Despite that, the difference in the predicted surface tensions of critical clusters and, consequently, in the predicted nucleation rates, is not tremendous, at least for the systems studied here. The reason is shown in the inset of Fig. 16: Incidentally, the curves predicted by DGT combined with PC-SAFT EoS intersect in the region corresponding to the critical clusters.

This study has been performed with the goal of possibly finding a reason for a disagreement of CNT with experimental data, which is particularly striking in the case of $n$-alkanes. In Figs. 3 to 6, we provide comparison of theoretical predictions with experimental nucleation rate data. The differences between the DGT/CW and CNT predictions is significant, reaching several orders of magnitude in the nucleation rate. Also the difference between predictions by the different EoSs is appreciable. Despite that, the huge deviation from the experimental data, particularly in the temperature dependence, was not remedied. To characterize the deviation from experiments, various approaches can be adopted. One can deduce the magnitude of the microscopic surface tension (corresponding to the critical cluster) from nucleation rate data. For example, Baidakov and SkripovBaidakov and Skripov (1982) found that for bubble nucleation in several liquified gases it is appropriate to reduce the macroscopic surface tension by a few percent. In the present case, this approach does not appear to be insightful. Due to the strong temperature dependence, the correction to the surface tension would change from large positive values at low temperature to large negative corrections at the high-temperature limit of the experimental window. In section II we noted that there are several modifications to CNT differing in the pre-exponential factor. The quantity $\vartheta$, the “translational volume scale”,Reiss, Kegel, and Katz (1998) remains, despite numerous theoretical contributions,Reguera and Reiss (2004) a rather uncertain quantity. For this reason, we consider as more appropriate to represent the deviation in terms of a multiplicative correction to the predicted nucleation rate. We provide this correction in form of Eq. (36) for all combinations of theories (DGT/CW, CNT) and equations of state (PC-SAFT, PR) to facilitate further theoretical analyses as well as practical computations. The dominant term is the temperature correction. We also included a correction to the supersaturation dependency which is weaker and difficult to asses, because it relies on comparison of various experimental methods, which might be subject to systematic offsets.

We presented partial improvements to the classical nucleation theory which bring the modeling closer to reality. However, the resolution of the temperature-dependent deviation,Malila et al. (2008) which is very strong in case of $n$-alkanes, remains unresolved.

Supplementary material

The data points of the figures presented in this article are listed in the Supplementary Material 1. Additional figures are included in Supplementary Material 2 that compare saturated properties of the four $n$-alkanes for both EoSs and justify superiority of PC-SAFT EoS over PR EoS. Further figures of nucleation rates are also included.

Acknowledgements.

The research leading to these results has received funding from the Norwegian Financial Mechanism 2009-2014 under Project Contract No. 7F14466, from the Czech Science Foundation grant No. GJ15-07129Y, and from the institutional support RVO: 61388998. Besides, B.P. thanks for a partial support from the grant No. SGS15/215/OHK4/3T/14.

Appendix A Parameters of equations of state

For PC-SAFT EoS, molecular parameters $m,\sigma$, and $\epsilon$ for selected $n$-alkanes used in this study were taken from the original work of Gross and SadowskiGross and Sadowski (2001).

The Peng–Robinson equation is given by

[TABLE]

where $R_{\mathrm{g}}=8.3144598\,\mathrm{J\,mol^{-1}\,K^{-1}}$ is the universal gas constant and $a,b,\alpha$ are Peng–Robinson constants given as follows

[TABLE]

where $p_{\mathrm{c}}$ is the critical pressure and $\omega$ is the acentric factor. Values fot he acentric factor were taken from the work by Lemmon and Goodwin.Lemmon and Goodwin (2000) Critical parameters $T_{\mathrm{c}}$ and $p_{\mathrm{c}}$ and acentric factor $\omega$ for four selected $n$-alkanes are listed in Tab. 2.

Appendix B Experimental surface tension

The influence parameter $c$ was computed using the experimental surface tension. The experimental values were taken from Jasper et al.Jasper, Kerr, and Gregorich (1953) for all the selected substances, GrigoryevGrigoryev et al. (1992) for $n$-heptane (C7) and $n$-octane (C8), Voliak and AndreevaVoliak and Andreeva (1961) for C7 and C8, Rolo et al.Rolo et al. (2002) for C7 and $n$-decane (C10), Okada et al.Okada et al. (1990) for C8.

The actual value of the surface tension was computed using a temperature correlation. For $n$-heptane and $n$-octane, we used approach by SomayajuluSomayajulu (1988),

[TABLE]

For $n$-nonane and $n$-decane we developed our own fit based on the scaling around the critical temperature,

[TABLE]

Values of parameters $A$, $B$, and $C$ from Eqs. (42) and (43) are listed in Tab. 3.