6884
Langmuir 1997, 13, 6884-6888
Articles A Simple Computer Simulation of Ostwald Ripening Yves De Smet,* Luc Deriemaeker, and Robert Finsy Faculteit Wetenschappen, Fysische en Colloı¨d Chemie, Vrije Universiteit Brussel, Pleinlaan 2, B-1050 Brussel, Belgium Received April 14, 1997. In Final Form: August 20, 1997X The coarsening of a polydisperse set of N particles is studied by computer simulation. The main steps of the simulation are the following: generation of a set of N (typically 20 000) particles according to an a priori model of the particle size distribution, and subsequent transfers of molecules from one particle to another according to predetermined growth rules. As an application, we consider the Ostwald ripening process. In this case, the exchange of molecules between small and large particles is governed by Kelvin’s equation. The evolution of the average particle radius and of the size distribution is studied both in the stationary and nonstationary regimes. The typical characteristics of the stationary regime predicted by the Lifshitz-Slyozov-Wagner theory are recovered. For a more realistic initial log-normal size distribution, it was found that during an increase in average radius of about a factor of 4.3, the stationary regime was not yet attained. Results are also reported for the combined effects of Ostwald ripening and membrane resistance.
1. Introduction In our everyday life, we use many dispersed-phase systems (food emulsions, ink, paints, pharmaceutical and cosmetic preparations, photographic emulsions, etc.). A common feature of these products is that they are unstable. Their thermodynamic equilibrium state is a demixed one. This is ultimately attained by different processes whereby the dispersed-phase particles or droplets grow with time. The growth rate is of technological and fundamental interest. Many theoretical models and numerical studies have been presented for different growing, aging, or ripening mechanisms.1-17 In this paper, an extremely simple game that allows us to simulate on a common PC the coarsening of a polydisperse set of N particles is presented. The transport of molecular material from one particle to another is described by predetermined growth rules, allowing us to study empirically the effect of different kinds of hypotheses * Telephone: 32-2-629 36 89. Fax: 32-2-629 33 20. E-mail:
[email protected]. X Abstract published in Advance ACS Abstracts, November 15, 1997. (1) Tadros, T.; Vincent, B. In Encyclopedia of emulsion technology; Becher, P., Ed.; Marcel Dekker: New York, 1983; Vol. 1, p 129. (2) Friberg, S.; Yang, J. In Emulsions and emulsion stability; Sjo¨blom, J., Ed.; Marcel Dekker: New York, 1996; p 1. (3) Stein, H. The preparation of dispersions in liquids; Marcel Dekker: New York, 1996; p 15. (4) Kahlweit, M. In Physical Chemistry; Eyring, E., Henderson, D., Jost, W., Eds.; Academic: New York, 1970; Vol. 10, p 719. (5) Kabalnov, A. S.; Shchukin, E. D. Adv. Colloid Interf. Sci. 1992, 38, 69. (6) Lifshitz, I.; Slyozov, V. J. Phys. Chem. Solids 1961, 19, 35. (7) Wagner, C. Z. Elektrochem. 1961, 65, 581. (8) Marqusee, J. A.; Ross, J. J. Chem. Phys. 1984, 80 (1), 536. (9) Enomoto, Y.; Kawasaki, K. Acta Metallurg. 1989, 37 (5), 1399. (10) Yao, J. H.; Laradji, M. Phys. Rev. E 1993, 47 (4), 2695. (11) Yao, J. H.; Elder, K. R.; Guo, H.; Grant, M. Phys. Rev. B 1993, 47 (21), 14110. (12) Zheng, X.; Bigot, B. J. Phys. 2 Fr. 1994, 4, 743. (13) Pflu¨gl, W.; Titulaer, U. M. Phys. A 1995, 214, 52. (14) Abinandanan, T. A.; Johnson, W. C. Mater. Sci. Eng. 1995, B32, 169. (15) Maksimov, L. A.; Ryazanov, A. I.; Heinig, K.-H.; Reiss, S. Phys. Lett. A 1996, 213, 73. (16) Venzl, G. Ber. Bunsenges. Phys. Chem. 1983, 87, 318. (17) Voorhees, P. W. J. Stat. Phys. 1985, 38 (1/2), 231.
S0743-7463(97)00379-X CCC: $14.00
for the growth mechanism. As a case study, the effect of the increasing solubility of smaller particles as predicted by the Kelvin (or Gibbs-Thomson) equation is investigated. The stationary growth regime, for which analytical results4-7 are well-known, and the transition from an initial situation toward the stationary regime, for which no general analytical results are known,16,17 are studied. The combined effects of Ostwald ripening and membrane resistance are also investigated. 2. Simulation Game In this study, the coarsening of a set of N spherical particles is simulated as follows. First, an initial polydisperse set of N particles is generated according to an a priori model for the size distribution; i.e., a set of particle radii {ai,0}, where i ) 1, ..., N, is defined. Each particle i is filled up with ni,0 molecules with fixed molecular volume. Subsequently, the particles exchange molecules according to growth rules, which may be size-dependent. In each step j of the simulation, the number of molecules ni, j in a particle i is given by
ni, j ) ni, j-1 + MPi, j-1(ai, j-1)
(1)
In eq 1, M is proportional to the total number of exchanged molecules and Pi, j embodies the growth law. A particle will grow if Pi, j is positive and will shrink if Pi, j is negative. The critical radius ac, j of a particle that is neither growing nor shrinking at step j is given by Pi, j ) 0. The mass balance requires that at each step
∑i Pi, j ) 0
(2)
Equation 2 allows us, in principle, to determine ac, j if the explicit dependence of Pi, j with particle size is known. When the smallest particle (with index i ) N) contains fewer molecules than the number it should lose according to eq 1, the number M is replaced with Mj ) nN, j /M, particle N is discarded from the set, and the number of © 1997 American Chemical Society
A Simple Computer Simulation of Ostwald Ripening
Langmuir, Vol. 13, No. 26, 1997 6885
particles becomes N - 1. At each step of the simulation, the size distribution and the number-average particle radius 〈a〉j, as well as ac, j are computed. Typical simulations start with N ) 20 000 and end when N has dropped to 200. Finally, it is assumed that time is proportional to the k Mj, where number of exchanged molecules; i.e., tk ∼ ∑j)1 j and k are step numbers and Mj can be either M or nN, j /M (see above). Therefore, the dimensionless time k
Tk )
Mj ∑ j)1
(3)
is introduced.
ac, j )
(Ra)
C(a) ) C(∞) exp
(4)
where C(a) is the solubility of a particle with radius a and C(∞) is the bulk solubility; R is a material-dependent parameter (“capillary length”) typically of the order of 1 nm. A qualitative, theoretical description of the coarsening process driven by eq 4 was given by Lifshitz and Slyozov and by Wagner (LSW theory4-7). The main results are as follows. 1. There exists a critical radius ac. 2. A stationary regime is attained. At this point, the coarsening or ripening rate v is described by
d〈a〉 4 ) RDmC(∞) dt 9
(5)
where Dm is the diffusion coefficient of the dispersed-phase molecules in the continuous phase. 3. A limiting stationary distribution is obtained. An analytical form of this LSW distribution is given by 2
81eu exp[1/(2u/3 - 1)] 3 3
ac x32(u + 3)7/3(1.5 - u)11/3 (0 e u e 1.5) (u > 1.5)
(6)
with a cut-off for u ) 1.5. A particular feature of this LSW distribution is that 〈a〉 ) ac. Note that the LSW theory only describes the coarsening in the stationary regime, not the transition from an initial situation (freshly prepared dispersion or emulsion) toward the stationary regime. The growth law corresponding to Ostwald ripening is derived with the aid of Fick’s first law and Kelvin’s equation, assuming R , a, yielding
(
N
) 〈a〉j
M in eq 1 can be considered to be proportional to the product DmC(∞)R. The conversion of the number of exchanged molecules per particle into time units follows from eqs 1, 7, and 8, yielding ∆t for one step:
)
dn a -1 ) 4πDmC(∞)R dt ac
Hence, our simulation is carried out with
(7)
(9)
Hence, the dimensionless time T can be converted into time t by
t ) T/[4πDmC(∞)R] Note that the simulation with eqs 1 and 8 now allows us to describe both the stationary regime and the transition from an arbitrary initial situation toward a stationary regime, as will be shown in section 4.1. 3.2. Ostwald Ripening with Membrane Resistance. The Ostwald ripening theory assumes that the molecular diffusion through the continuous medium is the rate determining factor. In many practical systems (e.g., oil-in-water emulsions stabilized by adsorbed surfactant layers), the molecular transfer is also limited by the resistance to penetration through the adsorbed layer. For the growth rate of a particle surrounded by a thin membrane, Kabalnov and Shchukin5 propose the following equation:
(
)/( )
dn a -1 ) 4πDmC(∞)R dt ac
3
W(u) ) 0
∑i ai, j
∆t ) M/[4πDmC(∞)R]
3.1. Ostwald Ripening. A particular coarsening process is Ostwald ripening. Here the particle size mainly changes by solubilization and condensation of molecules. According to Kelvin’s equation, the larger particles will grow by condensation of material diffused through the continuous phase, coming from the smaller droplets that solubilize more than the large ones. Quantitatively, Kelvin’s equation reads
W(a/ac) ) W(u) )
(8)
From the mass balance (eq 2), it follows that
3. Applications
v)
ai, j -1 ac, j
Pi, j )
1+
β a
(10)
with
β)
δDm KDf
where δ is the membrane thickness, Dm and Df are the molecular diffusion coefficients in the continuous medium and the membrane, and K is the distribution coefficient between membrane and medium. When the membrane resistance becomes rate determining (i.e., β/a . 1), the growth process is altered as follows: (1) In the stationary regime, the square of the average radius, 〈a〉2 (instead of 〈a〉3), increases linearly with time.4 (2) Another broader limiting distribution is obtained:4
()
W
(
)
u 6 a ) W(u) ) 3 exp 5 ac 2 u a (2 - u) c
(0 e u e 2) W(u) ) 0 (u g 2)
(11)
For intermediary situations, i.e., both Ostwald ripening and membrane resistance (β/a ≈ 1), no analytical results for growth laws or limiting distributions can be derived with the LSW theory. Our simulation of this problem is carried out with
6886 Langmuir, Vol. 13, No. 26, 1997
Pi, j )
(
De Smet et al.
)(
)
ai, j β -1 / 1+ ac, j ai, j
(12)
The critical radius is obtained from the mass balance,
ac, j )
(
ai, j
∑i 1 + β/a
)/(
)
1
∑i 1 + β/a
i, j
i, j
Again M in eq 1 can be considered proportional to DmC(∞)R. 4. Results 4.1. Ostwald Ripening. The Ostwald ripening of sets of initially 20 000 particles was simulated using eqs 1 and 8. The following initial particle size distributions were considered. Firstly, the limiting stationary LSW distribution given by eq 6 with 〈a〉 ) 30 nm. In this case, it is expected that the shape of the size distribution remains unaltered during the simulation. Secondly, three initial Gaussian distributions
W(a) )
1
x2πσ
[
exp -
Figure 1. Initial size distributions: LSW (label 1), Gaussian with σ ) 4.5 nm (label 2), Gaussian with σ ) 9 nm (label 3), and log-normal distribution with σ ) 9 nm (label 4).
]
(a - 〈a〉)2 2σ2
with 〈a〉 ) 30 nm and with three different widths: σ ) 1.5, 4.5, and 9 nm. Finally, since many practical preparation procedures of dispersions and emulsions yield positively skewed size distributions, instead of the negatively skewed limiting LSW distribution, an initial log-normal distribution
W(a) )
1
x2πσ
[
exp -
]
[log(a) - 〈log(a)〉]2 2
2σ
with 〈a〉 ) 30 nm and σ ) log(9 nm) was also considered (for the log-normal distribution, σ is the standard deviation of log(a)). The different initial size distributions are plotted in Figure 1. The scale factor M was 104. Computations were performed with a PC with a Pentium processor. Typical computation times were 30 min. The final (N ) 200) average sizes 〈a〉 were about 130 nm, or about 4.3 times the initial ones. 4.1.1. Evolution of the Size Distributions. As expected, the scaled form of the distribution of the set generated with the LSW distribution appeared to be invariant during the entire simulation. In Figure 2, the evolution of the distribution of the set with the Gaussian distribution with σ ) 1.5 nm is shown. Clearly the distributions come closer and closer together and finally overlap each other. Note that the most important change in shape of the distribution is observed in the beginning of the simulation (〈a〉 ) 30-35 nm). Similar results are obtained with the initial Gaussian distribution with σ ) 4.5 nm. The broader Gaussian (σ ) 9 nm) and the log-normal (skewed with a tail toward the larger particle sizes) evolve slower toward a limiting form. This slower evolution is due to the initial presence of the larger particles, with sizes above the cut-off size of the LSW distribution. In Figure 3, the final size distributions (obtained from the five different initial ones) at 〈a〉 ) 130 nm are shown. The set generated with the log-normal distribution has not reached the stationary (LSW) distribution yet: it is still broader with more large particles than predicted by the stationary distribution. This illustrates that for more
Figure 2. Evolution of the Gaussian size distribution with σ ) 1.5 nm toward the stationary LSW distribution, in the case of Ostwald ripening. The narrow initial Gaussian (not shown in Figure 1) broadens very fast (even before the average particle size has increased with 1 nm) to the distribution labeled with 1 and evolves toward the stationary LSW distribution (labeled with 2).
Figure 3. Final distributions for the different sets in the case of Ostwald ripening, ordered from highest to lowest peak: the set with initially the LSW distribution (label 1); the Gaussian distributions (σ ) 1.5, 4.5, and 9 nm) (almost undistinguishable); the log-normal distribution (label 2).
practical situations, the initial growth process (an increase of a few times the initial average size) is different from the stationary regime predicted by the LSW theory. 4.1.2. Growth of the Average Particle Size. In Figure 4, 〈a〉3 is plotted as a function of the dimensionless time T. Clearly when starting with an LSW distribution 〈a〉3 increases linearly (the correlation coefficient is more than 0.9999). The slope appeared to be proportional to the
A Simple Computer Simulation of Ostwald Ripening
Figure 4. Evolution of 〈a〉3 with time T for the set with initially the LSW distribution (label 1), the Gaussians with σ ) 1.5 nm (label 2), 4.5 nm (label 3), and 9 nm (label 4), and the lognormal distribution (label 5).
Figure 5. Close-up of the early stage of the Ostwald ripening process. Same legend as in Figure 4.
parameter M, i.e., to the product DmC(∞)R or proportional to each one of the three factors in that product. In the first stage of the growth process (Figure 5), a different and nonlinear behavior is observed for the other initial distributions. The limiting slope, representing the Ostwald ripening rate, for the three Gaussian distributions corresponds to the one obtained with the LSW distribution. Such a limiting slope has not yet been obtained with the lognormal distribution. The conclusion is that the more the initial distribution contains larger particles than the stationary one, the slower the stationary regime obtained is. 4.2. Ostwald Ripening with Membrane Resistance. 4.2.1. Only Membrane Resistance. The ripening was simulated with sets of 20 000 particles. The same initial Gaussian and log-normal distributions as used in the previous section were investigated. As the reference initial distribution, the LSW distribution was replaced by the limiting stationary distribution of eq 11. This corresponds to the case when the membrane resistance becomes the rate-determining factor. In this instance, eq 1 was used with eq 12 for β/a . 1; i.e.,
Pi, j )
(
)
ai, j ai, j -1 β ac, j
The final average sizes were about 130 nm (4.3 times the initial ones). It appeared that now all distributions, except the log-normal, approach the same limiting distribution
Langmuir, Vol. 13, No. 26, 1997 6887
Figure 6. Final distributions for five sets of particles, initially with respectively a distribution corresponding to (ordered from highest to lowest peak) eq 11 (label 1), a Gaussian with σ ) 1.5, 4.5, and 9 nm (almost undistinguishable), and a log-normal distribution (label 2).
Figure 7. Evolution of 〈a〉2 with time T for the set with initially the distribution corresponding to eq 11 (label 1), the Gaussians with σ ) 1.5 nm (label 2), 4.5 nm (label 3), and 9 nm (label 4), and the log-normal distribution (label 5).
(see Figure 6), which is the one predicted by eq 11. Note that this limiting distribution is broader (cut-off at a/ac ) 2, instead of 1.5 for the LSW distribution). Therefore, the broader Gaussian distribution with σ ) 9 nm was able to evolve faster toward the stationary one compared to the previous case. The log-normal distribution, however, has not completely reached the stationary one: it remains a little skewed with a few larger particles. To verify the growth law, 〈a〉2 is now plotted as a function of dimensionless time. The results are presented in Figure 7. Again, when starting with the predicted stationary distribution (eq 11), a linear dependence is observed, and again in the first stage of the growth process, a different and nonlinear behavior is observed for the other initial distributions. Also, the limiting slopes for the three Gaussian distributions correspond to the one obtained with the stationary distribution. This limiting slope has not yet been reached by the set initially with the lognormal distribution. Again the conclusion is that the more large particles present in the initial distribution, compared to the stationary one, the slower the stationary regime obtained is. 4.2.2. Ostwald Ripening and Membrane Resistance. This case was simulated with sets of initially 20 000 particles and using eqs 1 and 12. The following initial size distributions were used: the distribution corresponding to eq 11, the three Gaussian distributions (σ ) 1.5, 4.5, and 9 nm), and the log-normal distribution.
6888 Langmuir, Vol. 13, No. 26, 1997
De Smet et al.
Figure 8. Final distributions of the set with initially the narrow Gaussian distribution (σ ) 1.5 nm) when simulated with β ) 0 (Ostwald ripening) (label 1), 102 (label 2), 103 (label 3), and 105 (only membrane resistance) (label 4).
Figure 9. Evolution of 〈a〉2 with time T. Same legend as in Figure 8.
In order to vary the relative contributions of Ostwald ripening and membrane resistance, the following values of β were used: 0 (pure Ostwald ripening), 10, 102, 103, 105 (only membrane resistance). In these cases (except for the first and the last values of β), the stationary size distribution cannot be derived analytically. However, one can expect them to be intermediate between the LSW distribution and the one given by eq 11. Since the evolution toward a stationary distribution is the fastest with a narrow Gaussian (σ ) 1.5 nm), we report in Figure 8 the limiting distributions of the sets generated with this Gaussian distribution for the different values of β. Clearly the expected trend from the smaller LSW distribution (β ) 0) toward the broader distribution of eq 11 (β ) 105) is observed. In Figure 9, the growth of 〈a〉3 with dimensionless time T is presented for the different values of β. Clearly nonlinear plots are obtained and the growth rate decreases with increasing membrane resistance, i.e., with increasing β values.
An extremely simple game that allows us to simulate the coarsening of dispersed particle systems is presented. As a first application, the Ostwald ripening process (both the diffusion-controlled and the membrane resistance determined cases) was studied. The familiar growth laws and limiting (and analytically predicted) forms of the particle size distributions were recovered when the game was started with initial size distributions identical or smaller than the limiting stationary ones. For broader initial distributions, especially for the more realistic log-normal distribution with a tail toward the larger particle sizes, the stationary regime was not completely obtained, even for a growth to about 4.3 times the initial average particle size. Furthermore, results are reported for the Ostwald ripening process controlled by both molecular diffusion and membrane resistance, for which no analytical results are available.
5. Conclusions
LA970379B
