Kinetics of particle coarsening with allowance for Ostwald ripening and coagulation

In the absence of particle coagulation and integration, one can list three main stages of the phase transition process in metastable media. They are (i) an initial stage during which critical nuclei of the new phase appear and evolve when the system metastability level may be regarded as roughly constant, (ii) an intermediate stage during which crystal growth concurrently with continuing nucleation of particles and metastability reduction are equally important, and (iii) a final relaxation stage of Ostwald ripening in which larger particles grow at the expense of the dissolution of smaller ones.

The first stage was studied by a number of authors [1, 2]. The evolution of a metastable system at the second stage was analytically described in a complete form very recently for single-component and binary systems [3, 4]. The final stage of phase transformations was analyzed for cases when one of the particle coarsening mechanisms (Ostwald ripening or coagulation) predominates [5–8]. However, these mechanisms cannot always be identified separately [9, 10]. In addition, it is important to study their combined effect for the following reasons. The metastability level decreases at the concluding stage of phase transitions and tends to zero, whereas the radius of critical nuclei increases. As a result, fine crystals vanish, and large crystals continue to grow. It was previously assumed that this process predominates and it was described by means of the Lifshitz–Slyozov (LS) theory [11, 12] in which it was shown that (i) the cube of the average particle radius is a linear function of time t, and (ii) the cube of the supersaturation and the number of particles per unit volume represent inverse functions of time in the limiting case $t\to \infty$. Although qualitative behavior of the LS asymptotic theory has been confirmed [13, 14], there are a lot of discrepancies between experimental data and theoretical and numerical predictions [15, 16]. It has been suggested that this situation can be connected with the effect of competition for space by the growing crystals [17]. This study showed that allowance for this competition greatly affects the size distribution function [18]. It is significant that both mechanisms, Ostwald ripening and coagulation, make a contribution of the same order of magnitude to the coarsening of particles [16, 19].

From the physical point of view the present study is interesting as it combines two mechanisms of coalescence and coagulation simultaneously for a long-duration self-similar stage of the phase transition process. From the mathematical point of view the present analysis is of some interest because we discuss here how to solve analytically a modified kinetic equation involving nonlinear terms which are responsible for these mechanisms. In view of the fact that a general theoretical method to solving the problem containing the kinetic integral equation with a memory kernel is absent, the most popular approach is to apply different numerical methods [20–23]. Another widespread approach is concerned with situations in which one of the aforementioned mechanisms predominates. To solve such problems, some researchers have been proposing to look for different approximate solutions of the corresponding kinetic and balance equations (e.g. the conventional methods of Laplace transforms and expansion of solutions with respect to a small parameter) [24–26]. Therefore, the question remains open, of how to develop the appropriate theoretical method of solution of this combined nonlinear problem involving the mechanisms of coalescence and coagulation. In this paper, a new analytical approach taking into account both mechanisms simultaneously is developed, seemingly for the first time. An important point is that this approach can be used for the description of numerous phenomena in fundamental and applied research, ranging from applications relating to the condensation and crystallization of liquids, to chemical industry, geophysics and atmospheric sciences.

Consider a late stage of phase transformation (condensation or crystallization), when a metastable system is only slightly supersaturated/supercooled and nucleation has come to a halt because the critical nucleus size is large enough. In addition, the critical nucleus volume v_* is also large and comparable with the average volume of the growing particles. Crystals with dimensions smaller than the critical volume v_* vanish, whereas particles with larger dimensions evolve. An important point is that the moving particles collide and agglomerate so that the process becomes more complicated than in the case of Ostwald ripening. Namely, two small particles which should normally vanish (e.g. dissolve or evaporate) are capable of merging to grow into a larger particle able, perhaps, to evolve.

The polydisperse ensemble of particles is described by their size density distribution function f(t, v) normed to the particle denumerable concentration (v and t are the particle volume and time). Let us consider the mean-field approximation in which the notion of space is not present. The kinetic equation of Ostwald ripening corrected for coagulation can be written as [10, 11, 19, 27]

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}\frac{\partial f}{\partial t}+\frac{\partial}{\partial v}(gf)+{{I}_{c}}=0,~{{I}_{c}}=\int_{0}^{\infty}A\left({{v}_{1}},v\right)f\left(t,{{v}_{1}}\right)f(t,v)\text{d}{{v}_{1}}\\ \quad -\frac{1}{2}\int_{0}^{v}A\left({{v}_{1}},v-{{v}_{1}}\right)f\left(t,{{v}_{1}}\right)f\left(t,v-{{v}_{1}}\right)\text{d}{{v}_{1}},~t>0,~v>0.\end{array}\end{eqnarray} \tag{ 1 }$$

The first two terms in equation (1) represent the standard continuity equation and the integral I_c is the coagulation term obtained by Lifshitz and Slyozov [11]. The factor 1/2 in this term is needed to eliminate the collisions between two groups of crystals twice (for more details, see [11, 28]). Here g(t, v) is the growth rate of a particle with volume v and A(v₁, v) is the frequency of paired collisions per unit volume. For the sake of simplicity, one can assume an equal coagulation rate for all the particles irrespective of their size, i.e. $A\left({{v}_{1}},v\right)=A=\text{const}$. In the case of supersaturated solutions, the growth rate can be written in the form [12]

$$\begin{eqnarray}&&g(t,v)=\frac{\text{d}v}{\text{d}t}=\beta \left[w(t)\frac{{{v}^{1/3}}}{\alpha}-1\right],~w(t)=\frac{\Delta (t)}{{{C}_{\infty}}},\end{eqnarray} \tag{ 2 }$$

where $\beta =8\pi \sigma D{{v}^{\prime}}{{C}_{\infty}}/(kT)$ is the kinetic parameter, $\Delta (t)=C(t)-{{C}_{\infty}}$ is the system supersaturation, and $\alpha =2{{\left(4\pi /3\right)}^{1/3}}\sigma {{v}^{\prime}}/(kT)$. Here C(t) and ${{C}_{\infty}}$ represent the gradually reducing volume concentration of the dissolved substance and its value at a planar interface, D and ${{v}^{\prime}}$ are the diffusion coefficient and the volume per atom of the dissolved component, σ is the surface tension, k is the Boltzmann constant, and T is the absolute temperature. In the case of supercooled melts the dimensionless supersaturation w should be replaced by the dimensionless supercooling [29] $w=\left({{\theta}_{\infty}}-\theta (t)\right)/{{\theta}_{\infty}}$, where θ and ${{\theta}_{\infty}}$ stand for the liquid temperature and its equilibrium value above the planar interface. As is easily seen from expression (2), the critical nucleus volume is ${{v}_{\ast}}={{\left(\alpha /w\right)}^{3}}$. This expression assumes that the critical nucleus is a sphere (because of the definition of α). In many cases, this approximation is not good (see, among others, [30]) and the theory under consideration should be modified in accordance with more precise expressions. However, in order to demonstrate the analytical method let us use here the aforementioned expressions.

Note that the kinetic equation (1) should be supplemented by the following mass or heat balance equations

$$\begin{eqnarray}&&C(t)-{{C}_{\infty}}+\int_{0}^{\infty}vf(t,v)\text{d}v={{Q}_{1}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\theta}_{\infty}}-\theta (t)+P\int_{0}^{\infty}vf(t,v)\text{d}v={{Q}_{2}},\end{eqnarray} \tag{ 4 }$$

where P  =  L/C_m, L and C_m are the latent heat of the phase transition and the specific heat of the mixture, and Q₁ and Q₂ stand for the initial saturation and supercooling in the absence of crystals.

The integro-differential problem (1)–(4) is extremely complex. Note that even in the absence of coagulation (integral) terms in equation (1), the problem under consideration can only be solved asymptotically in the limit of large times [11, 31, 32]. We consider a final stage of the phase transition taking into account the integral contributions in equation (1). In this case, the system supersaturation $C(t)-{{C}_{\infty}}$ (or supercooling ${{\theta}_{\infty}}-\theta (t)$) is much smaller than the initial supersaturation Q₁ (or initial supercooling Q₂), which makes it possible to neglect these terms in the balance equations (3) and (4). Also, note that the problem under consideration must be supplemented with the boundary condition at infinity of the form $f\left(t,\infty \right)=0$.

Let us introduce the following self-similar variables

$$\begin{eqnarray}&&w(t)=\frac{K}{{{\tau}^{1/3}}},~\tau =\frac{\beta t}{{{\alpha}^{3}}},~\varphi (x)=\frac{A{{\tau}^{2}}{{\alpha}^{6}}f(t,v)}{{{K}^{3}}\beta},~x=\frac{{{K}^{3}}v}{{{\alpha}^{3}}\tau},\end{eqnarray} \tag{ 5 }$$

where parameter K characterizes the self-similar solution. Note that it cannot be determined by considering the concluding self-similar stage under consideration, and thus must be regarded as experimentally specified (this parameter characterizes the rate of metastability reduction). Now substituting expressions (5) into equations (1)–(4), one can get

$$\begin{eqnarray}&&\frac{\text{d}\varphi}{\text{d}x}=\frac{1}{g(x)}\left[q(x)\varphi (x)+F(x)\right],~\int_{0}^{\infty}x\varphi (x)\text{d}x=R,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}m=\int_{0}^{\infty}\varphi \left(\,y\right)\text{d}y,~q(x)=\frac{{{K}^{3}}}{3{{x}^{2/3}}}-2+m,~\\ F(x)=-\frac{1}{2}\int_{0}^{x}\varphi \left(\,y\right)\varphi (x-y)\text{d}y,\end{array}\end{eqnarray}$$

where $g(x)=x-{{K}^{3}}\left({{x}^{1/3}}-1\right)$ and $R={{K}^{3}}A{{Q}_{1}}/\beta$ in the case of supersaturated solutions and $R={{K}^{3}}A{{Q}_{2}}/\left(\beta P\right)$ in the case of supercooled liquids.

Note that at $K\geqslant {{K}_{0}}={{(27/4)}^{1/3}}$ equation (6) has no solution. This is due to the fact that the integral expression for R diverges at some point x where g(x) becomes zero. Note that the integro-differential equation (6) has been solved numerically [20, 33] at K  <  K₀ on the basis of approximate estimates of the asymptotic behavior of $\varphi (x)$ at $x\to \infty$. The classical solution describing the ripening problem without coagulation exists merely at K  =  K₀ [11, 12].

Numerical solutions of the integro-differential system (6) show that the dimensionless distribution function φ is of an order of magnitude less than 10⁻¹ for all values of variable x. Estimating this function as $\varphi \lesssim {{10}^{-1}}$ one can find from equation (6) that $F(x)\ll q(x)\varphi (x)$. Therefore the zero-order approximation ${{\varphi}_{0}}(x)$ of the distribution function can be found from the following differential equation $\text{d}{{\varphi}_{0}}/\text{d}x=q(x){{\varphi}_{0}}(x)/g(x)$. Its integration gives

$$\begin{eqnarray}&&{{\varphi}_{0}}(x)=\varphi (0)\exp \left(\int_{0}^{x}\frac{q\left(\,y\right)\text{d}y}{g\left(\,y\right)}\right).\end{eqnarray} \tag{ 7 }$$

Substituting equation (7) into the integral for F, one can obtain its zero-order approximation

$$\begin{eqnarray}&&{{F}_{0}}(x)=-\frac{1}{2}\int_{0}^{x}{{\varphi}_{0}}\left(\,y\right){{\varphi}_{0}}(x-y)\text{d}y.\end{eqnarray} \tag{ 8 }$$

Now substitution of equation (8) into the right-hand side of the first equation (6) leads to

$$\begin{eqnarray}&&\frac{\text{d}{{\varphi}_{1}}}{\text{d}x}=\frac{1}{g(x)}\left[q(x){{\varphi}_{1}}(x)+{{F}_{0}}(x)\right].\end{eqnarray}$$

Integrating this equation one can find the first-order solution ${{\varphi}_{1}}$ for the distribution function. Then the aforementioned expression for F determines its first-order approximation F₁. By iterating this procedure over and over, we obtain the following complete solution

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{\varphi}_{i}}(x)=\left[\varphi (0)+\int_{0}^{x}\frac{{{F}_{i-1}}\left(\gamma \right)}{g\left(\gamma \right)}\exp \left(-\int_{0}^{\gamma}\frac{q\left(\,y\right)\text{d}y}{g\left(\,y\right)}\right)\text{d}\gamma \right]\\ \qquad \quad \times \exp \left(\int_{0}^{x}\frac{q\left(\,y\right)\text{d}y}{g\left(\,y\right)}\right),\end{array}\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&{{F}_{i-1}}(x)=-\frac{1}{2}\int_{0}^{x}{{\varphi}_{i-1}}\left(\,y\right){{\varphi}_{i-1}}(x-y)\text{d}y.\end{eqnarray}$$

and i  =  1, 2, 3,... Thus ${{\varphi}_{0}}(x)$ from (7) represents the main contribution whereas $\varphi (x)={{\varphi}_{i}}(x)$ from (9) is the ith approximation of the complete solution for the distribution function $\varphi (x)$.

Figure 1 shows the complete analytical solution of the problem under consideration (here we chose K and $\varphi (0)$ on the basis of estimations given by Alyab'eva and Mansurov [33]). As is seen, the zero-order solution ${{\varphi}_{0}}(x)$ represents the main contribution to the distribution function. The first- and second-order solutions (${{\varphi}_{1}}(x)$ and ${{\varphi}_{2}}(x)$) lie in close proximity to ${{\varphi}_{0}}$. In other words, the two next contributions ${{\varphi}_{1}}$ and ${{\varphi}_{2}}$ converge to the complete solution so that ${{\varphi}_{2}}$ does not change significantly the resulting distribution function. Note that the infinite series (9) quickly converges due to smallness of the integral term F in comparison with the product $q(x)\varphi (x)$ in the right-hand side of the first equation (6). Figure 2 demonstrates that the distribution function φ becomes broader and flatter and its maximum point moves to the right with increasing parameter K. This is due to the fact that parameter K is directly proportional to the metastability level w in accordance with the first equation (5). In other words, smallness of K implies smallness of w.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By using self-similar variables (5) one can derive the following expressions for the mean particle radius $\bar{a}(t)$ and the total number of particles N(t):

$$\begin{eqnarray}&&\bar{a}(t)=\frac{H(t)}{I}\int_{0}^{\infty}{{x}^{1/3}}\varphi (x)\text{d}x,~N(t)=\frac{I}{At},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&H(t)={{\left(\frac{3\beta t}{4\pi {{K}^{3}}}\right)}^{1/3}},~I=\int_{0}^{\infty}\varphi (x)\text{d}x,\end{eqnarray}$$

which demonstrate that $\bar{a}$ increases as t^1/3 and N decreases in inverse proportion to time. This behavior corresponds to the classical case of Ostwald ripening without coagulation [7, 12].

Figure 3 demonstrates the distribution function $W(u)=b{{u}^{2}}\varphi \left(x(u)\right)$ plotted according to the complete analytical solution (9) in comparison with experimental data [16] (here $u=a/\bar{a}$, $x(u)={{l}^{3}}{{u}^{3}}$, $l=K\bar{a}{{\left[4\pi /\left(3\beta t\right)\right]}^{1/3}}$, and a represents the particle radius). Notwithstanding the fact that Ostwald ripening predominates in these experiments (in comparison with particle coagulation), these data were chosen because the theory under consideration and experimental data [16] describe the concluding stage of phase transitions where the self-similar regime establishes. In addition, experimental data [16] were chosen due to the fact that (i) experimental conditions and theoretical considerations are very similar, and (ii) experimental curves of the distribution function enable us to find all combinations of unknown theoretical parameters.

Zoom In Zoom Out Reset image size

As is easy to see, the particle size distribution is in quantitative agreement with conducted experiments. The obtained distribution function taking into account the processes of coalescence and coagulation substantially differs from the LS asymptotic solution [11, 12] found in the absence of coagulation (solid and dashed lines in figure 3). Analytically this case is described by equation (1) in the absence of integral terms (I_c  =  0). An important point is that the dashed curve in figure 1 (the LS solution) exists only at $K={{K}_{0}}={{(27/4)}^{1/3}}$ and equals zero for x  >  27/8. The theory of coalescence combined with coagulation, contrastingly, has a solution at K  <  K₀ on the right of the blocking point x  =  27/8. In other words, the coagulation process substantially changes the Ostwald ripening stage of phase transitions. From the physical point of view this is explained by the competition of two evolutionary mechanisms (coalescence and coagulation) when coalescence plays an important role and particles approach each other due to Brownian motion (or another reason) and stick together. This situation completely differs from the Ostwald ripening mechanism when larger particles grow at the expense of smaller ones.

In summary, we have presented a theoretical investigation of particle coarsening kinetics with allowance for Ostwald ripening and coagulation at the concluding stage of phase transition processes. A complete analytical solution of the problem is found within the framework of the mean-field approximation when the particles interact only via a common mean-field, and their coagulation rate is constant. The obtained analytical solution of integro-differential equations demonstrates that the particle distribution function decreases and becomes broader and flatter when increasing the dimensionless volume x and highly depends on metastability reduction dynamics.

This work was supported by the Ministry of Education and Science of the Russian Federation under the project 315 and by the contract 02.A03.21.0006 (act 211 of Government of the Russian Federation).