11778
J. Phys. Chem. B 2001, 105, 11778-11784
Temperature Dependence of a Kelvin Model for Binary Nucleation† Jin-Song Li and Gerald Wilemski* Department of Physics and Cloud and Aerosol Sciences Laboratory, UniVersity of Missouri-Rolla, Rolla, Missouri 65409-0640 ReceiVed: May 3, 2001; In Final Form: August 7, 2001
The temperature dependence of binary nucleation rates was studied using several different models for the cluster evaporation rates. The evaporation rates were found using the detailed balance relations in combination with three different models for the constrained equilibrium cluster concentrations: the classical Reiss model, the self-consistent classical model of Wilemski and Wyslouzil, and the new Kelvin model. The latter two models are modifications of the classical liquid droplet expression for the equilibrium cluster concentrations. The steady state binary kinetics equations were numerically solved using the technique of matrix inversion by patitioning. Numerous test calculations were made for the ideal ethanol-hexanol system. We found that the Kelvin model gave the most realistic temperature dependence of the three models examined. An analytical rate expression that accurately reproduces the numerical results is also presented.
I. Introduction Binary nucleation is important in many fields of science and technology ranging from atmospheric science to materials science.1 Since the pioneering work of Reiss2 and the major advances of Stauffer3 and Trinkaus,4 considerable progress has been made in understanding the kinetics of binary nucleation. Recent efforts include numerical and analytical calculations of binary nucleation rates for the steady state and transient cases.5-21 However, only a few papers address the experimental22-26 or theoretical27-30 temperature dependence of binary nucleation rates. The temperature dependence of unary homogeneous nucleation rates is known to be poorly predicted by classical nucleation theory,1,31 and a similar shortcoming is expected for the classical binary nucleation theory. In this paper we present an approximate Kelvin model32,33 to improve the predicted temperature dependence of the classical theory of binary nucleation. The essence of this model is the assumption that the evaporation rate of molecular species ν from a binary cluster is proportional to the equilibrium vapor pressure Pν of a binary liquid droplet of radius r and composition x as given by the Kelvin equation. The two Kelvin equations for a binary droplet have the form Pν/P∞ν (x) ) exp[σ(x)Vν(x)/(rkBT)]; P∞ν (x) σ(x), and Vν(x) are, respectively, the equilibrium vapor pressure of species ν, the surface tension, and the molecular volume of species ν for a bulk solution of composition x. Of course, the Kelvin equation is approximate.34-36 Despite this, Powles35 has found that for pure Lennard-Jones fluids the Kelvin equation with constant surface tension is surprisingly accurate even for drops containing only a few thousand atoms. A thousand atoms is still much larger than a typical critical droplet, but Viisanen et al.37,38 have found that the sizes of critical clusters with fewer than a hundred molecules are also well described by the Kelvin equation. Theoretical support for the ability of the Kelvin equation to predict small critical cluster sizes is found in the recent work of McGraw and Laaksonen,39-41 Talanquer,42 and Koga and Zeng.43 Although the magnitude of the predicted rate is not always improved using the Kelvin model for unary nucleation, the temperature dependence of the rate predicted †
Part of the special issue “Howard Reiss Festschrift”.
by the Kelvin32,33 and closely related44,45 models is better than that of other classical models. It is thus of interest to adapt this model to binary nucleation. We use the term “Kelvin model” for this adaptation because the earlier model for unary nucleation was also referred to as the “Kelvin approach” by Katz32 and the “Kelvin model” by Wilemski.33 We note that the Kelvin equation has also been used in a fundamentally different way by McGraw and Laaksonen39 to obtain a significant, purely temperature dependent correction to the classical work of critical cluster formation. In contrast to our approach, the magnitude of their correction term is not determined by the Kelvin equation and must be evaluated using other means.39,41-43 This can be advantageous, however, since it allows the inclusion of nonclassical effects. This approach has not yet been extended to binary systems. Our new binary model is still essentially classical and, thus, cannot be expected to predict a temperature dependence that is in full agreement with experimental measurements. Nevertheless, for unary systems the differences between the Kelvin model and the traditional classical model produce a beneficial temperature dependent correction to the classical rate. This is mainly a consequence of the asymptotic relationship46 between the classical model and the Kelvin model for the size dependence of the surface free energy term WS. The size variable here is i, the number of molecules in the cluster. In the classical capillarity approximation, WSC varies as i2/3 for all values of i, including the smallest clusters. In the Kelvin model, on the other hand, WS emerges as a sum of σV/r over discrete cluster sizes,33 WSK i ∝ ∑k)2 k-1/3, and WSK only approaches i2/3 for very large values of i. Even for large critical clusters with, for example, i ∼ 100, the difference between i2/3 and the Kelvin sum is significant (about 6%). Since this difference is multiplied by a temperature dependent factor, it is clear why the Kelvin model leads to a different temperature dependence for the rate. Hale has previously given a similar analysis.45 One other introductory comment is warranted. Some care is required to introduce the Kelvin equations into the binary kinetics formalism without running the risk of internal contradictions. For example, naive use leads immediately to violation of the principle of detailed balance, which should not be lightly disregarded.47,48 Since we rely on detailed balance to construct
10.1021/jp011690b CCC: $20.00 © 2001 American Chemical Society Published on Web 10/10/2001
Kelvin Model for Binary Nucleation
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11779
our evaporation rate coefficients, we are also bound to employ a physically reasonable constrained equilibrium cluster distribution. At a minimum, this requires our distribution function to obey the law of mass action which is a rigorous consequence of the second law of thermodynamics for ideal gas mixtures.
B. Kelvin Model. By applying detailed balance to the pairs of fluxes between any four nearest neighbor points forming a closed loop, we obtain the following rule,
EA(i, j - 1)
EB(i, j)
II. Binary Kinetics and Kelvin Equation
JA(i, j, t) ) ΓA(i, j)NAf(i, j, t) - EA(i + 1, j)f(i + 1, j, t) (1)
where i and j denote the numbers of the molecules in the cluster for A and B species, respectively; Nν is the number density of monomers of species ν ()A or B) in the vapor; f is a nonequilibrium cluster concentration, and Γν(i, j) and Eν(i, j) are the forward and reverse rate coefficients for the growth and decay of a cluster (i, j). The forward rate coefficient Γν(i, j) is determined by collision cross sections and mean molecular velocities, and it reads
x
2πmνm(i, j)
)
EA(i, j)
which is a special case of the kinetic product rules developed previously.12 This rule must be satisfied to eliminate the possibility of unidirectional currents around a loop at equilibrium.47 By using the Kelvin equation for the vapor pressure of a spherical drop to evaluate the Eν, we have the following two expressions for the evaporation rate coefficients:
JB(i, j, t) ) ΓB(i, j)NBf(i, j, t) - EB(i, j + 1)f(i, j + 1, t) (2)
kBT[mν + m(i, j)]
EB(i - 1, j)
1 (8)
A. Kinetics Equations. By assuming that the kinetic process involves only monomer addition or evaporation, the net fluxes for the conversion of clusters of size (i, j) to size (i + 1, j) or (i, j + 1) take the forms
Γν(i, j) ) π[dν + d (i, j)]2
ΓB(i - 1, j - 1) ΓA(i - 1, j)
ΓB(i, j - 1) ΓA(i - 1, j - 1)
(3)
where k is the Boltzmann constant, T is the temperature, dν is the diameter of a monomer of type ν, d (i, j) is the diameter of the cluster, mν is the molecular mass of the impinging species, and m(i, j) is the molecular mass of the cluster. To obtain the evaporation rates Eν(i, j), let us apply detailed balancing to eqs 1 and 2,
ΓA(i, j)NAN(i, j) ) EA(i + 1, j)N(i + 1, j)
(4)
ΓB(i, j)NBN(i, j) ) EB(i, j + 1)N(i, j + 1)
(5)
where N(i, j) is the equilibrium cluster concentration. Thus, once the equilibrium cluster concentration N(i, j) is known, by employing eqs 3-5, the evaporation rates Eν(i, j) can be obtained. The time dependence of f (i, j) is governed by
df (i, j, t) ) dt JA(i - 1, j, t) - JA(i, j, t) + JB(i, j - 1, t) - JB(i, j, t) (6) At steady state, df/dt ) 0, and JA and JB do not depend on time, so that eq 6 becomes
JA(i - 1, j) - JA(i, j) + JB(i, j - 1) - JB(i, j) ) 0 (7) As pointed by Wyslouzil and Wilemski,13 for the mixed dimer concentration, either JA(0, 1) or JB(1, 0) should be omitted from the right-hand-side of eq 6 in order to avoid doubling the mixed dimer formation flux, and EA(1, 1) ) EB(1, 1) ) E(1, 1). The nucleation rate can be obtained by summing all the fluxes JA and JB, which cross any arbitrary line joining the A and B axes. At steady state the value for the nucleation rate is a constant that does not depend on the particular line chosen.49 In what follows, we only consider the steady-state problem. Our goal is to solve the kinetics equations and obtain the nucleation rate at steady state.
RA(i, j) )
EA(i, j) ΓA(i - 1, j)NA
)
( (
) )
N∞A(i, j) 2σ(i, j)VA exp NA rkBT
N∞B (i, j) 2σ(i, j)VB ) exp RB(i, j) ) NB rkBT ΓB(i, j - 1)NB EB(i, j)
(9)
(10)
where N∞ν (i, j) denotes the vapor number density of species ν in equilibrium with a bulk mixture whose composition is that of the cluster (i, j). Substituting eqs 9 and 10 into eq 8, we obtain
[ [
] ]
[
]
N∞B (i, j) 2σ(i, j)VB N∞A(i, j - 1) 2σ(i, j - 1)VA exp exp ) NB N r (i, j)kBT r(i, j - 1)kBT A
[
2σ(i, j)VA N∞B (i - 1, j) 2σ(i -1, j)VB N∞A(i, j) exp exp NA NB r (i, j)kBT r(i - 1, j)kBT
]
(11)
The above equation is a consequence of detailed balance and the Kelvin equation, and it restricts the allowable values for the surface tension. For an ideal solution, we have
N∞A(xA) ) xAN∞A N∞B (xB) ) xBN∞B and eq 11 simplifies to
σ(i, j)VB r(i, j)
+
σ(i, j - 1)VA r(i, j - 1)
)
σ(i, j)VA r(i, j)
+
σ(i - 1, j)VB r(i - 1, j)
(12)
We have found by direct numerical evaluation that eq 12 is not satisfied if we use the properties of bulk mixtures to evaluate the various terms. Since the Kelvin equation and its thermodynamic basis are inadequate for molecular length scales, this is not surprising for small clusters. However, eq 12 is not satisfied even for very large droplets. Thus, this straightforward use of the Kelvin equation is inconsistent with the requirements of detailed balance. Suppose, however, we take a different point of view and simply use the Kelvin equation as a parametric form for the cluster evaporation rates. We can then conjecture that with proper boundary conditions, eq 12 can be solved to yield a set of self-consistent surface tension values that would satisfy detailed balance and would, presumably, have a significant dependence on the size of the cluster, in addition to its composition. For example, if the surface tensions for the pure liquid droplets and the mixed dimer are known as a function of size, eq 12 can give the binary surface tensions in the size space.
11780 J. Phys. Chem. B, Vol. 105, No. 47, 2001
Li and Wilemski
These “surface tension” values would then serve as a parametrization of the evaporation rates based on the form of the Kelvin equation. The detailed balance relations 4 and 5 may be iterated to obtain the following expression relating the equilibrium distribution N(i, j) and N(1, 1), j
N(i, j) ) N(1, 1)
i
1
∏ ∏ k)2 R (1, k) l)2 R B
1
(13)
A(l, j)
Equation 12 ensures that the above connection between N(1, 1) and N(i, j) is independent of the iteration sequence chosen. To get a self-consistent expression for N(1, 1), let us choose the harmonic mean of the product of the pure dimer number densities
N(1, 1) ) [NA(2, 0)NB(0, 2)]1/2
[
]
σAVA σBVB 2σ(1, 1)V(1, 1) ) + r(1, 1) r(2, 0) r(0, 2)
(15)
As noted elsewhere,12 the above expression for N(1, 1) does not appear to violate any physical or chemical principles, but the prefactor (N∞A N∞B )1/2 does not suffice in general, since the resulting expression for N(i, j)does not reduce properly to the unary distribution. Following the suggestion by Wilemski and Wyslouzil, we replace the prefactor (N∞A N∞B )1/2 with (N∞A)xA(N∞B )xB. Thus, by employing the Kelvin equation and eq 14 in eq 13, we obtain the following expression for N(i, j):
( )( )
NA i NB j(i + j)! 1 N(i, j) ) (N∞A)xA(N∞B )xB 2 N∞A N∞B i!j!
[
r(1, 1)kBT
i
-
2σ(l, j)VA
]
2σ(1, m)VB
∑ ∑ l)2 r(l, j)k T m)2 r(1, m)k T B
ΓA(i - 1, 0)
i)2
EA(i, 0)
∏
(17)
(Although we use the notation of the binary case, please remember that it is for the unary case.) The evaporation rate based on the Kelvin equation for a unary system is
B
(16)
Thus, with the proper surface tensions that are satisfied by eq 12, an equilibrium distribution based on this use of the Kelvin equation can be found, and the nucleation rate can be calculated. Unfortunately, at present we do not have data on the size dependent surface tension for pure liquids or, equivalently on the evaporation rates of pure liquid clusters, so that we cannot perform the calculation for the nucleation rate. Hence, we introduce a simplified Kelvin model in the next section that is a more practical working model at present. It should be noted that N(i, j) will not reproduce exactly the expressions for EA and EB given above, because of the change in the prefactor. However, this contradiction can be resolved by slightly modifying the expressions for EA and EB as proposed by Wilemski and Wyslouzil.12 Note also that this model result based on the Kelvin equation automatically includes the combinatorial contribution ((i + j)!)/ (i!j!), which ensures that for small clusters the ideal entropy of mixing term is not overestimated.50 III. Simplified Kelvin Model In this section we follow a different strategy that circumvents the difficulties of the direct approach. This strategy relies, first,
[
]
2σ(i, 0)VA r(i, 0)kBT
(18)
Substituting eq 18 into eq 17, setting r(i, 0) ) r0i1/3, and assuming that the surface tension is independent of the size, we obtain the following expression for the equilibrium cluster concentration:
N(i, 0) )
N∞A
() ( NA
i
N∞A
2σVA i exp k-1/3 kBTr0 k ) 2
∑
)
(19)
where r0 ) (3VA/4π)1/3. Equation 19 is the key result for the unary Kelvin model: The sum over cluster sizes is a consequence of treating the cluster formation kinetics discretely, and the k-1/3 size dependence follows directly from the Kelvin equation in eq 18. The sum in eq 19 can be well-approximated by the Euler-Maclaurin summation formula as33 i
j
-
j
(14)
where V(1, 1) ) (VA + VB)/2, and σ(1, 1)is defined by
2σ(1, 1)V(1, 1)
N(i, 0) ) (NA)i
EA(i, 0) ) ΓA(i - 1, 0)N∞Aexp
NA NB 2σ(1, 1)V(1, 1) ) (N∞A N∞B )1/2 ∞ ∞ exp r(1, 1)kBT NA NB
exp -
on developing an equilibrium cluster distribution that contains the characteristic size dependence of the Kelvin equation and, second, on using the detailed balance equations to derive a consistent set of cluster evaporation rates. Before introducing this simplified Kelvin model, let us first review the Kelvin model for unary nucleation.33 By employing the detailed balance condition, the equilibrium cluster concentration can be written as
3
1
1
∑ k-1/3 ) 2 (i2/3 - 1) - 2 (1 - i-1/3) + 36 (1 - i-4/3) k)2
(20)
Substituting the above results into eq 19 leads to the following equilibrium distribution for the Kelvin model:33
N(i, 0) ) N∞A
() { NA
N∞A
i
1 exp -Θ(i2/3 - 1) - Θ - (1 - i-1/3) + 3
[
1 (1 - i-4/3) 54
]} (21)
where Θ ) 4πr02σ/kBT. The above result differs from the classical equilibrium cluster concentration NC(i, 0) found by Courtney,51
[
NC(i, 0) ) N∞A exp -Θi2/3 + i ln
( )] NA
N∞A
(22)
in that the size dependence of the surface term i2/3 is replaced by
1 1 (i2/3 - 1) - (1 - i-1/3) + (1 - i-4/3) 3 54
(23)
The first term in expression 23 gives the correction included in the self-consistent classical equilibrium cluster concentration,32,33,44 and the second term is a correction that appears in Hale’s analysis of small cluster contributions in scaled nucleation theory.45 Hence, the Kelvin model automatically gives rise to
Kelvin Model for Binary Nucleation
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11781
some important corrections to the classical equilibrium distribution. Since these additional terms are multiplied by the temperature-dependent factor Θ, they give rise to the different and, in fact, improved temperature dependence of the Kelvin model compared to the usual classical model. Using eq 23, we may write the Kelvin surface free energy as
WSK ) WSC - (71/54)Θ + i-1/3(1 - i-1/18)Θ/3 which shows that the main correction to the classical result is a purely temperature-dependent factor.45 The other source of difference arises from the prefactors of the two models. We now extend this result to binary systems. The classical equilibrium cluster concentration for binary nucleation was proposed by Reiss,2
NR(i, j) ) (NA + NB)
( )( ) NA ∞ NA(xA)
i
NB ∞ NB (xB)
j
exp[-Θ(xB)n2/3] (24)
where n ) i + j, Θ(xB) ) σ(xB)s(xB)/kBT, and s(xB) ) (36π)1/3 (xAVA + xBVB)2/3. This expression violates the mass action law, and it does not reduce to the cluster concentration for pure A and B species when i or j are separately equal to zero.52 Recently, Wilemski and Wyslouzil12 proposed a self-consistent equilibrium binary cluster concentration that satisfies the law of mass and action and reduces to a self-consistent classical concentration for unary clusters. With the inclusion of the combinatorial expression that ensures that the internal ideal entropy of mixing is not overestimated,18 it reads
NSCC(i, j) ) (N∞A)xA(N∞B )xB
( )( ) NA
N∞AγA
i
NB
N∞B γB
(i + j)! i!j!
j
Figure 1. Nucleation rates for the Kelvin model as a function of temperature at different ethanol-hexanol gas-phase activities: aE ) NE/N∞E and aH ) NH/N∞H.
TABLE 1: Physical Properties of Ethanol and Hexanol53 density (g/cm3) FE ) 1.059 - 1.0044 × 10-3T + 2.9 × 10-7T2 FH ) 1.0948 - 9.5928 × 10-4T saturation pressure (Pa) log P∞E ) 11.885 - 2371.0/T log P∞H ) 12.196 - 3045.0/T surface tension (dyn/cm) σ(T, xB) ) σA(T) + [σB(T) - σA(T)]xB σE ) 47.188 - 0.085T σH ) 51.901 - 0.088471T
exp(xAΘA + xBΘB) exp[-Θ(xB)n2/3] (25) where γν is the liquid-phase activity coefficient of species ν in a bulk mixture. We notice that a n2/3 term appears in the surface term that is similar to the unary case. Accordingly, we propose to replace n2/3 in eq 25 with expression 23 to obtain a simplified Kelvin model,
{
]}
1 1 N(i, j) ) NSCC(i,j)exp Θ(xB) (1 -n-1/3) - (1 - n-4/3) 3 54 (26)
[
Since NSCC(i, j) already includes the contribution from the selfconsistent distribution, we exclude -1 from the (n2/3 - 1) term. Although the equilibrium distribution in eq 26 is obtained by more heuristic means than that in eq 16, it is readily implemented in practical calculations through eqs 4 and 5. IV. Results and Discussion Numerical solutions of eq 7 were carried out using the method described in the Appendix. We performed the calculations for the ideal ethanol/hexanol system. The temperature-dependent physical properties for this system are shown in Table 1. The nucleation rates are calculated over the temperature range 220 e T e 300 K. We performed the calculations for the classical Reiss distribution, eq 24, the self-consistent classical (SCC) model with combinatorial prefactor, eq 25, and the Kelvin model, eq 26. To isolate the effect of the Kelvin correction terms, we also included the combinatorial prefactor in the Reiss model. This makes the volume term in the reversible work the same in all three models.
Figure 2. Ratio of SCC and Reiss model nucleation rate to that of the Kelvin model as a function of temperature at different gas-phase activities of ethanol and hexanol.
Figure 1 shows the numerically derived nucleation rates for the Kelvin model. All three models show the same qualitative dependence of nucleation rate with temperature. The nucleation rates for the Kelvin model are always slightly higher than those of the SCC model. Since there are no SCC and Kelvin corrections in it, rates for the Reiss model are always the lowest for these three models. It is more instructive to examine the ratio of nucleation rates for the different models at various temperatures and supersaturations. Figure 2 contains some examples of the ratios J/JKelvin of nucleation rate J (of either the SCC or
11782 J. Phys. Chem. B, Vol. 105, No. 47, 2001
Li and Wilemski
Reiss model) to the nucleation rate JKelvin of the Kelvin model for different combinations of vapor activities. The SCC model gives nucleation rates that are only slightly more sensitive to temperature than those of the Kelvin model; the ratio of rates varies by about a factor of 5 over the full temperature range. By contrast, the rates from the Reiss model exhibit a much stronger temperature dependence, with the ratio varying by about 2.5 orders of magnitude over the temperature range shown. The qualitative trend is identical to that found by Zeng and Oxtoby54 when they compared classical nucleation rates for unary nucleation with nonclassical rates based on density functional theory. In the present case, however, since the Kelvin model is still inherently classical and is missing important nonclassical features, the size of the improvement is several orders of magnitude smaller than that found by Zeng and Oxtoby.54 Besides improving the predicted temperature dependence of the classical model, the Kelvin model also yields much larger nucleation rates than the Reiss model. The same terms responsible for the improved temperature dependence also raise the nucleation rate, although this tendency is obviously moderated at higher temperatures. This trend was already apparent in the earlier unary and binary SCC models.12,32,33,44 The Kelvin rates are roughly 3-6 orders of magnitude higher than the Reiss rates, depending on the temperature, but only about 3-40 times greater than the SCC rates. In view of the labor required to calculate nucleation rates by the numerical method described in the Appendix, it is desirable to have an accurate analytical expression with which to calculate rates more easily. On the basis of past experience, Stauffer’s rate expression should be a very satisfactory analytical approximation.13 To verify this, we have compared our numerical nucleation rates with those calculated analytically using Stauffer’s formula,3 which may be written as13
J ) RavZN(i*, j*)
(27)
where
Figure 3. Comparison of theoretical rates with experimental rates55 at T ) 260 K for different gas-phase activities of ethanol and hexanol: x ≡ aH/(aE + aH) and a )
xaE2+aH2.
W(1)(i, j) ) ln
[ ] ( ) ( )
NA NB i!j! - i ln ∞ - j ln ∞ + Θ(xB)n2/3 (i + j)! NAγA NB γB
and W(2) arises from the present Kelvin model,
1 1 W(2)(i, j) ) Θ(xB) - (1 - n-1/3) + (1 - n-4/3) 3 54
[
]
We found that the analytical results were within 5% of the numerical nucleation rates for the numerous cases that we examined. We also found that the nucleation rate can be well approximated by
[
J ) JSCC exp Rav )
ΓA(i*, j*)ΓB(i*, j*)NANB ΓA(i*, j*)NA sin2 φ + ΓB(i*, j*)NB cos2 φ
Z ) -(WAAcos2 φ + 2WAB cos φ sin φ + WBB sin2 φ)/ (WAB2 - WAAWBB)1/2 tanφ ) s + xs2 + rS 2s ) (rSWBB - WAA)/WAB rS )
ΓB(i*, j*)NB ΓA(i*, j*)NA
where (i*, j*) is the saddle point on the reversible work surface, and WAA, WAB, and WBB are the second derivatives of the reversible work at the saddle point. For N(i*, j*) we used eq 26. We used the following expression for the reversible work:
Wrev(i, j) ) W(1)(i, j) + W(2)(i, j) where W(1) was used in earlier work,18
(28)
]
W(2)(i(1)*, j(1)*) kBT
(29)
where JSCC is the nucleation rate for the SCC model, eq 25, and (i(1)*, j(1)*) is the saddle point of W(1)(i, j). This change simplies the computation of the saddle point. The nucleation rates predicted by eq 29 were found to be within 10% of the numerical nucleation rates. Up to now, there is no experimental work on the temperature dependence of nucleation rates for the ethanol-hexanol system, although Strey and Viisanen have measured rates at a fixed temperature of 260 K.55 Figure 3 shows a comparison of experimental rates55 at different gas-phase activities with theoretical rates calculated using eq 27 with eq 28. The theoretical rates are 3-6 orders of magnitude higher than the experimental rates. The slopes of ln J agree fairly well with those found experimentally indicating that the activity or supersaturation dependence of the rate is properly described by the model. The trends seen here are consistent with the critical activity plot published previously for the SCC model.12 Acknowledgment. This work was supported by the Engineering Sciences Program of the Division of Materials Sciences and Engineering, Basic Energy Sciences, U.S. Department of Energy. Appendix: Computation Methods The set of kinetic eqs 7 can be solved numerically subject to proper boundary conditions. Here, these conditions are that the
Kelvin Model for Binary Nucleation
( )
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11783
ck11
Ck )
(
ck22 ···
(33)
ckkk
mk11 mk12 mk22 mk23 Mk ) ··· ··· k mkkk mk,k+1
)
(34)
k ) ΓB(k - i - 1, i - 1)NB pi,i-1
c322 ) -[(ΓA(1, 1)NA + ΓB(1, 1)NB + E(1, 1)] ck11 ) -[(ΓA(k - 1, 0)NA + ΓB(k - 1, 0)NB + EA(k - 1, 0)] Figure 4. Variable labeling scheme and the boundary conditions used in solving the binary kinetics equations. The cluster concentrations are labeled in the order shown by the circled numerals on the lattice sites.
ckkk ) -[(ΓA(0, k - 1)NA + ΓB(0, k - 1)NB + EB(0, k - 1)]
pkii ) ΓA(k - i - 1, i - 1)NA monomer concentrations are constant and that no clusters larger than a maximum size are present, i.e., f (i, j)|i+j)nmax+1 ) 0. In our calculations, we generally used nmax ) 200. We ascertained in several cases that virtually identical solutions were found using nmax ) 400. The maximum critical cluster size in our calculations was about n* ) 106, so an outer boundary of 200 suffices to simulate an infinite grid. We employ the technique of matrix inversion to solve the kinetic equations. The arrangement for the variables, shown in Figure 4, determines the structure of the matrix, and proper organization of the variables may greatly reduce the work for the calculation. We use f kl to characterize the variables, where k ) i + j + 1 and l ) j + 1, i.e., the variable vector F is
F)
( ) () F3 F4 l
Fnmax+1
(
,
f k1 fk Fk ) 2 l f kk
(30)
)
As result of this arrangement, the kinetic coefficient matrix C has the following structure,
C3 M3 P4 C4 M4 P5 C5 M5 C) ··· ··· ··· nmax P Cnmax Mnmax Pnmax+1 Cnmax+1
( )
mkii ) EA(k - i + 1, i - 1) k mi,i+1 ) EB(k - i + 1, i - 1)
The blank cells in matrices C, Pk, Ck and Mk are zeros. Thus, our problem is to solve the following matrix equation,
C·F )B where
()
B3 B) 0 , l 0
(35)
()
b1 B3 ) b 2 b3
b1 ) -ΓA(1, 0)NA2 b2 ) -ΓB(1, 0)NANB b3 ) -ΓB(0, 1)NB2
(31)
where Pk, Ck, and Mk are k×(k - 1), k×k, and k×(k + 1) submatrices, respectively, and they have the following forms,
pk11 pk21 pk22 pk32 ··· Pk ) ··· pk k-1,k-1 i pk,k-1
ckii ) -[(ΓA(k - i, i - 1)NA + ΓB(k - i, i - 1)NB + EA(k i, i - 1) + EB(k - i, i - 1)],(k * 3; k > i > 1)
Equation 35 can be solved by calculation of the inverse matrix C-1. Most of the matrix elements of C are zero. It is inefficient to use general inversion techniques on such a problem because most of the arithmetic operations involve zero operands. Furthermore, it is also wasteful to reserve storage for the zero elements of a very big matrix C. However, inversion by partitioning is an efficient way to solve this problem.56 The matrix C can be partitioned into
(
C C C ) C 33 C 3,K-3 K-3,3 K-3,K-3 (32)
)
nmax+1 where K ) ∑l)3 l is the size of square matrix C; C33 is the same as the square matrix C3; and C3,K-3, CK-3,3, and CK-3,K-3 are the other parts of matrix C: i.e., they are matrices of size 3×(K - 3), (K - 3)×3, and (K - 3)×(K - 3), respect-
11784 J. Phys. Chem. B, Vol. 105, No. 47, 2001
Li and Wilemski
ively. CK-3,K-3 is supposed to be partitioned further as
(
C C4,K-3-4 CK-3,K-3 ) C 44 C K-3-4,4 K-3-4,K-3-4
)
... Cnmax+1,nmax+1 ) Cnmax+1 The inverse of C is also partitioned in the same manner:
(
T T C-1 ) T 33 T 3,K-3 K-3,3 K-3,K-3
)
then T33, T3,K- 3, TK-3,3, and TK- 3,K-3 have the same sizes as C33, C3,K-3, CK-3,3, and CK-3,K-3, respectively. We continue to partition the inverse of CK-3,K-3 as
(
T T4,K-3-4 CK-3,K-3-1 ) T 44 T K-3-4,4 K-3-4,K-3-4
)
By employing C ‚ C-1 ) 1, CK-3,K-3 ‚ CK-3,K-3-1 ) 1, ..., we can find T33, T44, ..., Tnmax+1,nmax+1 as
T33 ) (C33 - M3T44P4)-1 T44 ) (C44 - M4T55P5)-1 ... Tnmax,nmax ) (Cnmax,nmax - MnmaxTnmax+1,nmax+1Pnmax )-1 +1
Tnmax+1,nmax+1 ) Cnmax+1,nmax+1-1 Since Tnmax+1,nmax+1 can be obtained by directly inverting the matrix Cnmax+1,nmax+1, one successively obtains Tnmax,nmax, ..., T44, and T33. Then the nonequilibrium cluster concentrations can be found from the following formulas
F3 ) T33B3 F4 ) T44P4F3 ... Fnmax+1 ) Tnmax +1,nmax+1Pnmax+1 Fnmax Once the nonequilibrium cluster concentrations in the size space are determined, the nucleation rate can be calculated by summing all the fluxes JA and JB, which cross any arbitrary line joining the A axis to the B axis. References and Notes (1) Laaksonen, A.; Talanquer, V.; Oxtoby, D. W. Annu. ReV. Phys. Chem. 1995, 46, 489. (2) Reiss, H. J. Chem. Phys. 1950, 18, 840. (3) Stauffer, D. J. Aerosol Sci. 1976, 7, 319. (4) Trinkaus, H. Phys. ReV. B 1983, 27, 7372.
(5) Greer, A. L.; Evans, P. V.; Hamerton, R. G.; Shangguan, D. K.; Kelton, K. L. J. Cryst. Growth 1990, 99, 38. (6) Shi, G.; Seinfeld, J. H. J. Chem. Phys. 1990, 93, 9033. (7) Zitserman, V. Yu.; Berezhkovskii, L. M. J. Colloid Interface Sci. 1990, 140, 373. (8) Wu, D. T. J. Chem. Phys. 1993, 99, 1990. (9) Nishioka, K.; Fujita, K. J. Chem. Phys. 1994, 100, 532. (10) McGraw, R. J. Chem. Phys. 1995, 102, 2098. (11) Berezhkovskii, L. M.; Zitserman, V. Yu. J. Chem. Phys. 1995, 102, 3331. (12) Wilemski, G.; Wyslouzil, B. E. J. Chem. Phys. 1995, 103, 1127. (13) Wyslouzil, B. E.; Wilemski, G. J. Chem. Phys. 1995, 103, 1137. (14) Wyslouzil, B. E.; Wilemski, G. J. Chem. Phys. 1996, 105, 1090. (15) Li, J.-S.; Nishioka, K.; Maksimov, I. L. Phys. ReV. E. 1998, 58, 7580. (16) Li, J.-S.; Nishioka, K. Chem. Phys. Lett. 1998, 259, 211. (17) Demo, P.; Kozˇ´ısˇek, Z.; Sˇ a´sˇik, R. Phys. ReV. E. 1999, 59, 5124; 60, 4995. (18) Wyslouzil, B. E.; Wilemski, G. J. Chem. Phys. 1999, 110, 1202. (19) Wilemski, G. J. Chem. Phys. 1999, 110, 6451. (20) Kevrekidis, P. G.; Lazaridis, M.; Drossinos, Y.; Georgopoulos, P. G. J. Chem. Phys. 1999, 110, 8010. (21) Li, J.-S.; Maksimov, I. L.; Wilemski, G. Phys. ReV. E. 2000, 61, R4710. (22) Schmitt, J. L.; Whitten, J.; Adams, G. W.; Zalabsky, R. A. J. Chem. Phys. 1990, 92, 3693. (23) Wyslouzil, B. E.; Seinfeld, J. H.; Flagan, R. C.; Okuyama, K. J. Chem. Phys. 1991, 94, 6827. (24) Wyslouzil, B. E.; Seinfeld, J. H.; Flagan, R. C.; Okuyama, K. J. Chem. Phys. 1991, 94, 6842. (25) Anisimov, M. P.; Koropchak, J. A.; Nasibulin, A. G.; Timoshina, L. V. J. Chem. Phys. 1998, 109, 10004. (26) Doster, G. J.; Schmitt, J. L.; Bertrand, G. L. J. Chem. Phys. 2000, 113, 7197. (27) Lazaridis, M.; Drossinos, Y. J. Phys. A: Math. Gen. 1997, 30, 3847. (28) Noppel, M. J. Chem. Phys. 1998, 109, 9052. (29) Vehkama¨ki, H.; Ford, I. J. J. Chem. Phys. 2000, 113, 3261. (30) Vehkama¨ki, H.; Ford, I. J. J. Chem. Phys. 2001, 114, 5509. (31) Barrett, J. C. J. Phys.: Condens. Matter 1997, 9, L19. (32) Katz, J. L. Pure Appl. Chem. 1992, 64, 1661. (33) Wilemski, G. J. Chem. Phys. 1995, 103, 1119. (34) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Oxford University Press: New York, 1982. (35) Powles, J. G. J. Phys. A: Math. Gen. 1985, 18, 1551. (36) Reiss, H.; Koper, G. J. M. J. Phys. Chem. 1995, 99, 7837. (37) Viisanen, Y.; Strey, R.; Reiss, H. J. Chem. Phys. 1993, 99, 4680. (38) Viisanen, Y.; Strey, R. J. Chem. Phys. 1994, 101, 7835. (39) McGraw, R.; Laaksonen, A. Phys. ReV. Lett. 1996, 76, 2754. (40) Laaksonen, A.; McGraw, R. Europhys. Lett. 1996, 35, 367. (41) McGraw, R.; Laaksonen, A. J. Chem. Phys. 1997, 106, 5284. (42) Talanquer, V. J. Chem. Phys. 1997, 106, 9957. (43) Koga, K.; Zeng, X. C. J. Chem. Phys. 1999, 110, 3466. (44) Girshick, S. L.; Chiu, C.-P. J. Chem. Phys. 1990, 93, 1273. (45) Hale, B. N. Met. Trans. 1992, 23A, 1863. (46) Weakliem, C. L.; Reiss, H. J. Phys. Chem. 1994, 98, 6408. (47) Onsager, L. Phys. ReV. 1931, 37, 405; 38, 2265. (48) Tolman, R. C. The Principles of Statistical Mechanics; Oxford University Press: London, England, 1938. (49) Temkin, D. E.; Shevelev, V. V. J. Cryst. Growth 1984, 66, 380. (50) Wilemski, G. Unpublished. (51) Courtney, W. G. J. Chem. Phys. 1961, 35, 2249. (52) Wilemski, G. J. Chem. Phys. 1975, 62, 3763. (53) Fits to the ethanol properties are taken from Wegener, P. P.; Clumpner, J. A.; Wu, B. J. C. Phys. Fluids 1972, 15, 1869. The hexanol vapor pressure was fitted using data from the CRC Handbook of Chemistry and Physics, 61st ed.; CRC Press: Boca Raton FL, 1980, supplemented by the value at 260 K from ref 13. The hexanol surface tension curve represents a linear fit between the value at 260 K from ref 13 and the critical temperature 586.7 K. The hexanol density is a linear fit between the value at 260 K from ref 13 and the value at 293.2 K from the CRC Handbook of Chemistry and Physics. (54) Zeng, X. C.; Oxtoby, D. W. J. Chem. Phys. 1991, 94, 4472. (55) Strey, R.; Viisanen, Y. J. Chem. Phys. 1993, 99, 4693. (56) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes: The Art of Scientific Computing; Cambridge University Press: Cambridge, England, 1985.
