Effective binary theory of multi-component nucleation V. I. Kalikmanov Citation: The Journal of Chemical Physics 142, 124111 (2015); doi: 10.1063/1.4916356 View online: http://dx.doi.org/10.1063/1.4916356 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/12?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Nucleation in binary polymer blends: Effects of foreign mesoscopic spherical particles J. Chem. Phys. 121, 1105 (2004); 10.1063/1.1761053 Multicomponent dynamical nucleation theory and sensitivity analysis J. Chem. Phys. 120, 9133 (2004); 10.1063/1.1695323 Nucleation in binary polymer blends: Effects of adding diblock copolymers J. Chem. Phys. 118, 8997 (2003); 10.1063/1.1566941 From binary and ternary to multicomponent nucleation: Atmospheric aerosol formation J. Chem. Phys. 115, 2641 (2001); 10.1063/1.1385157 Thermodynamics of the curvature effect on ice surface tension and nucleation theory J. Chem. Phys. 106, 1921 (1997); 10.1063/1.473329

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

THE JOURNAL OF CHEMICAL PHYSICS 142, 124111 (2015)

Effective binary theory of multi-component nucleation V. I. Kalikmanova) Twister Supersonic Gas Solutions, Einsteinlaan 20, 2289 CC Rijswijk, Netherlands and Department of Geoscience and Engineering, Delft University of Technology, Stevinweg 1, 2628 CN Delft, Netherlands

(Received 5 February 2015; accepted 17 March 2015; published online 31 March 2015) Classical theory of multi-component nucleation [O. Hirschfelder, J. Chem. Phys. 61, 2690 (1974)] belongs to the class of the so-called intractable problems: it requires computational time which is an exponential function of the number of components N. For a number of systems of practical interest with N > 10, the brute-force use of the classical theory becomes virtually impossible and one has to resort to an effective medium approach. We present an effective binary model which captures important physics of multi-component nucleation. The distinction between two effective species is based on the observation that while all N components contribute to the cluster thermodynamic properties, there is only a part of them which trigger the nucleation process. The proposed 2D-theory takes into account adsorption by means of the Gibbs dividing surface formalism and uses statistical mechanical considerations for the treatment of small clusters. Theoretical predictions for binary-, ternary-, and 14-component mixtures are compared with available experimental data and other models. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4916356] I. INTRODUCTION

saddle-point is

Understanding multi-component nucleation is of great importance for atmospheric and environmental sciences. Vivid examples are polar stratospheric clouds, acid rains, and air pollution. All these phenomena occur because Earth’s atmosphere is a multi-component gaseous system in which nucleation leads to formation of droplets of complex composition.1 A rich field of applications of multi-component nucleation is associated with the natural gas industry since nucleation is the primary mechanism responsible for formation of mist during the expansion of natural gas.2,3 Theoretical description of multi-component nucleation pioneered by Hirschfelder4 and Trinkaus5 represents the extension of the phenomenological Binary Classical Nucleation Theory (BCNT) of Reiss6 and Stauffer7 to N-component mixtures. Calculation of the nucleation rate requires determination of the critical cluster which is the saddle-point of the Gibbs free energy ∆G(n1, . . . , n N ) necessary to form an arbitrary Ncomponent cluster. The steady-state nucleation rate is given by4

aN , (2) ν where ν is the computer performance in FLOPS (number of floating point operations per second). This is an example of a so-called “intractable problem,” for which the calculation time grows exponentially with N. The fastest modern supercomputer Cray XK7 Titan in the Oak Ridge National Laboratory (USA)8 has a performance of νCray = 17.6 petaFLOPS (1 petaFLOP = 1015 FLOPS). Recall that for a single-component nucleation, a number of molecules in the critical cluster corresponding to appreciable nucleation rates lie in the range between several and several tens of molecules; a reasonable estimate of a can thus be a ≈ 50. For the multicomponent case, one has to probe the same range [0, a] along each of the N axis of the cluster size space. Let us make some estimates. In a typical natural gas composition, the number of components which trigger the nucleation process is about 20 or more.9 Setting in (2) N = 20, a = 50, we find that calculation on the Cray Titan supercomputer of just one nucleation point would require

J = J0e−∆G

∗/k T B

,

(1)

where T is the temperature, kB the Boltzmann constant, J0 is the kinetic prefactor, and the nucleation barrier ∆G∗ is the value of ∆G at the saddle-point. To find ∆G∗, one has to explore the N-dimensional space of cluster sizes. Let 0 ≤ nk ≤ a be a typical range of values of {nk } for appreciable nucleation rates. In order to determine the critical cluster for the given external conditions (pressure, temperature, and gas-phase composition), one has to calculate ∆G in a N points. Let us assume that calculation of ∆G at each point of the N-dimensional space requires one floating point operation. Then, the time required for determination of the a)Electronic address: [email protected]

0021-9606/2015/142(12)/124111/13/$30.00

t sp =

t sp =

5020 ≈ 5 × 1017 s, νCray

which is of the order of the age of the Universe t Univ ≈ 13.6 × 109 yr ≈ 4.3 × 1017 s.10 These estimates show that the effective-medium approach becomes a necessary (and unavoidable!) approximation. Within this approach, the original N-component system is replaced by a system of a lower dimensionality constructed in such a way that its nucleation behavior is close to the system in question. Having said this, one has to formulate a criterion which will ensure this proximity. In view of the exponential dependence of the nucleation rate on ∆G∗ (see Eq. (1)), an effective system to be defined should possess the

142, 124111-1

© 2015 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-2

V. I. Kalikmanov

same energy barrier as the original one. The next question to be answered is what is the minimal number of components in the effective system which correctly captures the physical behavior of the original system? In an attempt to answer this question, one might notice that the roles of various components in the nucleation behavior of the system can be qualitatively different. While some of them are sufficiently supersaturated to trigger nucleation, the others only influence the cluster properties being entrained into clusters due to unlike intermolecular interactions with the supersaturated species. In the early attempts to tackle this problem (see, e.g., Ref. 11), the simplest route was chosen: a multi-component system was replaced by an effective single-component model (1D approach). By doing so, the free energy of a cluster formation in the effective system becomes a function of only one variable—the total number n of effective particles in the cluster; the nucleation barrier corresponds to the maximum of ∆G(n). Constructing the 1D model, one makes quite a strong assumption that the composition of the N-component critical cluster coincides with the equilibrium composition of the bulk N-component liquid at the given nucleation conditions (pressure pv and temperature T). The intensive thermodynamic properties (surface tension, liquid density, etc.) are given by the corresponding properties of the original system at (pv,T)-equilibrium. It is important to emphasize that within the 1D approach, all components of the mixture are treated on the equal footing. Obviously, the 1D model tremendously enhances calculations of the nucleation rate making them practically analytical. The model demonstrates good agreement with experiment for sufficiently low pressures. Meanwhile, a number of nucleation experiments with binary mixtures clearly demonstrate that the composition of the critical clusters resulting from the experimental data by using the nucleation theorems12 can differ substantially from the equilibrium bulk liquid composition at the given (pv,T)13–16 leading to the failure of BCNT for these systems. The discrepancy between BCNT and experiment was vividly demonstrated in the of series of studies of n-nonane/methane nucleation carried out by Luijten et al.,17 Peeters et al.,18 and Labetski.19 It was shown that at sufficiently high pressures, the fraction of methane molecules in the critical cluster is much higher than the equilibrium bulk liquid fraction of methane at the same pv and T. The authors attribute this disagreement to a possible cluster structure consisting of a cluster core with predominantly nonane molecules and a surface layer containing predominantly methane molecules. The latter conjecture indicates that in order to capture physics of multicomponent nucleation, it is necessary to take into account adsorption effects which result in the surface enrichment. This observation goes back to the early work of Wilemski20 who proposed “a revised BCNT,” in which surface enrichment is taken into account under the capillarity approximation by distinguishing between the bulk and surface molecules of each species and using the Gibbs formalism of dividing surface. For a binary (and more generally, a multi-component) system, the Gibbs construction is not just a useful tool but a necessity, since it is impossible to choose a dividing surface for a mixture in such a way that adsorption terms for all species vanish.21

J. Chem. Phys. 142, 124111 (2015)

These considerations imply that for an appropriate description of nucleation in multicomponent system, one should ease the assumption of the fixed (a priori known) composition of the critical cluster. Speaking about an effective model, this means that the minimal number of components in the effective system should be two. Moreover, in view of the preceding discussion, it is necessary to discriminate between the bulk and surface composition of the effective species in the cluster. Another requirement for the effective model stems from the fact that classical theory of nucleation is based on the phenomenological capillarity approximation which assumes that the surface energy of a cluster can be described in terms of the plane layer surface tension. Obviously, for small clusters, this concept loses its meaning. This difficulty shows up already in the single-component case. In the recently formulated mean-field kinetic nucleation theory (MKNT),22 this problem is tackled by formulating an interpolative model between small clusters treated using statistical mechanical considerations and big clusters described by the capillarity approximation. Later, these considerations were extended to the case of binary mixtures.23 The paper is organized as follows. In Sec. II, we recall the main aspects of physics of multi-component nucleation and its classical treatment. In Sec. III, we formulate the criterion for construction of the effective binary system and discuss its thermodynamics. Section IV is devoted to statistics of effective binary clusters studied by means of the coarsegraining procedure; statistical mechanical consideration yields the cluster distribution function valid for all cluster sizes and arbitrary composition. Nucleation rate in the effective system and application of the model to various multi-component systems is discussed in Sec. V.

II. PHYSICS OF MULTICOMPONENT NUCLEATION A. Non-equilibrium state of N -component mixture: Metastability parameters

Consider a mixture of N components at an initially gaseous state characterized by the initial pressure p0, temperature T0, and mole composition y = ( y1, . . . , y N ), N k=1 yk = 1. After fast expansion, the system is brought to a nonequilibrium state—characterized by the pressure pv and temperature T, where nucleation starts. The nucleation point (pv,T; y) is located inside the vapor-liquid coexistence region; therefore, for the same pv and T, there exists the vapor-liquid equilibrium characterized by the equilibrium liquid and vapor compositions xeq = (x 1,eq, . . . , x N,eq) and yeq = ( y1,eq, . . . , y N,eq), respectively. In what follows, the subscript “eq” refers to the (pv,T)-equilibrium. We characterize the state of component i in the mixture by the metastability parameter Si = yi /yi,eq, i = 1, . . . , N,

(3)

which is an experimentally controllable quantity. We emphasize two important features: (i) Si depends on the behavior of all components in the mixture through the

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-3

V. I. Kalikmanov

J. Chem. Phys. 142, 124111 (2015)

equilibrium fractions yi,eq(pv,T; y) and (ii) Si does not depend on the (arbitrary) cluster composition.

(termed also the K-surface25) satisfying n 

l nexc i vi = 0,

(8)

i=1

B. Kinetics

Following classical considerations,4 we assume that formation of the multi-component n-cluster, n = (n1, . . . , n N ), results from attachment or loss of a single monomer E j of one of the species j of the cluster E(n1, . . . , n j , . . . , n N ), i.e., we consider the reactions of the type

where vil is the partial molecular volume of component i in the liquid phase. This particular choice of the dividing surface implies that ∆G(n) takes the form,26 Chap. 14 ∆G(n) = −

N 

ni ∆µi + Gsurf (n),

(9)

i=1 k

E(n1, . . . , n j , . . . , n N ) + E j k ′j E(n1, . . . , n j + 1, . . . , n N ). j

(4) Here, k j (n1, . . . , n N ) and k ′j (n1, . . . , n N ) are the forward and backward reaction rates, respectively. In the chemical kinetics, k j is proportional to the impingement rate ν j of component j and the surface area A(n) of the cluster: k j = ν j A(n). For the gas-liquid transition, ν j is given by the gas kinetics expression24 νj = 

y j pv 2πm j kBT

,

j = 1, . . . , N,

(5)

where m j is the mass of the molecule of component j. At each point of the N-dimensional cluster space, there is a nucleation flux J = (J1(n), . . . , JN (n)), each component of which is the net flux along the n j -axis resulting from the reactions of type (4), J j (n) = k j ϱ(n1, . . . , n j , . . . , n N ) − k ′j ϱ(n1, . . . , n j + 1, . . . , n N ), where ϱ(n1, . . . , n j , . . . , n N ) is the number density of (n1, . . . , n j , . . . , n N )-clusters. The backward rates k ′j are determined from the detailed balance condition at equilibrium when J = 0. For future purposes, it is useful to write down the impingement rates for (pv,T)-equilibrium conditions, ν j,eq = 

y j,eq p

where ∆µi = µvi − µli is the difference in the chemical potential of component i in the bulk mother phase (vapor outside the cluster) and in the bulk liquid phase having the same composition as the bulk composition of the cluster  N b surf x bi = nbi / i=1 ni . G is the surface free energy of the cluster, which for the sake of generality, we do not specify at the moment. Using the ideal gas approximation for the vapor and incompressibility of the liquid, we present ∆µi as1   1 yi pv , β= , β∆µi ≈ ln coex coex b k BT yi p (x ,T) where pcoex(xb,T) is the pressure of the saturated N-component vapor over the bulk N-component liquid with the molar composition xb and yicoex(xb,T) is the vapor fraction of component i at pcoex(xb,T) (the superscript “coex” emphasizes that the corresponding quantity refers to (xb,T) -equilibrium, rather than the (pv,T)-equilibrium). Equation (9) takes the form   N  y i pv g(n) = − ni ln coex coex b + g surf (n), (10) y p (x ,T) i i=1 where g = β ∆G(n) and g surf = β Gsurf (n). The terms in the square brackets can be written as ( yi pv)/( yicoex pcoex) = Si Ψi , where the dimensionless quantity Ψi (xb,T) =

v

2πm j k BT

,

j = 1, . . . , N.

(6)

i = 1, . . . , N.

(11)

describes the effect of the given cluster composition on ∆G. Equation (10) becomes

Then, the metastability parameters can be written as Si = νi /νi,eq,

yi,eq pv coex b yi (x ,T) pcoex(xb,T)

g(n) = gnoneq + geq,

(7)

C. Energetics of N -component cluster formation

An arbitrary n-cluster of the new phase (liquid) is immersed in the “mother phase” (supersaturated gas) at pv and T. Energetics of cluster formation determines the minimum reversible work ∆G(n) needed to form the n-cluster in the surrounding vapor. We introduce an arbitrary located Gibbs dividing surface21 distinguishing between the bulk (superscript “b”) and excess (superscript “exc”) molecules of each species in the cluster: ni = nbi + nexc i . Each of the quantities in the right-hand side depends on the location of the dividing surface, while their sum can be assumed independent of this location. Therefore, only ni ≥ 0 are observable physical properties; the model quantities nbi and nexc i can be of either sign. We locate the dividing surface at the equimolar surface for the mixture

(12)

where gnoneq(n) = −

N 

ni ln Si ,

(13)

ni ln Ψi + g surf .

(14)

i=1

geq(n) = −

N  i=1

geq is the free energy of an arbitrary n-cluster formation at the (pv,T)-equilibrium (when all Si = 1). In gnoneq, the metastability parameters Si do not depend on the cluster composition, while in geq, all Ψi and g surf are composition dependent. Only for the clusters, having the bulk composition coinciding with xeq(pv,T), all Ψi = 1, implying that geq = g surf (note that this is always the case for single-component nucleation). For all other compositions, Ψi , 1 and the sum in (14) can be of either sign.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-4

V. I. Kalikmanov

J. Chem. Phys. 142, 124111 (2015)

D. Steady-state nucleation rate

The total steady-state nucleation rate J is largely determined by the flux through the saddle-point n∗ of the Gibbs free energy;4 the latter corresponds to the critical cluster. Derivation of J is based on the “direction of principal growth” approximation, pioneered by Reiss,6 stating that the direction of the flux (not its absolute value) remains constant in the entire saddle-point region being equal to the direction at n∗. Applying this approximation, one obtains4 ∗ ∗ J = νav A∗ Z ∗ ϱcons eq (n ), ∗ where νav is the average impingement rate of the molecules on the critical cluster (per unit surface), A∗ is the surface area of the critical cluster, and Z ∗ is the Zeldovich factor; all quantities depending on the critical cluster have the superscript “*.” The function −g(n) ϱcons eq (n) = C e

is the number density of clusters at the so-called constrained equilibrium and C is its normalizing factor. The concept of constrained equilibrium, used in classical theory, refers to an equilibrium, which would exist for a vapor at the same pv, T, and metastability parameters Si as the vapor in question. In an alternative procedure formulated by Katz and coworkers27,28 for unary systems (and developed for mixtures in Ref. 26), called a “kinetic theory of nucleation,” the evaporation rate is obtained from the detailed balance condition at the (true) stable equilibrium of the saturated vapor (all Si = 1) at the same (pv,T). Within the kinetic approach, N  n∗ ∗ J = νav A∗ Z ∗ * Si i + ϱeq(n∗), , i=1 -

(15)

ϱeq(n) = Ceq e−geq(n)

(16)

where is the cluster distribution function at (pv,T)-equilibrium and Ceq is its normalizing factor. In particular, for a binary system, the average impingement rate is6 ν1 ν2 ∗ , (17) νav = ν1sin2 φ∗ + ν2cos2 φ∗ where the angle φ∗ denotes the direction of the principal growth in the (n1, n2)-space.

not contribute to the driving force (i.e., to gnoneq). At the same time, all components influence cluster properties in geq. As the simplest example, one can think here about a binary mixture in which one of the species is a carrier gas which at sufficiently high pressures cannot be considered inert. This observation suggests that the mixture components can be divided into two groups based on the values of Si : 1. “group µ”: Si > 1 + ϵ S , i = 1, . . . , M and 2. “group ν”: S j ≤ 1 + ϵ S , j = M + 1, . . . , N, where 0 < ϵ S ≪ 1 is a small parameter. For convenience, components belonging to the group µ will be called “supersaturated,” while those from the group ν will be called “saturated.”29 Consider an arbitrary n-cluster, in which first M components are supersaturated and the rest N − M components are saturated. The fractions of molecules in the cluster N are x k = nk /n, n = k=1 n . The number of molecules  k  within groups is n µ = µ ni , nν = ν n j , where we denote  M  N µ ≡ ν ≡ i=1, j=M+1. The fractions of components within each group can be expressed in terms of the equilibrium liquid fractions xi ni l l = = w µ,i + ∆w µ, i = 1, . . . , M, (18) i, nµ x µ xj nj l l = = wν, j = M + 1, . . . , N, (19) j + ∆wν, j , nν xν   where x µ = µ x i , x ν = ν x j , and the equilibrium probability weights are  x i,eq l w µ,i = , x µ,eq = x i,eq, (20) x µ,eq µ  x j,eq l , x = x j,eq. (21) wν, = ν,eq j x ν,eq ν l l and ∆wν, Equations (18) and (19) define the quantities ∆w µ,i j which are alternating in sign and satisfy M 

l ∆w µ,i =

i=1

The steady-state nucleation rate is an exponential function of g ∗, which implies that a plausible criterion to construct an effective system is the matching of its nucleation barrier to that of the original system.11 Nucleation barrier is the saddlepoint of ∆G(n1, . . . , n N ). Therefore, the implementation of the above formulated criterion requires mapping of the Ndimensional cluster to the cluster of lower dimensionality which (approximately) preserves ∆G. One can note that the roles of various components in the process of nucleation can be qualitatively different. While components with relatively high values of Si trigger nucleation, components with Si close to unity practically do

l ∆wν, j = 0.

(22)

j=M+1

Consider an arbitrary quantity f having the value f i for component i in the cluster. Constructing an effective binary model, we define averages of f within each group M  l l + w µ,i f i + * ∆w µ, i fi , i=1 , i=1 N N   + * l l ⟨ f ⟩νl = wν, ∆wν, j f j/ . j fj + . j=M +1 , j=M +1 In view of (22), expressions in the round brackets contain mutually compensating terms. Therefore, it is feasible to approximate averaging by using equilibrium weights (20) and (21)

⟨ f ⟩lµ =

III. THERMODYNAMICS OF THE EFFECTIVE SYSTEM

N 

⟨ f ⟩lµ

≈

M  i=1

M 

l w µ,i

f i,

⟨ f ⟩νl

≈

N 

l wν, j f j.

(23)

j=M+1

The superscript “l” indicates that averaging is with respect to the liquid fractions. For future purposes, we define in a similar

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-5

V. I. Kalikmanov

J. Chem. Phys. 142, 124111 (2015)

way the averages with respect to the vapor fractions ⟨ f ⟩vµ ≈

M 

v w µ, i f i,

i=1

N 

⟨ f ⟩νv ≈

v wν, j f j,

(24)

j=M+1

with v w µ, i =

yi,eq , y µ,eq

v wν, j =

y j,eq , yν,eq

(25)

  and y µ,eq = µ yi,eq, yν,eq = ν yi,eq. Using these considerations, gnoneq can be approximated as gnoneq ≈ −n µ ln Sµ − nν ln Sν ,

(26)

where Sµ and Sν are, respectively, the effective metastability parameters of µ- and ν-particles derived through Eq. (23), ln Sµ ≡ ⟨ln Si ⟩lµ ,

ln Sν ≡ ⟨ln S j ⟩νl .

(27)

The outlined procedure results in the dimensionality reduction of gnoneq: it becomes the function of only two variables—n µ and nν , while Sµ and Sν are fully determined by the input conditions in the original system and by the equation of state (EoS). This mapping preserves the total number of particles N in the cluster n ≡ i=1 ni = n µ + nν . Equilibrium part of the free energy can be written as (cf. Eq. (26)) geq ≈ −n µ ln Ψµ − nν ln Ψν + g surf ,

(28)

ln Ψµ ≡ ⟨ln Ψi ⟩lµ ,

(29)

where ln Ψν ≡ ⟨ln Ψ j ⟩νl .

Calculation of Ψµ and Ψν is not straightforward since Ψk ’s depend on the given n-cluster composition (as opposed to Sµ and Sν ). Moreover, according to Eq. (11), Ψk ’s, as well as other intensive properties, depend on the bulk (rather than overall) composition of the cluster. This observation brings us to the problem of inhomogeneity of particle distribution inside the cluster due to adsorption.

Here, µlµ and µlν are the chemical potentials of the µ- and ν-particles inside the cluster where the pressure is pl, γ(x νb ) is the plane-layer surface tension at (x νb ,T)-coexistence, and A is the cluster surface area. The solution of Eqs. (30) and (31) is given by Eqs. (B9) and (B10). Partial molecular volumes can be written as vkl = (1/ρl) η lk , k = µ, ν, where η lk is given by (B3). Equation (30) implies that A and the volume Vcl of the cluster can be expressed in terms of the total as well as the bulk quantities 1 1 (n µ η lµ + nν η νl ) = l (nbµ η lµ + nνb η νl ), (32) ρl ρ ( ) 2/3 ( ) 2/3 A = s1 n µ η lµ + nν η νl = s1 nbµ η lµ + nνb η νl , (33) Vcl =

where s1 = (36π)1/3(ρl)−2/3 is the average surface area of the effective monomer. The quantities Ψµ and Ψν result from averaging of ln Ψk within each group according to Eq. (29). In line with our approach, we present ykcoex using the equilibrium probability weights in the vapor phase, v yicoex ≈ w µ,i y µcoex,

i = 1, . . . , M,

y jcoex

j = M + 1, . . . , N.

≈

v wν, j

yνcoex,

Then, within each group, Ψk are the same for all components, yielding yν,eq pv y µ,eq pv Ψµ = coex b coex b , Ψν = coex b coex b . y µ (x ν ) p (x ν ) yν (x ν ) p (x ν ) (34)

B. Surface energy: Classical form and modified Parachor method

Within the classical phenomenological approach based on the capillarity approximation (valid for big clusters), the reduced surface energy of the N-component cluster takes the form surf gclass (n) = β γ(xb,T) A(n).

(35)

The surface tension γ(x ,T) can be derived by means of the modified Parachor method. The standard Parachor method reads30 N 
γ 1/4 = P[k] x coex ρl − ykcoex ρv , k b

A. Adsorption

Applying the formalism of the Gibbs dividing surface to the (n µ , nν )-cluster, we discriminate between bulk and excess particles of each (effective) species: n µ = nbµ + nexc µ , nν = nνb + nνexc. We characterize the bulk composition of a cluster by x νb = nνb /(nbµ + nνb ). For a given x νb , the solution of the (x νb ,T)-coexistence problem (outlined in Appendix A) yields the coexistence properties—pressure pcoex(x νb ,T), vapor composition y µcoex(x νb ,T), and liquid and vapor densities ρl(x νb ,T), ρv(x νb ,T). The form of Gibbs energy (12)-(14) requires that for the effective system, we use the K-dividing surface (cf. Eq. (8)) l exc l nexc µ v µ + nν vν = 0,

(30)

where v µl and vνl are partial molecular volumes of µ- and exc ν-particles in the liquid phase. To find nexc µ , nν for the given b b n µ , nν , we combine (30) with the Gibbs adsorption equation21 l l b exc l l b b nexc µ dµ µ (p , x ν ) + nν dµν (p , x ν ) + A dγ(x ν ) = 0. (31)

k=1

where P[k] is the Parachor of component k. We present it as 1/4 γ 1/4 = γ 1/4 µ + γν , where  M x coex  l b   * i + γ 1/4 = ρ x P[i] µ µ b  i=1 , x µ - ( coex )   M yi  − ρv y µcoex  P[i] coex  , yµ  i=1     N  x coex j γν1/4 = ρl x νb  P[ j] * b +  j=M +1 , x ν -    N y jcoex  − ρv yνcoex  P[ j] * coex + ,  j=M +1 , yν -

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-6

V. I. Kalikmanov

J. Chem. Phys. 142, 124111 (2015)

  and µ x coex = x bµ , ν x coex = x νb . Applying the averaging i j procedure to the terms in the square brackets, we introduce the generalized Parachors of components µ and ν in the liquid- and vapor phases according to Eqs. (23) and (24). The surface tension then becomes   γ 1/4(x νb ,T) ≈ ⟨P[i]⟩lµ x bµ ρl − ⟨P[i]⟩vµ y µcoex ρv   + ⟨P[ j]⟩νl x νb ρl − ⟨P[ j]⟩νv yνcoex ρv . (36) This expression is similar, but not identical, to the standard model. In the latter, P[k] are empirical material constants depending only on the chemical structure of a molecule. In (36), each effective component is characterized by the two generalized Parachors: one corresponding to the liquid and the other—to the vapor phase; moreover, both of them depend on the (pv,T)-equilibrium of the original system and thus effectively take into account all components in the mixture. Introducing the reduced surface tension θ ∞ = βγ s1, we present (35) as surf gclass = θ ∞ (n µ η lµ + nν η νl )2/3.

(37)

When clusters are small, containing several tens of particles, the phenomenological notion of the surface tension loses its meaning and one has to invoke statistical mechanical considerations.

IV. BINARY CLUSTER STATISTICS AT (pv , T )-EQUILIBRIUM

Having replaced the N-component system by the effective binary mixture, we consider vapor at (pv,T)-equilibrium as a collection of binary clusters. The partition function of an arbitrary (n µ , nν )-cluster is Z n µ n ν ≡ Zn =

1 3n

Λ µ µ Λν3nν

qn µ nν ,

energy of this gas using Stirling’s formula reads ) ( Kn . F (n) = Kn kBT ln Zn e Then, the chemical potential of an n-cluster µn = ∂F (n)/∂Kn is 3n  V Λ µ µ Λν3nν   . µn = kBT ln  ϱeq(n µ , nν ) (40) qn   Here, ϱeq(n µ , nν ) = Kn/V is the equilibrium cluster distribution function—the quantity we are aiming to determine in order to find nucleation rate (15). Equilibrium between the cluster and surrounding vapor requires µn = n µ µvµ,eq + nν µvν,eq.

(41)

Here, µvi,eq(pv,T) is the chemical potential of a particle of component i in the equilibrium vapor. Using (40) and (41), we find (q ) n [z µ,eq]n µ [zν,eq]nν , (42) ϱeq(n µ , nν ) = V β µv

where z i,eq = e i,eq/Λ3i . It is important to note that the chemical potential of a particle inside an arbitrary binary cluster is not the same as in the surrounding bulk equilibrium vapor—it depends on the cluster composition. Among the two effective species, ν refers to the group of saturated components which are usually in abundance in the original N-component mixture. Following considerations of Ref. 23, we rewrite qn in the coarse-grained form by integrating out the positions of ν-particles,  1 CG qn = dRn µ e−β H , (43) n µ ! cl where

(38)

H CG = Uµ µ ({Rn µ }) + Fν/µ ({Rn µ }; n µ , nν ,T)

where Λi is the thermal de Broglie wavelength of a particle of component i and qn µ nν is the configuration integral of the cluster in a domain of volume V ,  1 qn µ nν (T) = qn(T) = dRn µ drnν e−βUn, (39) n µ ! nν ! cl

is the Hamiltonian of the equivalent single-component cluster, and Fν/µ ({Rn µ }; n µ , nν ,T) = −kBT ln qν/µ is the free energy of ν-particles in the instantaneous environment of µ-particles,  1 qν/µ ({Rn µ }) ≡ drnν e−β(Uνν +U µν ). nν ! cl

where Rn µ and rnν are locations of particles µ and ν in the cluster, and Un = Uµ µ (Rn µ ) + Uνν (rnν ) + Uµν (Rn µ , rnν ) is the full interaction energy of the cluster  which includes µ-µ, ν-ν, and µ-ν interactions. The symbol cl denotes that integration is only over those configurations that belong to the cluster. The cluster as a whole can move through the entire volume V , while the particles inside it are restricted to the configurations about cluster’s center of mass that is consistent with a chosen cluster definition. The most probable configurations of the gas at low densities and temperatures are isolated clusters. It is then plausible to represent the equilibrium gaseous state of the binary mixture as a system of noninteracting clusters.31 The partition function of the gas of Kn of such clusters in the volume V is factorized: Z (n) = (1/Kn!) ZnKn, where the prefactor 1/Kn! takes into account the indistinguishability of clusters viewed as independent entities. The Helmholtz free

By construction, H CG guarantees that qn and hence all thermodynamic properties of the cluster are preserved. The diagrammatic expansion of ln qν/µ in the Mayer functions of µ − ν and ν − ν interactions results in the expansion of Fν/µ in m-body effective interactions between µ-particles,32,33 Fν/µ = F0(n µ , nν ,T) + U2({Rn µ }; x ν ,T) + · · ·.

(44)

The zeroth order contribution F0(n µ , nν ,T), called the volume term, does not depend on positions of the particles but is important for thermodynamics since it depends on cluster composition. The first-order (one-body) term in (44) vanishes in view of translational symmetry33 resulting in H CG = F0(n µ , nν ,T) + U CG({Rn µ }; x ν ,T) with the total coarse-grained interaction energy U CG = Uµ µ ({Rn µ }) + U2({Rn µ }; x ν ,T) + · · · . This yields the formally exact result qn = e−β F0 qnCG , µ

(45)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-7

V. I. Kalikmanov

J. Chem. Phys. 142, 124111 (2015)

where qnCG (x ν ,T) = µ

1 n µ!

Φν (x νb ),



CG

dRn µ e−βU .

Φi (x νb ) =

(46)

cl

Interpretation of Eqs. (45) and (46) is straightforward: by tracing out the degrees of freedom of ν-particles, we are left with the single-component cluster of pseudo-µ particles (its configuration integral is qnCG ) with the interaction energy U CG µ which implicitly depends on x ν . The volume term is F0 = Fν,id + ∆F0, where βFν,id = nν f ν,id is the free energy of ideal gas of ν-particles in the cluster volume Vcl and ∆F0 is the excess (over ideal) contribution due to ν − ν and µ − ν interactions; f ν,id(x νb ,T) depends only on intensive quantities. Calculation of ∆F0 in terms of interaction potentials is a very challenging task. Fortunately, for our purposes, we do not need to know its exact form. Instead, we make use of the general statement that ∆F0 is a homogeneous function of the first order in n µ and nν ,34 n µ nν f 1, β∆F0 = Vcl where f 1 is some unknown function of x νb and T. Using (32), we present it as β∆F0 = nν f 0, where f 0(x νb ,T) again depends only on intensive quantities. Thus, the volume term becomes e−β F0 = Φνnν ,

(48)

(51)

Substituting (51) into (49) and taking into account that   yi,eq pv z i,eq   = Ψi (x b ), ≈ ν z i,coex(x νb )  yicoex(x νb ) pcoex(x νb )  we obtain ϱeq(n µ , nν ) p v,T = Ceq e−geq,

(52)

geq = −n µ Ψµ (x νb ) − nν Ψν (x νb ) + g surf .

(53)

It comes as no surprise that the first two terms in (53) coincide with those of Eq. (28). The third term gives the desired surface contribution specified in Eq. (50). The phenomenological expression (37) can be written as surf gclass = θ ∞, µ n2/3 µ , where ( ) 2/3 xν l θ ∞, µ = θ ∞ η lµ + ην . xµ

(54) (55)

As it follows from Appendix C, Eq. (50) recovers the classical result (54) for big clusters (n µ → ∞).

To accomplish the thermodynamic description, it is necessary to identify the quantities entering g surf in Eq. (50). θ micro, µ is found from Eq. (C2) in which the saturation pressure can be identified with the partial pressure of the µ-component at (x νb ,T)-coexistence, psat = y µcoex(x νb ,T) pcoex(x νb ,T).

qnCG µ

can be analyzed using the formalism of the singleNow component MKNT22 outlined in Appendix C. Substitution of (C1) into (48) yields  n µ  nν surf ϱeq(n µ , nν ) p v,T = Ceq Φ µ (x νb )z µ,eq Φν (x νb )zν,eq e−g , (49) where   g surf (n µ ; x νb ) = θ micro, µ (x νb ) n sµ (n µ ; x νb ) − 1 .

i = µ, ν.

A. Parameters of the pseudo-µ system

(47)

where Φν (x νb ,T) is another unknown function. Substituting (45) and (47) into (42), we find qnCG   µ + [z µ,eq]n µ Φν (x b ) zν,eq nν . ϱeq(n µ , nν ) = * ν , V -

1 , z i,coex(x νb )

(50)

Here, θ micro, µ (x νb ) is the so-called reduced microscopic surface tension of the µ-system and n sµ is the average number of surface particles in the n µ -cluster (for details see Appendix C). In order to identify the unknown intensive quantities Φ µ (x νb ) and Φν (x νb ), let us consider (x νb ,T)-coexistence, or, equivalently, equilibrium at the pressure pcoex(x νb ,T) and temperature T. The cluster distribution function for this state has the form of Eq. (49) in which z i,eq is replaced by z i,coex = e β µ i,coex/Λ3i , where µi,coex(x b ) is the chemical potential at (x νb ,T)-coexistence. For the clusters falling on the x νb -equilibrium line, i.e., those whose bulk numbers satisfy, nνb = nlµ (x νb /x bµ ), the chemical potential of the molecule inside the cluster is equal to its value in the surrounding vapor at the pressure pcoex. The Gibbs formation energy of such a cluster will contain only the (positive) surface term: ϱ(n µ , nν ) p coex,T = Cµ,coex exp[−g surf ], yielding the identification of Φ µ (x νb ) and

(56)

Since the intermolecular potential in the pseudo-µ fluid is not known, one has to introduce an approximation for the second virial coefficient, which must depend on x νb . The simplest form satisfying the pure components limit is given by the mixing rule B2 = B2, µ µ (x bµ )2 + 2B2, µν x bµ x νb + B2,νν (x νb )2.

(57)

The quantities B2, µ µ , B2, µν , B2,νν , referring to the effective components, can be found from the virial coefficients of the pure components of the original system B2,ii (T) and the cross-terms B2,i j (T) (calculated using the standard combining rules30), B2, µ µ =

M 

l l w µ,i w µ, j B2,i j ,

B2,νν =

N 

l l wν,i wν, j B2,i j ,

i, j=M+1

i, j=1

B2, µν =

M  N 

l l wν, w µ,i j B2,i j .

i=1 j=M +1

Then, using (C2),  B y coex(x b ,T) pcoex(x b ,T)  2 µ ν ν  . θ micro, µ = − ln − k BT  

(58)

Calculation of n sµ requires the coordination number in the liquid phase N1. Let us introduce the average molecular

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-8

V. I. Kalikmanov

J. Chem. Phys. 142, 124111 (2015)

diameter in the pseudo-µ-system using the van der Waals one-fluid approximation M 

σ 3µ =

l l 3 w µ,i w µ, j σi j ,

i, j=1

where σii is the molecular size of component i and σi j = (σii + σ j j )/2. The average packing fraction of µparticles is η = (π/6) σ 3µ ρl(x νb ,T). Then, N1 is given by N1 ≈ 5.5116 η 2 + 6.1383 η + 1.275. Finally, n sµ is found from Eqs. (C4)-(C6) with the parameters given by Eq. (C7) in which θ ∞ must be replaced by θ ∞, µ and θ micro—by θ micro, µ . The prefactor Ceq is found from Eq. (C3), Ceq = β y µcoex(x νb ,T) pcoex(x νb ,T).

(59)

V. NUCLEATION RATE IN THE EFFECTIVE SYSTEM. RESULTS AND DISCUSSION

The definition of Si and the averaging procedure of Sec. III imposes the form of the impingement rates in the effective system. From (7) and (27),  l   νj l νi ln Sµ = ln , ln Sν = ln . νi,eq µ,eq ν j,eq ν,eq These expressions together with Eq. (5) suggest the introduction of the effective vapor fractions ln Yµ ≡ ⟨ln yi ⟩lµ,eq, ln Yν l ≡ ⟨ln y j ⟩ν,eq and the effective masses ln m µ ≡ ⟨ln mi ⟩lµ,eq, ln l mν ≡ ⟨ln m j ⟩ν,eq. Then, the effective impingement rates of components µ and ν are νµ = 

Yµ pv 2π m µ kBT

, νν = √

Yν pv . 2π mν kBT

(60)

The average impingement rate νav is given by Eq. (17). Combining the results of Secs. II–IV, we obtain   ∗ ∗ ∗ J = νav A∗ Z ∗ Ceq e−g , (61) where g ∗ is the saddle point in the (n µ -nν ) space of the Gibbs free energy     g(n µ , nν ) = −n µ ln Sµ Ψµ (x νb ) − nν ln Sν Ψν (x νb )   (62) + θ micro, µ (x νb ) n sµ (n µ ; x νb ) − 1 .

A. Binary mixture n-nonane/methane

An obvious “simplest” test of the proposed 2D-model is a binary mixture; in this case, each group contains just one component implying that all probability weights are equal to unity. We consider nucleation in the n-nonane/methane mixture extensively studied in shock tube experiments.17–19 N-nonane is a supersaturated (µ-)component, while the more volatile methane is the saturated (ν-)component. Experiments were carried out for the nucleation temperature T = 240 K, and the total nucleation pressures pv = 10, 25, 33, and 40 bars. Figure 1 shows the experimental results together with theoretical predictions of 3 models: BCNT, the coarse-grained theory of binary nucleation (CGNT) of Ref. 23, and the present 2D-model. It comes as no surprise that predictions of the two latter models are very close—for binary mixtures, the 2D-model by construction recovers the CGNT: both models take into account adsorption using the K-surface and treat small clusters without invoking the capillarity approximation. Both models agree with experiments within 1-2 orders of magnitude for most of the conditions except for the extremely low values of SC9 < 4.5 at the highest pressure of 40 bars. Classical theory gives reasonable results (within 2 orders of magnitude accuracy) for low pressures 10 and 25 bars, but drastically fails at higher pv. Analysis of the critical cluster content presented in Fig. 2 shows that at low pressures (pv = 10 bars), there are no methane molecules in it. At higher pv, the critical cluster has a core-shell structure, where almost all nonane molecules are bulk, while almost all methane molecules are excess; the number of bulk methane molecules does not exceed 2. One can also observe that at a fixed pv, the bulk,∗ ∗ methane content is insensitive to SC9, while nC9 and nC9 are decreasing functions of SC9. This results in the increase bulk,∗ ∗ of the methane fractions in the critical cluster, x C1 and x C1 , with SC9. A binary system has two thermodynamic degrees of freedom. Fixing pv and T, we completely determine all phase equilibrium properties, in particular x C1,eq(pv,T). Meanwhile, the composition of the critical cluster (and therefore its thermodynamic properties) is not only fixed by pv and T but also depends on SC9. In 1D-models of multicomponent nucleation (see, e.g., Ref. 11), it is assumed that clusters have composition (and properties) corresponding to (pv,T)equilibrium. One can see from the presented results that

FIG. 1. N-nonane/methane nucleation at T = 240 K and various total pressures. Nucleation rate as a function metastability parameter of n-nonane, SC9. Triangles: experiment17–19 and squares: 2Dmodel; the 2D-model points are calculated at the same SC9 values as the experimental points. Solid lines—BCNT and dashed lines—CGNT of Ref. 23. (a) p v = 10 bars and 25 bars; (b) p v = 33 bars and 40 bars.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-9

V. I. Kalikmanov

J. Chem. Phys. 142, 124111 (2015)

FIG. 2. Composition of the critical cluster according to the 2D-model: (a) ∗ , n∗ ; total numbers of molecules n C9 C1 (b) bulk numbers of molecules bulk,∗ bulk,∗ n C9 , n C1 .

this assumption is violated. At low pressures, the effect of this violation can be small, while at higher pressures, it is substantial, leading to large discrepancies between the predicted and measured nucleation rates. B. Ternary mixture n-nonane/propane/methane

Labetski19 performed a series of nucleation experiments in the ternary mixture n-nonane/propane/methane (C9/C3/C1). The nucleation rate data were obtained for pv = 40 bars and T = 240 K; in most experiments, the C3 vapor mole fraction was yC3 = 0.01. The vapor mole fraction of n-nonane was in the range 50 × 10−6 ≤ yC9 ≤ 78 × 10−6. In our calculations, we set the group division parameter ϵ S = 0.01. Group µ consists of n-nonane and group ν contains propane and methane. Figure 3 shows the experimental results of Ref. 19 and predictions of the 2D-model. Except for the lowest values of SC9 < 3.2, theoretical predictions are within 2 orders of magnitude accuracy compared to experiment. The model predicts the nucleation rate in the ternary mixture C9/C3/C1 to be higher than in the binary mixture C9/C1 at the same nucleation conditions pv = 40 bars, T = 240 K, and the same values of SC9 (cf. Fig. 1(b))—in agreement with

experimental observations of Ref. 19. Thus, the presence of propane enhances nucleation. Composition of the critical cluster and the nucleation barrier as a function of SC9 are shown in Fig. 4. As in the C9/C1 mixture at the same conditions, most of the supersaturated (nonane) molecules are bulk while most of the saturated (ν-) particles are excess. However, in the ternary mixture, the amount of bulk ν-particles is significantly higher than in the binary system. It is instructive to compare the predictions of the 2Dmodel to the 1D-model in which there is no distinction between the effective particles (no adsorption), and the cluster intensive properties are those of the N-component bulk liquid at (pv,T) equilibrium. The effective supersaturation is given by11 S1D =

N 

x

Si i,eq

(p v,T )

.

(63)

i=1

In Eq. (14), all Ψi = 1 and the Gibbs free energy reads   g1D (n) = −n ln S1D + θ micro n s (n) − 1 .

(64)

Nucleation barrier is given by the maximum of g1D (n). For the experimental conditions shown in Fig. 3, the effective supersaturation lies in the range 1.88 < S1D < 2.29. ∗ The critical cluster size is in the range 500 < n1D < ∗ 1200 corresponding to high barriers 170 < g1D < 350 and as a result—the unphysically low nucleation rates J1D ∼ 10−40 cm−3 s−1. These results reveal the crucial role in the analysis of the nucleation behavior played by the cluster composition. The neglect of composition and adsorption effects can lead to unphysical predictions.

C. 14-component mixture (natural gas sample)

FIG. 3. N-nonane/propane/methane nucleation at T = 240 K, p v = 40 bars, and yC3 = 0.01. Nucleation rate as a function metastability parameter of n-nonane, SC9. Closed triangles: experiment of Labetski;19 semi-closed squares: 2D-model.

Muitjens et al.2 studied nucleation of natural gas in a Wilson cloud chamber. Initially dry gas at the pressure p0 = 60 bars and the temperature T0 = 283 K was expanded to a state pv = 40 bars and the temperature T = 273 K where nucleation starts. Table I shows the composition of the natural gas sample used in Ref. 2 and individual metastability parameters of components at pv = 40 bars and T = 273 K. The estimated experimental nucleation rate is Jexpt ∼ 1010 − 1011 cm−3 s−1.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-10

V. I. Kalikmanov

J. Chem. Phys. 142, 124111 (2015)

FIG. 4. N-nonane/propane/methane nucleation at T = 240 K, p v = 40 bars, and yC3 = 0.01. Composition of the critical cluster according to the 2D-model. (a) Total numbers of molecules n ∗µ, nν∗ ; (b) bulk numbers of molecules n bulk,∗ , nνbulk,∗. Figure (a) µ also shows the nucleation barrier g ∗ as a function of SC9.

TABLE I. Composition of the 14-component mixture of Ref. 2 and metastability parameters at p v = 40 bars, T = 273 K. Component Nitrogen Methane Ethane Hexane Heptane Methylcyclohexane Toluene Octane Nonane Decane Undecane Dodecane Tridecane Tetradecane

yi

Si

1.58 × 10−1 8.13 × 10−1 2.84 × 10−2 6.53 × 10−5 8.62 × 10−5 2.49 × 10−5 3.19 × 10−5 5.39 × 10−5 4.22 × 10−5 2.41 × 10−5 8.02 × 10−6 1.73 × 10−6 3.98 × 10−7 1.00 × 10−7

1.000 1.000 1.000 1.008 1.026 1.024 1.033 1.080 1.243 1.715 3.122 6.915 18.017 41.444

Applying the same group division criterion ϵ S = 0.01, we find that at pv = 40 bars and T = 273 K, the upper four species in Table I—nitrogen, methane, ethane, and hexane—belong to the group of saturated (ν-) components, while the rest 10 hydrocarbons—from heptane down to tetradecane—belong to the group of supersaturated (µ-) components. The 2Dmodel predictions—effective metastability parameters, critical cluster composition, nucleation barrier, and nucleation rate— are shown in Table II. The predicted nucleation rate agrees with experiment within 2.5 orders of magnitude. As one can see, only 2 of 8 ν-particles are in the bulk of the critical cluster, while almost all µ-particles belong to the bulk. Application of 1D-models (63) and (64) to the same system ∗ ∗ yields S1D = 1.527, n1D = 1429, and g1D = 409.4 resulting in the unphysically low nucleation rate J1D ∼ 10−153 cm−3 s−1. Although S1D is only slightly different from Sµ , the difference in the critical cluster size and nucleation barrier is huge emphasizing the role of composition effects in multicomponent nucleation. Finally, note that brute-force calculations of the original 14-component mixture using the classical theory4 on the supercomputer Cray Titan8 would require about 1 yr computer time (for a single (pv,T)-point), while the effective 2D-model calculations on PC take several microseconds.

TABLE II. 2D-model results for the 14-component mixture of Table I at p v = 40 bars, T = 273 K; nucleation rate J is in cm−3 s−1. Sµ

Sν

n ∗µ nν∗ n bulk,∗ µ

1.719 1.000 18

8

19

n exc,∗ µ

nνbulk,∗

nνexc,∗

g∗

log10 J

−1

2

6

21.06

13.49

VI. CONCLUDING REMARKS

We have proposed an effective 2D-model of Ncomponent homogeneous vapor-liquid nucleation. The model introduces a division of mixture components into two groups: “supersaturated” and “saturated” species. While components from both groups contribute to the cluster thermodynamic properties, only the supersaturated components trigger the nucleation process. Mapping of an arbitrary N-component cluster onto the effective binary cluster takes into account adsorption in the effective system within the Gibbs dividing surface formalism resulting in the surface enrichment of saturated component. Treatment of effective binary clusters of arbitrary size is based on statistical mechanical considerations: applying the coarse-graining procedure to the configuration integral of the cluster, we trace out the degrees of freedom of the effective saturated component; the resulting coarsegrained cluster is treated without invoking the capillarity approximation using the formalism of the mean-field kinetic nucleation theory.22 The model does not contain adjustable parameters. For all systems studied, the 2D-model demonstrates a good agreement with available experimental data. It captures important physics of multi-component nucleation and opens the opportunity of solving the nucleation problem for systems with large number of components. ACKNOWLEDGMENTS

Support of Twister BV is acknowledged. I thank Twan Verweij for performing equilibrium flash calculations for multi-component mixtures using the “Multiflash” software and Bob Prokofiev for helpful and stimulating discussions. APPENDIX A: (xνb, T )-COEXISTENCE IN THE EFFECTIVE SYSTEM

Implementation of the nucleation model requires the solution of the (x νb ,T)-coexistence problem (known also as

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-11

V. I. Kalikmanov

J. Chem. Phys. 142, 124111 (2015)

bubble point calculations) in the effective system. The most widely used EoS (van der Waals, Peng-Robinson, and RedlichKwong-Soave) contains the energy parameter of the mixture a m (originating from intermolecular interactions) and the size parameter bm (responsible for excluded volume effects). If ak k and bk k are individual a and b parameters of components, then am =

N  N 

x i x j a i j , bm =

i=1 j=1

N  N 

x i x j bi j ,

(A1)

i=1 j=1

where x k are fractions either in the bulk liquid or in the bulk vapor phase. The unlike terms ai j , bi j are constructed by means of the mixing rules,30 ai j =

√

aii a j j (1 − k i j ),

bi j =

bii + b j j , 2

(A2)

where k i j is the so-called binary interaction parameter (empirical quantity); k ii = 0. In the standard equilibrium calculations parameters ai j , bi j , i, j = 1, . . . , N are the same in both coexisting phases. In the effective system, one has to discriminate between the bulk liquid and bulk vapor phase. Consider first the liquid phase, then from (A1) and (A2), the parameter bm in the bulk liquid reads   blm = x bi bii + x bj b j j ≈ x bµ blµ + x νb bνl , (A3) µ

ν

where blµ = ⟨bii ⟩lµ , bνl = ⟨b j j ⟩νl . In order to find a m for the effective system, let us first set all k i j = 0, which would correspond to the Lorentz-Berthelot mixing rules. Then, from (A1), 2

 √ √ x bi aii + x bj a j j +/ µ ν  b 2 l b b ≈ (x µ ) a µ + 2 x µ x ν alµ aνl + (x νb )2aνl ,

(alm )LB = *. ,



(A4)

  √ √ where alµ = ⟨ aii ⟩lµ , aνl = ⟨ a j j ⟩νl . Effective binary cluster contains two types of particles, and one can introduce the binary interaction parameter k lµν between the two species. It is then plausible to modify (A4) as  alm ≈ (x bµ )2alµ + 2 x bµ x νb alµ aνl (1 − k lµν ) + (x νb )2aνl . (A5) Parameter k lµν (pv,T) can be chosen from the following considerations. Let i lµ be the component which at given (pv,T) has the maximum liquid bulk fraction x i,eq among the supersaturated (µ-) components. Similarly, let jνl be the component which at given (pv,T) has the maximum liquid bulk fraction x j,eq among the saturated (ν-) components. Then, we set k lµν = k i lµ, jνl . The same considerations apply to the vapor phase with the composition y µcoex, bvm

≈

y µcoex bvµ

+ 

yνcoex bνv ,

(A6)

avm ≈ ( y µcoex)2avµ + 2 y µcoex yνcoex avµ aνv (1 − k vµν ) + ( yνcoex)2aνv , (A7)

where bvµ = ⟨bii ⟩vµ , bνv = ⟨b j j ⟩νv ,   √ √ avµ = ⟨ aii ⟩vµ , aνv = ⟨ a j j ⟩νv .

(A8) (A9)

The binary interaction parameter k vµν (pv,T) is found using the same considerations as outlined above but now applied to the equilibrium vapor phase: k vµν = k i vµ, jνv . Thus, calculation of (x νb ,T)-coexistence for the effective binary system is carried out using the same EoS as for the original Ncomponent mixture, but with EoS-parameters which depend on the phase and on the nucleation point (pv,T): bvµ , bνv , avµ , aνv , k vµν (vapor); blµ , bνl , alµ , aνl , k lµν (liquid).

APPENDIX B: K-SURFACE FOR THE EFFECTIVE BINARY CLUSTER: DETERMINATION OF EXCESS QUANTITIES

Combination of the K-surface Eq. (30) with the Gibbs adsorption Eq. (31) results in the set of two linear equations l exc l nexc µ v µ + nν vν = 0,

(B1)

l l b exc l l b nexc µ dµ µ (p , x ν ) + nν dµν (p , x ν ) + A dγ = 0.

(B2)

Partial molecular volume reads vil = (1/ρl) η li , where η li stands for the deviation of vil from its average value 1/ρl and ρl(x νb ,T) being the liquid number density. Standard thermodynamic considerations yield26 η lµ = 1 + x νb

∂ ln ρl ∂ ln ρl , η νl = 1 − x bµ , b ∂ xν ∂ x νb

(B3)

where x bµ = 1 − x νb . From the incompressibility of the liquid phase and the Laplace equation pl − pv = 2γ/r (where r is the cluster radius), we find µli (pl, x νb ) = µli (pv, x νb ) +

2γ η li , i = µ, ν. r ρl

(B4)

Excluding nexc µ from (B1), we obtain from (B2) and (B4)   * η νl +  2γ η νl * η νl +  nνexc  − dµ + dµ + d ln l  µ ν l  ηl  + A dγ = 0, r ρ η µµ - , , (B5) where for brevity, the notation µi ≡ µli (pv, x νl ) is used. Combining (B5) with the Gibbs-Duhem equation  x bi dµi = 0 i=µ,ν

and taking into account vil = (1/ρl) η li , we find nexc µ = −A

 l l l  −1 ∂γ  1 ∂ µ µ 2γ η µ ∂ ln(η µ /η ν )  +   , (B6) r ρl ∂ x νb  x νb η νl ∂ x νb ∂ x νb 

nνexc = −A

 l l  −1 ∂γ  1 ∂ µν 2γ η νl ∂ ln(η ν /η µ )  +   . (B7) r ρl ∂ x νb  x bµ η lµ ∂ x νb ∂ x νb 

The derivative ∂ µli /∂ x νb can be specified using a correlation for activity coefficients.30 To a good approximation, the activity

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-12

V. I. Kalikmanov

J. Chem. Phys. 142, 124111 (2015)

coefficients in the liquid can be set to unity, resulting in ∂ µlµ ∂ x νb

= −k BT

1 ∂ µlν 1 ; = kBT b . b b x µ ∂ xν xν

It is easy to see that these expressions satisfy the Gibbs-Duhem equation. The logarithmic derivatives in (B6) and (B7) are found from (B3) ∂ ln(η lµ /η νl ) ∂ x νb

=

1 (τ2 − τ12), η νl

η lµ

(B8)

where ∂ ln ρl ∂ 2 ln ρl , τ2 = . b ∂ xν ∂(x νb )2 Substituting (B8) into (B6) and (B7), we obtain

and ncore = 0 showing that such clusters do not possess liquid-like properties. For n ≥ N1 + 1, an n-cluster has a core containing on average ncore(n) particles and characterized by a radius Rcore and the number density ρl. The surface layer, containing on average n s (n) particles, is characterized by a thickness λr l, λ ≥ 1 (r l is the intermolecular distance in the liquid) and the constant density ξ ρl with 0 < ξ ≤ 1. Thus, for n ≥ N1 + 1, the actual smooth density profile is replaced by a two-step function with the middle step corresponding to the surface layer(s). Introducing the dimensionless core radius f = Rcore/r l, we have n s (n) = n − [ f (n)]3,

τ1 =

  −1 2γ ∂γ  kBT 2  (τ2 − τ1 ) . (B10) = −A b  b b l − r ρl η lµ  ∂ x ν  x µ x ν η µ

Here, the cluster surface area A(nbµ , nνb ) is given by Eq. (33). APPENDIX C: MAIN FEATURES OF THE MKNT22

MKNT22 is a semi-phenomenological single-component nucleation theory, nonperturbative in the cluster size, and as such is valid for arbitrary clusters. It is an interpolative model between small clusters treated using statistical mechanical considerations and big clusters described by the capillarity approximation. An arbitrary n-cluster is considered to consist of the two groups of molecules: the core ncore and the surface molecules n s : n = ncore + n s . The physical idea behind this distinction is that the core, if present, should possess the liquidlike structure which can be characterized by the coordination number in the liquid phase N1. Both ncore and n s fluctuate around their mean values ncore(n; T) and n s (n; T). Within the mean-field approximation, the configuration integral of the n-cluster in the physical volume V at the temperature T reads qn s = Ceq Φn e−θmicro [n (n)−1]. (C1) V Here, Φa = 1/zsat, zsat(T) is the fugacity at saturation; θ micro is the reduced free energy per surface particle, termed in Ref. 22 the “reduced microscopic surface tension.” MKNT allows to express the parameters in Eq. (C1) in terms of the equilibrium properties. It is shown that θ micro originates from nonideality of the vapor, ) ( B2 psat (C2) θ micro = − ln − kBT and psat . (C3) kBT Here, B2(T) is the second virial coefficient and psat(T) is the saturation pressure. For small clusters (n ≤ N1), Ceq =

n s (n) = n,

n ≤ N1

(C5)

where f (n) is the real positive root of the cubic equation

  −1 kBT 2γ ∂γ  exc 2  n µ = −A b − b b l + (τ2 − τ1 ) , (B9)  ∂ x ν  x µ x ν η ν r ρl η νl nνexc

for n ≥ N1 + 1,

(C4)

f 3 = −3ω f 2 − 3ωλ f + (n − ωλ 2).

(C6)

The parameters, λ and ω = ξ λ, of order unity, are found from the relationships  1 θ∞ N1 3 3 ω= , λ= − − . (C7) 3 θ micro ω 4 2 N1 is expressed in terms of the molecular packing fraction 3 in the liquid η = (π/6) ρl d hs , where d hs(T) is the effective hard sphere diameter of the molecule in the theory of liquids (see, e.g., Ref. 35). For 4 ≤ N1 ≤ 8, which is a typical range of coordination numbers in the liquid phase, N1 is well approximated by N1 ≈ 5.5116 η 2 + 6.1383 η + 1.275. By construction, for n → ∞, MKNT recovers the classical expression θ micro [n s (n) − 1] → θ ∞ n2/3. 1H.

Vehkamäki, Classical Nucleation Theory in Multicomponent Systems (Springer-Verlag, Berlin, 2006). 2M. J. Muitjens, V. I. Kalikmanov, M. E. H. van Dongen, A. Hirschberg, and P. Derks, Revue de l’Institut Français du Pétrole 49, 63 (1994). 3V. Kalikmanov, J. Bruining, M. Betting, and D. Smeulders, 2007 SPE Annual Technical Conference and Exhibition, Anaheim, California (USA), 11-14 November 2007 Paper No.: SPE 110736. 4O. Hirschfelder, J. Chem. Phys. 61, 2690 (1974). 5H. Trinkaus, Phys. Rev. B 27, 7372 (1983). 6H. Reiss, J. Chem. Phys. 18, 840 (1950). 7D. Stauffer, J. Aerosol Sci. 7, 319 (1976). 8http://www.olcf.ornl.gov/titan/. 9C. Luijten, R. van Hooy, J. Janssen, and M. E. H. van Dongen, J. Chem. Phys. 109, 3553 (1998). 10Planck collaboration, 2013. 11V. I. Kalikmanov and M. E. H. van Dongen, Phys. Rev. E 55, 1607 (1997). 12D. Kashchiev, Nucleation. Basic Theory with Applications (Butterworth, Oxford, 2000). 13P. Mirabel and J. L. Katz, J. Chem. Phys. 67, 1697 (1977). 14J. L. Schmitt, J. Witten, G. W. Adams, and R. A. Zalabsky, J. Chem. Phys. 92, 3693 (1990). 15R. Strey, P. E. Wagner, and Y. Viisanen, in Nucleation and Atmospheric Aerosols, edited by N. Fukuta and P. E. Wagner (Deepak, Hampton, 1992), p. 111. 16R. Strey and Y. Viisanen, J. Chem. Phys. 99, 4693 (1993). 17C. C. M. Luijten, P. Peeters, and M. E. H. van Dongen, J. Chem. Phys. 111, 8535 (1999). 18P. Peeters, J. Hruby, and M. E. H. van Dongen, J. Phys. Chem. B 105, 11763 (2001). 19D. G. Labetski, Ph.D. thesis, Eindhoven University, 2007. 20G. Wilemski, J. Chem. Phys. 80, 1370 (1984). 21J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity (Clarendon Press, Oxford, 1982). 22V. I. Kalikmanov, J. Chem. Phys. 124, 124505 (2006). 23V. I. Kalikmanov, Phys. Rev. E 81, 050601(R) (2010).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21

124111-13 24L.

V. I. Kalikmanov

D. Landau and E. M. Lifshitz, Statistical Physics (Pergamon, Oxford, 1969). 25A. Laaksonen, R. McGraw, and H. Vehkamaki, J. Chem. Phys. 111, 2019 (1999). 26V. I. Kalikmanov, Nucleation Theory (Springer, Dordrecht, 2013). 27J. L. Katz and H. Wiedersich, J. Colloid Interface Sci. 61, 351 (1977). 28J. L. Katz and M. D. Donohue, Adv. Chem. Phys. 40, 137 (1979). 29Group notations µ and ν should not be confused with the chemical potential and impingement rate notations.

J. Chem. Phys. 142, 124111 (2015) 30R.

C. Reid, J. M. Prausnitz, and B. E. Poling, The Properties of Gases and Liquids, 4th ed. (McGraw-Hill, NY, 1987). 31J. Frenkel, J. Chem. Phys. 7, 200 (1939); W. Band, ibid. 7, 324 (1939). 32M. Dijkstra, R. van Rooij, and R. Evans, Phys. Rev. Lett. 81, 2268 (1998). 33C. N. Likos, Phys. Rep. 348, 267 (2001). 34V. I. Kalikmanov, Phys. Rev. E 68, 010101 (2003). 35V. I. Kalikmanov, Statistical Physics of Fluids. Basic Concepts and Applications (Springer-Verlag, Berlin, 2001).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Mon, 18 May 2015 02:18:21