Homogeneous droplet nucleation modeled using the gradient theory combined with the PC-SAFT equation of state

In this work, we used the density gradient theory (DGT) combined with the cubic equation of state (EoS) by Peng and Robinson (PR) and the perturbed chain (PC) modification of the SAFT EoS developed by Gross and Sadowski [1]. The PR EoS is based on very simplified physical foundations, it has significant limitations in the accuracy of the predicted thermodynamic properties. On the other hand, the PC-SAFT EoS combines different intermolecular forces, e.g., hydrogen bonding, covalent bonding, Coulombic forces which makes it more accurate in predicting of the physical variables. We continued in our previous works [2,3] by solving the boundary value problem which arose by mathematical solution of the DGT formulation and including the boundary conditions. Achieving the numerical solution was rather tricky; this study describes some of the crucial developments that helped us to overcome the partial problems. The most troublesome were computations for low temperatures where we achieved great improvements compared to [3]. We applied the GT for the n-alkanes: nheptane, n-octane, n-nonane, and n-decane because of the availability of the experimental data. Comparing them with our numerical results, we observed great differences between the theories; the best results gave the combination of the GT and the PC-SAFT. However, a certain temperature drift was observed that is not satisfactorily explained by the present theories.


Introduction
Homogeneous nucleation of droplets plays an important role in processes in the atmosphere, since it describes forming of the secondary aerosols, and in many technological processes such as nucleation in the crude nature gas during the gas cleaning, nucleation in steam turbine, etc.
Despite many attempts that resolved partial subproblems of the nucleation, there is no complete theory which would give quantitatively correct predictions.
Studying the nucleation has its roots in the work of Gibbs and Volmer and Farkas; however, it was mainly developed by Bercker end Döring.Their theory is nowadays called the classical nucleation thery (CNT).Its main idea is based on the capillarity approximation of the phase interface of the droplet which assumes a sharp phase interface.However, in reality the thickness of the interface is often comparable with the size of the droplet and the sharp interface assumption becomes a source of errors.Therefore, it is more realistic to consider a smoother change of the interface.
This approach was developed by Cahn and Hilliard [4,5] although the main ideas were given already by van der Waals [6] (English translation by Rowlinson [7]).Main idea of their theory (now known as DGT) was to extend the Helmholtz free energy of the system containing the droplet by the term of gradient of the particle density.Theoretical foundations of the gradient theory and the related statistical physics background are given by Rowlinson and Widom [8], Davis [9], and Kalikmanov [10].A recent study of the gradient theory with comparison to molecular simulations was given by Baidakov et al.
Correspondence to: e-mail: barbora.plankova@gmail.com In this work, we continue in our approach [3] to concentrate on the nucleation of one component using the DGT, which mathematically leads to the problem how to solve the boundary value problem for the ordinary differential equation of the first order.Although this problem seemed simple, many technical problems arose during solution; in this paper, they will be briefly mentioned.For more information regarding this issue, see [11].
The DGT combined with the suitable equation of state (EoS) provides rates at which new droplets condensate, nucleation rates, the main property of the nucleation phenomena.In this work, we compare the widely used cubical Peng-Robinson (PR) EoS and the modern PC-SAFT EoS with experiments.
Cubic PR EoS provides a qualitatively good picture of the thermodynamics, however, it is quantitatively inaccurate compared to the PC-SAFT.
In this work we compare both EoSs, PC-SAFT EoS with PR EoS, and the DGT with the CNT with experiments.We performed DGT numerical computations for nalkanes: n-heptane, n-octane, n-nonane, and n-decane because of the presence of the experimental data.
We observed that the combination of PC-SAFT EoS with DGT is closest to the experimental data.However, we encountered a systematic temperature drift of both nucleation theories to the experimental data which is not explained by present theories.pressure p and chemical potential µ of the liquid and vapor are same.This state will be denoted by lower index ∞ .If the density of the gas is greater than the saturated value, ρ G > ρ G∞ , this thermodynamic state is then called supersaturated.Quantity describing this supersaturation is defined as where k B is the Boltzmann constant.Supersaturated system is in metastable state: it is stable to small fluctuations but unstable to large ones.Droplet nucleation is the process of forming liquid droplets in the supersaturated vapor phase.
Microscopically, a new cluster is created when at least two monomers collide.This cluster can vanish by spontaneous abandoning one monomer, or grow by joining another monomer.By this stochastic process a relative large cluster can be formed.When its size exceeds the so-called critical size, depending on the supersaturation, the probability of addition of a molecule exceeds the probability of evaporation and the cluster can grow to a macroscopic droplet.The number of droplets formed in unit of volume per unit of time is called nucleation rate, J.Besides the supersaturation of the vapor phase, also temperature T is another parameter of nucleation rate.
In classical way, nucleation rate can be expressed as [12] where ρ G is the density of the gaseous state, ρ L∞ is the density of liquid of the saturated state, N A is the Avogadro constant, σ ∞ is surface tension of the saturated state, M is molecular mass, k B is the Boltzmann constant, and ∆Ω is the work of formation of the critical cluster.According to the CNT, the work of formation of the critical cluster reads where A s = 4πr 2 s is a surface of the droplet.Work of formation according to the GT will be given in next section.
Although Eq. ( 2) is convenient as an expression of the nucleation rate, it has one flaw: inconsistency to the monomers.Internally consistent (IC) expression was proposed by Girshick and Chiu [13], where is the dimensionless surfaces tension.Using the definition of the supersaturation S Eq. ( 1) expression (3), and the Young-Laplace equation the dimensionless supersaturation and the work of formation can be linked as Let us emphasize that last equality is valid only for the CNT.

Density gradient theory
The DGT uses the density of the molecules ρ = N/V instead of their exact number.When working with the spherical droplets, it is natural to use spherical coordinates.Since the spherical symmetry is assumed, the only variable of the density is the distance from the center of the droplet, r.The Helmholtz free energy density φ is according to the DGT given by First term is the Helmholtz free energy density of the homogeneous fluid at actual density ρ.Second term corresponds to the inhomogeneity given by the interface.Parameter c in this term is called the influence parameter.
It is generally temperature and density-dependent and is a property of the given substance.In this work (as done usually), we will consider c as density-independent.For the so-called critical cluster, the chemical potential is homogeneous throughout the system; therefore, the grand potential density ω is given by where µ G is the chemical potential of the gaseous state.
The work of formation of the cluster, needed for evaluation of nucleation rates, can be expressed as a difference between the grand potential of the system including the cluster and the grand potential of the homogeneous system not constrained to contain the cluster, where ∆ω hom is given by ∆ω hom can be obtained using the suitable EoS.However, the work of formation depends on the density profile of the droplet ρ(r) which has to be, therefore, computed first.In a very similar way it is possible to obtain the surface tension of the saturated state, where z is the coordinate perpendicular to the interface.This surface tension can be measured in an experiment, therefore, we used Eq. ( 12) for computation of the influence parameter c.Applying variation methods onto Eq. ( 10), a differential equation for the density profile of critical droplet can be derived.Critical droplet has such a size among all the possible sizes of the droplet, that has the same probability to shrink as to grow.Density profile of the critical droplet is a non-trivial solution of the Euler-Lagrange equation where ∆µ = µ(ρ) − µ G .This equation must be solved for boundary conditions

Numerical solution
Equation ( 13) with boundary conditions ( 14) forms a boundary value problem (BVP).This simply looking problem has two big difficulties: the density profile near the gaseous phase has a very sharp shape; its slope changes abruptly from the very steep decline to an almost constant profile.Second problem is that for large droplets the density profile in the interior of the droplet changes only negligibly and is almost constant.This causes very significant cumulation of numerical errors.In previous work [2], this BVP was solved using the matlab procedure BPV4 [14].This procedure uses the finite difference scheme: it devides the original interval of the problem into subintervals, where the right-hand side is integrated.This procedure requires the initial net of intervals and initial approximation of the solution.However, this initial approximation had to be very close due to the problems mentioned above.
This work is based on results of [11], in which the shooting method instead of the BVP4 procedure was used.To overcome many difficulties that arose during the solution process, several original numerical methods were developed.All the procedure and the algorithm is described in detail in [11].

Radius of surface of tension and equimolar radius
Radius of surface of tenstion r s is a such radius for which the Young-Laplace equation has a form of Eq. ( 6).On the other hand, equimolar radius r e relates to the number of molecules.Let us denote N G nad N L as number of gas molecules and number of liquid molecules, respectively.The rest of the molecules then would be N S and corresponds to the number of molecules absorbed by the interface.These three numbers depend on the choice of the interface position.Equimolar radius r e then correspond to the choice such that N S = 0.
Excess number of molecules ∆N is the difference between the number of molecules of the actual droplet and number of molecules of the bulk liquid of same volume.It is the property of the interface.Its value is given by

PC-SAFT equation of state
This physically based EoS was introduced by Chapman et al. [15,16].It is based on the Statistical Associating Fluid Theory (SAFT) combining important intermolecular forces, such as hydrogen bonding, covalent bonding, Coulombic forces and can be used for very different shapes of molecules.
Due to the fact that the SAFT EoS works directly with the molecular structure of substances, it allows modeling fluids in the metastable region, which is our case.
The PC-SAFT model extends the SAFT model using the three additional substance parameters: the number of segments per chain m, the segment diameter s, and the segment energy parameter .
Molar mass, critical temperature and pressure along with the PC-SAFT parameters for four substances considered in this work (n-heptane, n-octane, n-nonane, n-decane) are summarized in table 1.
The influence parameter c was computed using Eq. ( 12) of the surface tension of the saturated state.We computed c for each temperature based on experimental values of the surface tension from Jasper and et al. [22] for all the four substances, Grigoryev [23] for C7 and C8, Voliak and Andreeva [24] for C7 and C8, Rolo et al. [25] for C7 and C10, Okada et al. [26] for C8.

Nucleation rates
For the nucleation rate we used combination of classical (2) and internally consistent (4) approach, for the GT values and, also using the expression (7), for CNT values.Reason for neglecting the term θ in the exponent was the would-be theory-inconsistency in (16) in which the combination of GT and CNT would be used.Justification of this negligence lies in the fact that θ in the exponent has no physical meaning; it is merely empirical term.
Figures 1 -6 show nucleation rates computed using the improved CNT and GT by the internally consistent term 01076-p.3expressed in Eqs. ( 16) and ( 17) for all four substances nheptane, n-octane, n-nonane, and n-decane and for temperatures stated earlier in this section.Using the PR and PC-SAFT EoSs, nucleation rates were computed for various supersaturations to cover the range of experimental data.Unlike the numerical data that have been computed for the temperatures stated above, temperature of the experimental data was slightly different and had to be interpolated.
For this purpose, one-dimensional Taylor expansion up to the linear term of the logarithmical value of supersaturation ln S around the experimental temperature T exp in the point of computed temperature T 0 was used.Logarithm of the experimental value of the nucleation rate ln J exp remained constant as a parameter.Therefore, interpolated value of the new ln S 0 ≡ ln S (T 0 , ln J exp ) reads It was more accurate to use the logarithms of supersaturation S and nucleation rate J because of the exponential character of their dependence.We did not have any explicit expression of ln S as a function of T , but we could easily get expression ln J as a function of temperature T and logarithm of supersaturation ln S , ln J = ln J(T, ln S ) and due to the mathematical rule of changing the derivative variable, we rewrote the derivative in (18) into Using CNT-IC expression for nucleation rate (17) and supersaturation approximation S =ρ G /ρ G∞ , derivatives on the right-hand side of ( 19) can be expressed as According the definition of θ (5), its derivative in (20) reads Remaining derivatives were in this work computed numerically, using differences.Marker colours for nucleation theories and line types for EoSs are same as in figure 1, experimental data are given by Wagner and Strey [21] (circles) for 220 K and (left-oriented triangles) for 230 K, by Hung et al. [18] (diamonds), Luijten [19] (right-oriented triangles), and Viisaanen et al. [20] (squares) for 230 K.
Figs. 1 -6 show that for all of the components and most of their temperatures the combination of PC-SAFT and GT-IC best matched the experimental values.Main improvement brought the use of the PC-SAFT; difference between GT-IC and CNT-IC was not so pronounced.However, this was caused by the range of the experimental data; for higher supersaturations, their values are more distinct.Still, this work made great improvement to the nucleation rates prediction.

Nucleation rates ratio
Another way how to compare numerical and experimental data gives the ratio of their nucleation rates.Figure 7 shows the ratio of nucleation rates computed using the GT-IC, J GT−IC , and the experimental data J exp as a function of the inverse reduced temperature T c /T .T c is different for each substance and are listed in table 1.
Numerical value J GT−IC is interpolated in two steps to compute nucleation rates corresponding to the experimental temperature T exp and experimental supersaturation S exp for each experimental value.
First, nucleation rate is interpolated with respect to supersaturation.Similarly to Sec. 5.1, logarithms ln S , ln J are used instead of original values.The procedure is as follows: first, for given experimental temperature T exp the lowest higher and the highest lower numerical isotherms T low , T high are found.Then, nucleation rates are interpolated with respect to the supersaturation to match the value ln S exp on both isotherms with outputs ln J low , ln J high .
Second, nucleation rate is interpolated with respect to temperature between the two values ln J low and ln J high to match T exp .  .Nucleation rates J of n-nonane as a function of supersaturation S computed for 240 K (circles), 247.6 K (diamonds), 257 K (triangles), and 267 K (squares).Marker colours for nucleation theories and line types for EoSs are same as in figure 1, experimental data are given by Hung et al. [18] for 240 K (circles), 247.6 K (diamonds), 257 K (hexagrams), and 267 K (squares), by Wagner and Strey [21] (pentagrams) for 240 K, and by Rudek et al. [17] for 257 K (triangles).
Figure 8 shows the ratio of nucleation rates computed using the GT-IC J GT−IC and using the CNT-IC J CNT−IC .Both numerical values were interpolated using the procedure described above to match the experimental values T exp , S exp .
Range of the nucleation rate ratios in figure 7 is quite large which points to not very good temperature dependence of the GT-IC.
Values of nucleation rate ratios in figure 8 shows that numerical values computed using the GT-IC are rather higher than CNT-IC values.Points corresponding to Wagner's and Viisanen's measurements are out of the bulk cloud.This is caused by the fact, that these measurements were done for higher supersaturations than the others, where CNT and GT differ significantly.

Temperature dependence and number of molecules
Figures 9 and 10 show the temperature dependence: dependence of the supersaturation S on the inverse reduced temperature T c /T of all the substances C7-C10.GT-IC and CNT-IC are computed using the PC-SAFT EoS.Parameter for all the values is the reference nucleation rate J ref = 10 12 m −3 s −1 for figure 9 and J ref = 10 6 m −3 s −1 for figure 10.These two values were chosen as a reasonable values corresponding to the experimental data: the former for [19][20][21], the latter for the rest.
Numerical values of ln S are interpolated for for each temperature T with respect to ln J to match the value ln J ref .
From Figs. 1 -6, linear dependence of ln J exp on ln S exp was assumed (both axis on all the figures are in logarithmic Least squares method was used to find the coefficients b   .Ratio of nucleation rates J GT−IC evaluated using the expression (16) and experimental values as a function of inverse reduces temperature T/T c , where critical temperatures T c are listed in table 1. Experimental data are given by Rudek et al. [17] for nheptane (circles), n-octane (squares), n-nonane (pentagrams), and n-decane (diamonds), by Hung et al. [18] for n-nonane (crosses), by Luijten [19] for n-nonane (stars), by Viisaanen et al. [20] for n-nonane (triangles), and by Wagner and Strey [21] for n-nonane (rotated crosses).  .Ratio of nucleation rates J GT−IC evaluated using the expression (16) and nucleation rates J CNT−IC by Eq. ( 17).Temperatures and supersaturations correspond to the listed experimental data.Markers are same as in figure 7.
tor f such as 1 ln S exp (1) 1 ln S exp (2) . . . . . .   1) for nucleation rate J = 10 12 m −3 s −1 .Markers of the experimental data are same as in figure 7. Lines correspond to values computed using the GT-IC (solid lines) and CNT-IC (dashed lines) for n-heptane (C7), n-octane (C8), n-nonane (C9), and ndecane (C10).As EoS was used only PC-SAFT.Then, the least squares solution of the over-determined system of equations Fb = f gives vector and for the reference nucleation rate Figures 11 and 12 show the number of molecules in critical droplets depending on the inverse reduced temperature T c /T corresponding to the reference nucleation rates  J ref = 10 12 m −3 s −1 for figure 11 and J ref = 10 6 m −3 s −1 for figure 12, respectively.Numerical values are, as in figure 9, computed for all the substances C7-C10 using the GT-IC and CNT-IC combined with the PC-SAFT EoS.
Theoretical lines are based on the definition of the excess number of molecules.The whole number of molecules N in the critical droplet is the sum of the bulk volume number and the excess number [2], where V N * is volume of the droplet consisting of N * molecules.A rational option is to take the volume of the sphere of tension for the volume of the cluster, V N * = 4/3πr 3 s,N * and ∆N N * is the excess number of molecules of N * -molecule droplet.V N * and ∆N N * are connected with ρ G via (15).

01076-p.7
Experimental data marks uses the nucleation theorem proofed by Kaschiev [27]: therefore, in our case, number of molecules in critical cluster is given by N * = b 2 − 1, where b 2 is the second component of solution (25).Figures 9,10,11,and 12 show that the theoretical data describe the experiments quite satisfactorily; because of the range of the experimental data, Figs. 9, 11, correspond more to the data [19][20][21], Figs. 10, 12 to the rest.
Despite the agreement with the experimental data, there is a clear disagreement between their slopes.Temperature derivative of supersaturation and number of molecules is wrong especially for low temperatures.Therefore, there is certain temperature-dependent deviation which is still not satisfactorily resolved by the recent theory.

Conclusions
We computed nucleation rates of four n-alkanes: n-heptane, n-octane, n-nonane, and n-decane.We used traditional cubic PR EoS and modern accurate PC-SAFT EoS.We evaluated nucleation rates using the CNT and we solved the work of formation minimization problem using the GT.This led to numerical boundary value problem which we determined using a modified shooting method.During the process of solution, we had to develop several numerical enhancements because the problem was numerically unstable.We enhanced the nucleation rates expressions by the internally consistent correction given by [13].
We compared the nucleation rates theoretical data with the experiments.Best results were given by the combination of the GT with the PC-SAFT EoS; however, they were still not very accurate.
Using the Taylor expansion and simple linear interpolation we computed the ratio of GT-PC-SAFT nucleation rates and experimental values.This showed a temperaturedependence inconsistency between the theoretical and experimental data.Figures 9,10,11,and 12 showed this discrepancy even better.We argued that although the numerical and experimental data were quite in accord in these figures, there is a temperature drift especially to the low temperatures.
In future work, we would like to investigate and improve the temperature dependence of the theory and extend our work to the binary mixtures.

Fig. 1 .
Fig.1.Nucleation rates J of n-heptane as a function of supersaturation S computed for four different temperatures.Each symbol corresponds to one temperature.Nucleation rates were computed using the GT-IC (solid lines) and CNT-IC (dashed lines), compared to the experimental data (black symbols) by Hung et al.[18].Two equations of state are used: PR (white symbols) and PC-SAFT EoS (gray symbols).

Fig. 2 .
Fig. 2. Nucleation rates J of n-octane as a function of supersaturation S computed for seven different temperatures.Nucleation theories, EoSs, experimental data reference, line types and markers are same as in figure 1.

Fig. 6 .
Fig. 6.Nucleation rates J of n-decane as a function of supersaturation S computed for seven different temperatures.Nucleation theories, EoSs, experimental data reference, line types and markers are same as in figure 1.

Fig. 8
Fig.8.Ratio of nucleation rates J GT−IC evaluated using the expression(16) and nucleation rates J CNT−IC by Eq.(17).Temperatures and supersaturations correspond to the listed experimental data.Markers are same as in figure7.

Fig. 9 .
Fig.9.Dependence of the supersaturation S on the inverse reduced temperature T/T c (critical temperatures T c are listed in table 1) for nucleation rate J = 10 12 m −3 s −1 .Markers of the experimental data are same as in figure7.Lines correspond to values computed using the GT-IC (solid lines) and CNT-IC (dashed lines) for n-heptane (C7), n-octane (C8), n-nonane (C9), and ndecane (C10).As EoS was used only PC-SAFT.

Fig. 11 .
Fig. 11.Number of molecules in critical droplets.Markers of the experimental data are same as in figure 7, line types are same as in figure 9.

Table 1 .
Molar mass, critical temperature and pressure, and PC-SAFT parameters of considered substances Substance M [kg mol −1 ] T c [K] p c [bar] 1, b 2 : first, data was divided into groups corresponding to the measurements series.Let us concentrate on i-th series consisting of n measurements.Let us define matrix F and vec- )