Transient nucleation kinetics of crystal growth at the intermediate stage of bulk phase transitions

Transitions from metastability to thermodynamic stability occur when particles of the new phase grow on fluctuation-generated nuclei or on impurities serving as nucleation sites. Whereas at the early stages of the process it is possible to consider the individual particles as being independent of one another, as the phase transition progresses, the effect of particle interaction on the metastability level (supercooling or supersaturation) becomes substantial [1–3]. The dynamics of the evolving two-phase mixture at these stages is important from both theoretical and applied viewpoints. This applies particularly to crystallization from supercooled melts or supersaturated solutions, where quite frequently the bulk of growing nuclei is supplied by fluctuations, whereas the contribution of impurities as solidification sites is relatively small [4]. These processes control the size distribution in the products obtained in certain types of crystallizers and granulators [5–8]. The rate of formation of the solid particles has a decisive effect on the structure of the two-phase zone [9–11] and thus also on the physical and mechanical properties of the produced materials [12, 13].

The bulk crystallization processes from supercooled melts or supersaturated solutions can be theoretically divided into four basic stages [14]. A supercooled (supersaturated) state develops during the preliminary stage. The generation of new phase nuclei occurs during the initial stage caused by the homogeneous or heterogeneous nucleation. Then nucleation, crystal growth and desupercooling (desupersaturation) develop during the intermediate stage, whereas the late stage describes the relaxation processes such as Ostwald ripening and agglomeration.

As far as the initial stage is concerned, a relevant theory can be put forward by following common trends specific to the theory of nucleation in molecular systems [14]. Studies of this stage are usually based on the kinetic theory of nucleation and on the assumption of independence of individual nuclei [2, 15–17]. The last stage (recondensation of particles of the new phase) is analyzed on the basis of the Lifshitz–Slyozov model, which leans markedly on the assumption that the appearance of new nuclei can be neglected [1, 3, 18–21].

This study is concerned with the intermediate stage, when the effect of the growing particles on the state of the ambient medium and the ongoing nucleation rate are of equal importance. A reliable theory of this stage is incomplete, although it comprises a major part of the whole phase transition process. Such a regrettable state of art appears undoubtedly to be due to insurmountable mathematical difficulties which have presented a challenge until quite recently. Let us especially note that the pertinent information on the metastability reduction and about the evolution of new phase nuclei can be obtained through an analysis of just this stage. Those difficulties connected with unknown methods for solving integro-differential equations represent a probable reason why some growth models entirely ignore changes in metastability [22, 23]. The necessity of allowing simultaneously for the nucleation and the growth of new phase elements has been acknowledged in [24, 25]. Nevertheless, some models developed recently rely, in fact, on a concept of the steady-state approximations [26, 27], what may actually be true merely at the initial stage or in the case of certain process parameters.

The first attempt to overcome the aforementioned mathematical difficulties has been made in [28] on the basis of the saddle-point method for the Laplace integral with memory kernel. Next to that, a similar approach has been applied to the problems of evolution of aggregates in terrestrial magma oceans [29] and metastable colloids and magnetic fluids [30, 31]. The solutions and approaches developed in [28] have been used for the description of protein crystallization [32] and for the verification of the von Weimarn rules governing the average size of crystals [33]. It is important to note that the analytical solutions of [28] have been found only for a principal (main) contribution to the Laplace integral. This contribution corresponds to the zero-order approximation of the Laplace integral in the vicinity of a saddle point when the dimensionless Gibbs number p entering in the nucleation rate I = I_*exp ( − p/w²) formally tends to infinity (here I_* and w stand for the pre-exponential factor and dimensionless supercooling/supersaturation). Since p ranges from ∼10⁻¹ to ∼10³ for actual supercooled melts and supersaturated solutions this theory should be generalized to a broad range of metastable systems. It is significant that even in the case of sufficiently large p ∼ 10²–10³ the zero-order approximation (one-term expansion) obtained in [28] significantly differs from a complete solution containing next terms of infinite series (section 3).

This paper is organized as follows. In section 2 we briefly describe the rates of nucleation and particle growth as well as the integro-differential model of the intermediate stage of phase transitions, in section 3 and in the appendices we describe a complete analytical solution of the model under consideration and discuss some theoretical results, and conclusions are presented in section 4.

To be able to calculate coupled nucleation and coarsening of nuclei in the supercooled or supersaturated system one first of all must know the nucleation rate, which is the main parameter of the kinetics of decomposition of the metastable state in dispersed systems [34]. The molecular-kinetic theory of nucleation and methods of calculation of nucleation rates for single-component systems were established by many authors (see, among others, [15, 16, 35–37]). This theory is based on the assumption that micronuclei of the new phase can be formed in metastable melts or solutions as a result of density fluctuations. Then these nuclei become able to grow further when their characteristic dimension exceeds some critical value of r_*.

It is known from thermodynamics that the inception of a spherical particle of a new phase with radius r in an originally single-component system requires performing work, which is equal to

$$\begin{equation} W(r) = 4\pi r^2 \gamma _i - \case{4}{3}\pi r^3\rho _s(\mu _l -\mu _s) , \end{equation} \tag{ 1 }$$

where γ_i is the surface tension, ρ_s is the density of the stable solid phase and μ_l and μ_s are the chemical potentials of the original liquid and stable solid phases. For the metastable liquid system μ_l > μ_s, for which reason W(r) attains a maximum corresponding to some critical particle size r_*. Since W(r) is equal to the change in the Gibbs thermodynamic potential, which is at minimum in the stable state, particles with radii r < r_* are unstable and vanish, whereas those with r > r_* are unstable and grow. Nuclei with r = r_* are in a state of unstable equilibrium with the matrix medium. The work of their formation W_* follows from equation (1) at r = r_*

$$\begin{equation} W_* = \frac{16\pi \gamma _i^3}{3\rho _s^2(\mu _l -\mu _s)^2}. \end{equation} \tag{ 2 }$$

The difference between chemical potentials can be calculated from the Gibbs–Helmholtz equation and is equal to [38]

$$\begin{equation} \mu _l -\mu _s = \frac{L\Delta \theta }{\theta _p},\ \Delta \theta = \theta _p - \theta _l, \end{equation} \tag{ 3 }$$

where L is the latent heat of phase transition per unit mass of substance and θ_l and θ_p are the temperatures of the melt and the phase transition. The above indicates that the inception of nuclei of a new phase in a metastable system can be treated as the passage through an energy barrier of height W_* that interferes with nucleation. In this case the nucleation rate can be represented as an exponential function of the height of the energy barrier (see, among others, [14, 37, 39]). Hence, using equations (2) and (3), we have

$$\begin{equation} I=I_*\exp \left(\frac{-W_*}{k_B\theta _l} \right) =I_*\exp \left(\frac{-16\pi \gamma _i^3\theta _p^2}{3\rho _s^2L^2\Delta \theta ^2k_B\theta _l} \right), \end{equation} \tag{ 4 }$$

where k_B is the Boltzmann constant. Note that the pre-exponential factor I_* depends little on the supercooling/supersaturation [16, 35, 36, 40]. This factor can be a function of r [26] but in the theory under consideration where the nucleation rate enters into the boundary condition, I_* can be regarded as constant. In addition, since Δθ ≪ θ_p, equation (4) can be rewritten as

$$\begin{equation} I=I_*\exp \left(\frac{-16\pi \gamma _i^3\theta _p}{3\rho _s^2L^2\Delta \theta ^2k_B} \right). \end{equation} \tag{ 5 }$$

The aforementioned model is valid for single-component melts and can be modified in the case of binary melts, since the change in the free energy stemming from nucleation depends on the composition of the melt. However, as was shown in [12, 41], the surface diffusion of impurity of the growing nuclei can be neglected, and the expression for the nucleation rate can be written similarly to equation (5). In this case, the phase transition temperature θ_p entering in Δθ depends on the concentration C_l of impurity and thus θ_p should be replaced by a function θ_p(C_l).

For convenience in subsequent presentation we shall introduce some characteristic supercooling Δθ₀ and shall write equation (5) in the form

$$\begin{equation} I=I_*\exp \left[\frac{-p}{(\Delta \theta /\Delta \theta _0)^2} \right] ,\quad p=\frac{16\pi \gamma _i^3\theta _p}{3\rho _s^2L^2\Delta \theta _0^2k_B}, \end{equation} \tag{ 6 }$$

where p has the meaning of the dimensionless Gibbs number, corresponding to supercooling Δθ₀ and equation (6) expresses the nucleation rate as a function of the relative supercooling w = Δθ/Δθ₀.

Note that in the case of supersaturated solutions the expression analogous to equation (6) is [42]

$$\begin{equation} I=I_*\exp \left[\frac{-p}{\ln ^2 (C_l /C_p)} \right] ,\quad p=\frac{16\pi \gamma _i^3M_s^2}{3\rho _s^2R_g^2 \theta _s^3k_B}, \end{equation} \tag{ 7 }$$

where C_p is the concentration at saturation, R_g is the universal gas constant, M_s is the molecular weight and θ_s is the temperature of solution.

The growth of a particle appearing as a result of nucleation involves surface processes, the removal of heat from the phase interface or the supply of crystallizing component to this interface. This problem represents a Stefan-type heat or mass transfer problem with a moving boundary of the phase transition. Since there is no exact analytical solution of this problem, we use its traditional steady-state approximation. In this case the temperature (concentration) field θ around growing particle and its growth rate dr/dτ are determined from the following boundary value problem

$$\begin{equation*} \nabla ^2\theta =0,\quad \rho > r(\tau ), \end{equation*}$$

$$\begin{equation} \frac{{\rm d}r}{{\rm d}\tau }=-\frac{\lambda _l}{\rho _sL}\frac{\partial \theta }{\partial \rho }=\beta _* (\theta _p-\theta ),\quad \rho = r(\tau ), \end{equation} \tag{ 8 }$$

$$\begin{equation*} \theta \rightarrow \theta _l,\quad \rho /r(\tau ) \gg 1. \end{equation*}$$

Here ρ is the spherical coordinate, τ is the growth time, λ_l is the temperature conductivity coefficient, and β_* is the kinetic parameter.

The solution of equations (8) determines the rate of particle growth in the form

$$\begin{equation} \frac{{\rm d}r}{{\rm d}\tau }=\frac{\beta _*\Delta \theta }{1+\beta _* qr},\quad q=\frac{\rho _s L}{\lambda _l}. \end{equation} \tag{ 9 }$$

Equation (9) demonstrates that the growth rate is a linear function of the supercooling. For very fine crystals (r ≪ (β_*q)⁻¹) the rate of growth does not depend on their size. This mode of particle growth is known as 'kinetic', since it is fully defined by surface processes. The rate of growth for crystals whose dimensions exceed markedly the value (β_*q)⁻¹ is controlled by the rate of heat removal. In the case of supersaturated solutions expression (9) transforms as

$$\begin{equation} \frac{{\rm d}r}{{\rm d}\tau }=\frac{\beta _*\Delta C}{1+\beta _* qr},\quad q=\frac{C_p}{D}, \end{equation} \tag{ 10 }$$

where ΔC and D are the supersaturation and diffusion coefficient.

Note that there exist the growth theories in which expressions different from equations (9) and (10) are used (see, among others, [4, 39]). So, for example, when necessary, the dependence of phase transition temperature on the surface curvature should be taken into account. The aforementioned expressions for the nucleation rate and for the rate of particle growth have been checked against experimental data and in many cases were found to adequately describe the corresponding growth process [4, 13, 39, 40, 43]. The use of these expressions makes it possible to study the interaction between nucleation and particle coarsening in metastable media. Note that the theory developed below is easily adjusted for use with a large variety of expressions for the growth rate.

Let us consider a suspension of solid spherical particles in a macroscopically homogeneous single-component supercooled (supersaturated) system. The evolution of such a two-phase medium is described by a Fokker–Plank type kinetic equation for density f(τ, r) of particle radius distribution function and by the heat balance for the supercooling Δθ. Neglecting random fluctuations in the rate of particle growth, and assuming that the temperature distribution is uniform, we have

$$\begin{equation} \frac{\partial f}{\partial \tau }+\frac{\partial }{\partial r}\left( \frac{{\rm d}r}{{\rm d}\tau }f\right) =0,\quad r>r_*, \end{equation} \tag{ 11 }$$

$$\begin{equation} \rho _mC_m\frac{{\rm d}\Delta \theta }{{\rm d}\tau }=-4\pi \rho _s L\int_{r_*}^{\infty }r^2\frac{{\rm d}r}{{\rm d}\tau }f {\,\rm d}r, \end{equation} \tag{ 12 }$$

where ρ_m and C_m are the density and specific heat of the mixture. We assume for simplicity that ρ_m and C_m are independent of the temperature and also of the volume fraction of the solid phase. We assume that the liquid, which previously did not contain any crystals, is virtually instantaneously cooled at zero time to below the crystallization temperature by the amount Δθ₀. Then the boundary conditions that should be imposed onto solutions of equations (11) and (12) are

$$\begin{equation} \, f=0,\quad \Delta \theta = \Delta \theta _0,\quad \tau =0, \end{equation} \tag{ 13 }$$

$$\begin{equation} \frac{{\rm d}r}{{\rm d}\tau }f=I(\Delta \theta ),\quad r=r_*. \end{equation} \tag{ 14 }$$

The last condition defines the flux of nuclei that have overcome the critical barrier. Expressions (4)–(7) and (9)–(14) comprise the total system of equations for describing the intermediate stage of bulk phase transitions.

Designating I₀ = I(Δθ₀), we introduce the following dimensionless variables and parameters

$$\begin{equation*} t=\frac{\tau }{\tau _0},\quad s=\frac{r}{l_0},\quad w=\frac{\Delta \theta }{\Delta \theta _0}, \ F=l_0^4f, \end{equation*}$$

$$\begin{equation} \tau _0=\big( \beta _*^3\Delta \theta _0^3I_0\big)^{-1/4}, \ l_0=( \beta _*\Delta \theta _0)^{1/4}I_0^{-1/4}, \end{equation} \tag{ 15 }$$

$$\begin{equation*} \alpha _*=\frac{\beta _*\rho _sLl_0}{\lambda _l},\quad b_1=\frac{4\pi \rho _s L}{\rho _mC_m\Delta \theta _0}. \end{equation*}$$

In terms of variables (15) the problem under consideration have the form (r_* → 0)

$$\begin{equation} \frac{\partial F}{\partial t}+w\frac{\partial }{\partial s}\left( \frac{F}{1+\alpha _* s}\right) =0,\quad s>0, \end{equation} \tag{ 16 }$$

$$\begin{equation} \frac{{\rm d}w}{{\rm d}t}=-b_1w\int _{0}^{\infty }\frac{s^2 F(t,s)}{1+\alpha _*s}{\rm d}s, \end{equation} \tag{ 17 }$$

$$\begin{equation} F=0,\quad w=1,\quad t=0, \end{equation} \tag{ 18 }$$

$$\begin{equation} F=\frac{1}{w}\exp[ pg(w)] ,\quad s=0. \end{equation} \tag{ 19 }$$

Note that the dimensionless model (16)–(19) is valid for particle growth from supercooled melts and supersaturated solutions. The function g(w) entering in equation (19) defines the law of crystallization kinetics.

In the case of supercooled melts (subscript 'sm') with the Weber–Volmer–Frenkel–Zel'dovich (WVFZ) kinetics, we have from equation (6)

$$\begin{equation} g(w)=g_{{\rm sm}}=1-w^{-2}. \end{equation} \tag{ 20 }$$

In the case of an empirical law for the nucleation rate I = I_*(Δθ)^p which is frequently used in analyzing many industrial processes in various crystallizers [44], we get

$$\begin{equation} g(w)=g_{{\rm sm}}=\ln w. \end{equation} \tag{ 21 }$$

We will point out that in the case of supersaturated solutions (subscript 'ss') w represents the relative supersaturation, whereas the dimensionless parameters are

$$\begin{eqnarray} &&w=\frac{\Delta C}{\Delta C_0},\quad \alpha _*=\frac{\beta _*C_pl_0}{D},\quad b_1=\frac{4\pi C_p}{\Delta C_0},\nonumber\\ &&\tau _0=\big( \beta _*^3\Delta C_0^3I_0\big)^{-1/4}, \quad l_0=( \beta _*\Delta C_0)^{1/4}I_0^{-1/4}, \end{eqnarray} \tag{ 22 }$$

and ΔC₀ is the initial value of supersaturation.

For supersaturated solutions with the WVFZ kinetics [14], we have from equation (7)

$$\begin{equation} g(w)=g_{{\rm ss}}=\ln ^{-2}\left( 1+\frac{1}{w_p}\right) -\ln ^{-2}\left( 1+\frac{w}{w_p}\right) , \end{equation} \tag{ 23 }$$

where w_p = C_p/ΔC₀.

In the case of supersaturated solutions with the Mier kinetics I = I_*(ΔC)^p [39], we obtain

$$\begin{equation} g(w)=g_{{\rm ss}}=\ln w. \end{equation} \tag{ 24 }$$

Note that the boundary condition (19) can be used for homogeneous and heterogeneous nucleation [45].

The kinetics of metastability reduction is governed by the strongly non-linear integro-differential equations (16) and (17). We demonstrate below the analytical method of complete solution of this model based on the saddle-point technique for the Laplace type integral [46].

The exact solution of equation (16) with the initial and boundary conditions (18) and (19) has the form

$$\begin{equation} F(t,s)=( 1+\alpha _* s) \varphi( x(t) - y(s)) \eta ( x(t) - y(s)) , \end{equation} \tag{ 25 }$$

where we use the notations

$$\begin{equation*} x(t)=\int _{0}^{t}w(t){\,\rm d}t,\quad y(s)=\int _{0}^{s}(1+\alpha _* s){\,\rm d}s =s+\frac{\alpha _* s^2}{2}, \end{equation*}$$

$$\begin{equation*} \varphi (x(t))=\frac{1}{w(t)}\exp [ pg(w(t))] , \end{equation*}$$

and η(x) designates the Heaviside function.

Integrating the dimensionless growth rate

$$\begin{equation*} \frac{{\rm d}s}{{\rm d}t}=\frac{w(t)}{1+\alpha _*s}, \end{equation*}$$

which is obtained from equation (9) or (10) in the case of supercooled melts or supersaturated solutions when respectively using variables (15) or (22), with the initial condition that is s = 0 at t = ν, we come to

$$\begin{equation} \alpha _*s =\sqrt{1+2\alpha _*[ x(t)-x(\nu )]}-1 . \end{equation} \tag{ 26 }$$

Expression (26) defines, at ν = 0, the radius s_m(t) of the crystals that appeared at zero time, as of time t, in dimensionless form.

Now, replacing the variable s at any constant t by the new variable ν

$$\begin{equation*} x(\nu ) =x(t)-y(s), \end{equation*}$$

we see that the limits s = 0 and s = s_m(t) of integration with respect to s have corresponding to them limits ν = t and ν = 0 of integration with respect to ν. Thus, substituting expressions (25) and (26) into the integral equation (17), we get

$$\begin{equation} \frac{{\rm d}w}{{\rm d}t}=-b_1w\int _{0}^{t} h(\nu ,t)\exp [ pg(\nu )] {\,\rm d}\nu , \end{equation} \tag{ 27 }$$

where

$$\begin{equation} h(\nu ,t)=\frac{\lbrace \sqrt{1+2\alpha _*[ x(t)-x(\nu )]} -1\rbrace ^2}{\alpha _*^2\sqrt{1+2\alpha _*[ x(t)-x(\nu )]}}, \end{equation} \tag{ 28 }$$

and g(ν) = g(w(ν)) is determined by equations (20), (21), (23) or (24).

We further simplify equation (27) by making use of the fact that usually the dimensionless Gibbs number p is much greater than unity. This corresponds to rather high nucleation activation energy. Therefore it is possible to evaluate equation (27) by the familiar saddle-point method for the Laplace integral.

It follows from equations (20), (21), (23) and (24) that dg/dν < 0 for all kinetic mechanisms under consideration (so, for example, in the case of the WVFZ and Mier kinetics for supercooled melts, we have dg/dν = 2w⁻³dw/dν < 0 and dg/dν = w⁻¹dw/dν < 0, respectively). This clearly demonstrates that function g(ν) attains its maximum at boundary ν = 0 of the interval of its definition. Calculating the derivatives of w with respect to ν by the same equation (27), we see that the first three derivatives of g with respect to ν at point ν = 0 become zero for all aforementioned kinetic mechanisms, whereas the nth derivative g⁽ⁿ⁾|_{ν = 0} ≠ 0 at n ⩾ 4. The first four (nonzero) of them are listed in table 1, where

$$\begin{equation*} \chi = \left( \frac{6}{b_1}\right)^{1/4},\quad \kappa = \frac{1}{( 1+w_p)\ln ^3\big( 1+w_p^{-1}\big)}. \end{equation*}$$

Thus, the saddle-point method for the Laplace integral (27) gives the following complete solution [46]

$$\begin{equation} \int _{0}^{t} h(\nu ,t)\exp [ pg(\nu )] {\,\rm d}\nu \approx \sum \limits _{k=0}^{\infty }p^{-(k+1)/4}a_k(x(t)), \end{equation} \tag{ 29 }$$

$$\begin{equation} a_k(x(t))=(-1)^{k+1}\frac{4^k}{k!}\Gamma \left( \frac{k+1}{4}\right)H^k(0) \frac{{\rm d}^k}{{\rm d}\nu ^k}[ h(\nu ,t)H(\nu )]_{\nu =0}, \end{equation} \tag{ 30 }$$

where Γ is the Euler's gamma function and

$$\begin{equation*} H(\nu )=\frac{[ g(\nu )]^{3/4}}{g^\prime (\nu )}. \end{equation*}$$

An exact form of the first four coefficients a_k used for calculations is given in appendix A.

Table 1. Parameters and functions of the analytical solution for different kinetic mechanisms. Here z = z(x(t)) = 1 + 2α_*x(t) and f₀, f₁, f₂ and f₃ are given in appendix A.

	Supercooled melts		Supersaturated solutions
	WVFZ kinetics $I=I_*\exp \left[ \frac{-p}{\left(\Delta \theta /\Delta \theta _0\right)^2}\right]$	Mier kinetics I = I*(Δθ)p	WVFZ kinetics $I=I_*\exp \left[ \frac{-p}{\ln ^2 \left( C_l/C_p\right)}\right]$	Mier kinetics I = I*(ΔC)p
Variables	(6) and (15)	(15)	(7) and (22)	(22)
g(w)	equation (20)	equation (21)	equation (23)	equation (24)
g(IV)|ν = 0	−4b1	−2b1	−4κb1	−2b1
g(V)|ν = 0	24b1α*	12b1α*	24κb1α*	12b1α*
g(VI)|ν = 0	$-180b_1\alpha _*^2$	$-90b_1\alpha _*^2$	$-180\kappa b_1\alpha _*^2$	$-90b_1\alpha _*^2$
g(VII)|ν = 0	$1680b_1\alpha _*^3$	$840b_1\alpha _*^3$	$1680\kappa b_1\alpha _*^3$	$840b_1\alpha _*^3$
H|ν = 0	−χ/4	−21/4χ/4	−κ−1/4χ/4	−21/4χ/4
H'|ν = 0	−3α*χ/20	−3α*21/4χ/20	−3κ−1/4α*χ/20	−3α*21/4χ/20
H''|ν = 0	$9\alpha _*^2\chi /50$	$9\alpha _*^22^{1/4}\chi /50$	$9\kappa ^{-1/4}\alpha _*^2\chi /50$	$9\alpha _*^22^{1/4}\chi /50$
H'''|ν = 0	$-57\alpha _*^3\chi /100$	$-57\alpha _*^32^{1/4}\chi /100$	$-57\kappa ^{-1/4}\alpha _*^3\chi /100$	$-57\alpha _*^32^{1/4}\chi /100$
w0(z)	1 − f0(z)	1 − 21/4f0(z)	1 − κ−1/4f0(z)	1 − 21/4f0(z)
w1(z)	w0(z) + f1(z)	$w_0(z)+\sqrt{2}f_1(z)$	w0(z) + κ−1/2f1(z)	$w_0(z)+\sqrt{2}f_1(z)$
w2(z)	w1(z) + f2(z)	w1(z) + 23/4f2(z)	w1(z) + κ−3/4f2(z)	w1(z) + 23/4f2(z)
w3(z)	w2(z) + f3(z)	w2(z) + 2f3(z)	w2(z) + κ−1f3(z)	w2(z) + 2f3(z)

Now substituting

$$\begin{equation*} \frac{{\rm d}w}{{\rm d}t}=\frac{{\rm d}w}{{\rm d}x}\frac{{\rm d}x}{{\rm d}t}=w\frac{{\rm d}w}{{\rm d}x} \end{equation*}$$

into the left-hand side of equation (27) and taking into account expression (29), we come to the differential equation for dimensionless supercooling/supersaturation

$$\begin{equation*} \frac{{\rm d}w}{{\rm d}x}=-b_1\sum \limits _{k=0}^{\infty }p^{-(k+1)/4}a_k(x) . \end{equation*}$$

Its integration gives the complete solution for w

$$\begin{equation} w(x)= -b_1p^{-1/4}\int _{0}^{x}\sum_{k=0}^{\infty }p^{-k/4}a_k(x){\,\rm d}x +1. \end{equation} \tag{ 31 }$$

To be able to define w as a function of time t one must find the expression for t(x). Substituting w(x) from expression (31) into w = dx/dt, we finally obtain

$$\begin{equation} t(x)=\int _0^x \frac{{\rm d}x}{w(x)}. \end{equation} \tag{ 32 }$$

Thus, expressions (31) and (32) determine the relative supercooling/supersaturation w as a function of t in a parametric form. This function is demonstrated in figure 1 for two kinetic mechanisms of particle growth from supercooled melts (physical parameters used for calculations are typical for supercooled melts [47–50]: γ_i = 0.3 J m⁻², θ_p = 1273 K, L_V = 7 × 10⁹ J m⁻³, ρ_l = ρ_m = 7 × 10³ kg m⁻³, ρ_s = 7.8 × 10³ kg m⁻³, λ_l = 63 J m⁻¹ K⁻¹ s⁻¹, β_* = 0.5 × 10⁻⁷ m s⁻¹ K⁻¹, C_m = 840 J kg⁻¹ K⁻¹, I₀ = 10¹¹ m⁻³ s⁻¹). Here w = w₀ represents the main contribution to the sum in equation (31), whereas w = w_i (i = 1, 2, 3) is the ith approximation of the complete solution for w (these functions are written out in table 1 and appendix A), i.e.

$$\begin{equation*} w=w_i= -b_1p^{-1/4}\int _{0}^{x}\sum \limits _{k=0}^{i}p^{-k/4}a_k(x){\,\rm d}x +1. \end{equation*}$$

Zoom In Zoom Out Reset image size

Our calculations show that the main contribution w₀ (the limiting case previously analyzed by Buyevich and Mansurov [28]) is in rather poor agreement with the complete solution of equation (31). As is seen from the solution, the resulting supercooling/supersaturation (w = w₀ and w = w₃) may vary several-fold. In other words, the next summands of the oscillating series (31) make an important contribution to the complete solution for w. Figure 1 demonstrates that two next contributions w₁ and w₂ with alternating signs converge to the solution so that w₃ does not significantly change the resulting supercooling. Our calculations show that this conclusion is also valid for sufficiently small values of p ∼ 10⁻¹. In this case, the infinite series (31) quickly converges due to small values of α_* (which is of the order of 10⁻³–10⁻⁴ for actual metastable melts) at any fixed x and t.

Figure 1 illustrates that the complete solution for w contains the information about the rate of metastability reduction. If the liquid is initially supercooled to a greater degree (figures 1(a) and (b)) then its desupercooling rate is greater as well (mathematically this is explained by an inflexion point in figure 1(b) for small t). On the other hand, w₂ and w₃ in figures 1(b) and (c) show that the initially rapid metastability reduction becomes lower (plateau in figure 1(c)) after a lapse of time when the nucleation decelerates and newly born particles continue to grow.

The above results make it possible to find an explicit expression for the crystal-size density distribution function directly from analytical solutions (25) and (31). Taking into account that the maximum size of crystals that appeared at zero time is

$$\begin{equation} s_m(x(t))=\frac{\sqrt{1+2\alpha _*x(t)}-1}{\alpha _*}, \end{equation} \tag{ 33 }$$

we have

$$\begin{equation} F(x(t),s)=( 1+\alpha _* s) E_1^{-1}(t,s)\exp \left[ -\frac{p(1-E_1^2)}{E_1^2}\right] , \end{equation} \tag{ 34 }$$

where

$$\begin{equation*} E_1(x(t),s)=1-b_1\int _{0}^{x(t)-s-\alpha _*s^2/2}\sum \limits _{k=0}^{\infty }p^{-(k+1)/4}a_k(x)\,{\rm d}x. \end{equation*}$$

Here s is bounded by s_m and x(t) is determined by equation (32) in the form of its inverse function. Note that expression (33) follows from equation (26) at ν = 0 or from setting equal to zero the argument x(t) − y(s) of the Heaviside function in equation (25).

Taking into account the first four coefficients a_k(x) calculated in appendix A, we demonstrate the time evolution of dimensionless density distribution function F in figure 2. As is easy to see, weak inflection points appear with increasing the growth time. This is due to the fact that the high rate nucleation occurs at the initial stage of the process under consideration when the relative supercooling/supersaturation w = 1. Then the nucleation process will be accompanied by the particle growth reducing the level of supercooling/supersaturation in the system. This suppresses the nucleation mechanism and leads to the appearance of inflection points in the density distribution function. Note that such distributions for different materials have been previously reported in [51–54].

Zoom In Zoom Out Reset image size

An important parameter of the intermediate stage of phase transitions is the total number of nuclei (number of particles per unit volume) in the metastable system

$$\begin{equation} N(x(t))=\frac{1}{l_0^3}\int _0^{s_m(x(t))}F(x(t),s)\,{\rm d}s, \end{equation} \tag{ 35 }$$

where x(t) and s_m(x) are determined by equations (32) and (33). Another important parameter of the disperse system is the mean particle size which is

$$\begin{equation} {\bar{L}}(x(t))=l_0\int _0^{s_m(x(t))}sF{\rm d}s\left( \int _0^{s_m(x(t))}F{\rm d}s\right)^{-1}. \end{equation} \tag{ 36 }$$

Expressions (35) and (36) represent the complete solutions in a parametric form.

The relative (dimensionless) number of nuclei n_i(t) and mean particle size u_i(t) corresponding to the ith order solution in series of the distribution function (34) are demonstrated in figure 3 (n_i and u_i are determined in appendix B). One can readily see that the solution quickly converges so that the second-order approximations n₂ and u₂ practically coincide with the third-order solutions n₃ and u₃. Of fundamental importance is the stabilization of these system characteristics at the intermediate stage with increasing the process time. Let us especially emphasize that the main contributions n₀ and u₀ (zero-order solutions) previously obtained in [28] represent a rough approximation of the complete solution for actual metastable systems.

Zoom In Zoom Out Reset image size

In this paper, a complete solution of an integro-differential model describing nucleation and growth of particles at the intermediate stage of phase transitions is found. The initial non-linear system of governing equations is transformed to an exact solution of the Fokker–Plank type equation for the density distribution function F and to the integro-differential equation for the relative supercooling/supersaturation in the case of different kinetic mechanisms of particle evolution in a metastable medium. This equation is solved analytically on the basis of the saddle-point method for the Laplace integral with memory kernel. Eventually, a complete solution of the problem under consideration is written in a parametric form by means of infinite asymptotic series. The model and its analytical solutions are written in a form that is convenient for calculations of different nucleation kinetics in metastable supercooled melts and supersaturated solutions.

The zero-order terms of infinite series of our complete solution represent the asymptotic theory earlier developed by Buyevich and Mansurov [28]. Our calculations show that the zero-order solution is a very rough approximation of the complete solution in the case of phase transitions in actual metastable systems (the solution may vary several-fold). In other words, the next two summands of oscillating series make an important contribution to the complete solution for supercooling/supersaturation, density distribution function, relative number of growing particles and their mean size. As is evident from the foregoing account, the infinite series converges quickly so that its fourth term does not change the third-order solution even in the case of sufficiently small dimensionless Gibbs numbers p ∼ 0.1–1. This is because the solution quickly converges also due to small values of α_* at any fixed dimensionless coordinate and time. The main four contributions of our complete solution are found in an explicit form (table 1) for different nucleation kinetics in supercooled melts and supersaturated solutions.

It is theoretically demonstrated that the density distribution function and dimensionless supercooling/supersaturation have inflection points. These points are responsible for the different process phases at the intermediate stage of phase transitions. At the beginning, the high rate nucleation occurs. After a lapse of time this process will be accompanied by the particle growth reducing the level of metastability. The mechanism of particle coarsening becomes subsequently predominant at small supercoolings/supersaturations. Note that the theory under consideration can be used only at the intermediate stage of phase transformations. This is due to the fact that the mathematical model does not take into account of the Gibbs–Thomson effect, which plays a crucial role in determining the transition from the nucleation-growth stage to the ripening stage [1, 55]. It is important to keep in mind that the analytical solutions obtained herein for the intermediate stage can be used as initial conditions for the description of concluding (ripening) stages of the phase transition process.

As a final note, it should be mentioned that the present theory of the intermediate stage of the phase transition can be used in investigating nucleation, coarsening and agglomeration of particles in colloids, magnetic fluids and other metastable systems.

Partial support from the Ministry of Education and Science of the Russian Federation within the framework of the Federal Target Program 'Scientific and Academic—Teaching Staff of Innovative Russia' in 2009–2013 (contract No 14.A18.21.0858) and from the Russian Foundation for Basic Research (project Nos 11-01-00137, 12-01-97514, 13-01-96013-Ural and 13-08-96013-Ural) is acknowledged.

The first four coefficients in expression (30) have the form (z = z(x(t)))

$$\begin{eqnarray*} &&a_0(z)=-\Gamma \left(\frac{1}{4} \right) h(0,t)H(0),\nonumber\\ &&a_1(z)=4\Gamma \left(\frac{1}{2} \right)H(0)\left[ H(0)\left( \frac{\partial h}{\partial \nu }\right)_{\nu =0}+h(0,t)H^\prime (0)\right] ,\nonumber\\ &&a_2(z)=-8\Gamma \left(\frac{3}{4} \right) H^2(0)\left[ 2H^\prime (0)\left( \frac{\partial h}{\partial \nu }\right)_{\nu =0} + H(0)\left( \frac{\partial ^2 h}{\partial \nu ^2}\right)_{\nu =0} +h(0,t)H^{\prime \prime }(0)\right] ,\nonumber\\ &&a_3(z)=\frac{32}{3}H^3(0)\left[ 3H^{\prime \prime }(0)\left( \frac{\partial h}{\partial \nu }\right)_{\nu =0} +h(0,t)H^{\prime \prime \prime }(0)\right.\nonumber\\ &&\left. +3H^\prime (0)\left( \frac{\partial ^2 h}{\partial \nu ^2}\right)_{\nu =0} +H(0)\left( \frac{\partial ^3 h}{\partial \nu ^3}\right)_{\nu =0} \right] , \end{eqnarray*}$$

where H and its derivatives at point ν = 0 are listed in table 1 and

$$\begin{equation*} h( 0,t) =\frac{( \sqrt{z}-1)^2}{\alpha _*^2\sqrt{z}}, \end{equation*}$$

$$\begin{equation*} \left( \frac{\partial h}{\partial \nu }\right)_{\nu =0}=\frac{( \sqrt{z}-1)^2}{\alpha _*z^{3/2}}-2\frac{\sqrt{z}-1}{\alpha _* z}, \end{equation*}$$

$$\begin{equation*} \left( \frac{\partial ^2 h}{\partial \nu ^2}\right)_{\nu =0}=\frac{3}{z^{5/2}}-\frac{1}{z^{3/2}}, \end{equation*}$$

$$\begin{equation*} \left ( \frac{\partial ^3 h}{\partial \nu ^3}\right)_{\nu =0}=\frac{15\alpha _*}{z^{7/2}}-\frac{3\alpha _*}{z^{5/2}}. \end{equation*}$$

Functions f₀(z), f₁(z), f₂(z) and f₃(z) used in table 1 have the form

$$\begin{equation*} f_0=\Gamma \left(\frac{1}{4}\right)\frac{b_1\chi ( \sqrt{z}-1)^3}{12\alpha _*^3 p^{1/4}}, \end{equation*}$$

$$\begin{equation*} f_1=\frac{\sqrt{6\pi b_1}( \sqrt{z}-1)^2( \sqrt{z}-z+5)}{20\alpha _*^2\sqrt{pz}}, \end{equation*}$$

$$\begin{equation*} f_2=\Gamma \left(\frac{3}{4}\right)\frac{b_1\chi ^3}{200\alpha _*p^{3/4}z^{3/2}} ( 25+5z-66z^{3/2}+48z^2-18z^{5/2}+6z^3) , \end{equation*}$$

$$\begin{equation*} f_3=\frac{75+20z-99z^2+77z^{5/2}-111z^3+57z^{7/2}-19z^4}{100pz^{5/2}}. \end{equation*}$$

Note that

$$\begin{equation*} z=z(x(t))=1+2\alpha _*x(t),\quad x(t)=\int_{0}^{t}w(t){\,\rm d}t. \end{equation*}$$

The number of nuclei N_i and mean particle size ${\bar{L}}_i$ corresponding to the ith order solution in series of the distribution function (34) have the form

$$\begin{equation*} N_i(x(t))=\frac{1}{l_0^3}\int _0^{s_m(x(t))}F_i(x(t),s){\,\rm d}s, \end{equation*}$$

$$\begin{equation*} {\bar{L}}_i(x(t))=l_0\int _0^{s_m(x(t))}sF_i{\,\rm d}s\left( \int _0^{s_m(x(t))}F_i{\,\rm d}s\right)^{-1}, \end{equation*}$$

where i = 0, 1, 2, ... and

$$\begin{eqnarray*} &&F_i(x(t),s)=( 1+\alpha _* s) E_{1i}^{-1}(t,s)\exp \left[ -\frac{p(1-E_{1i}^2)}{E_{1i}^2}\right] ,\\ &&E_{1i}(x(t),s)=1-b_1\int _{0}^{x(t)-s-\alpha _*s^2/2}\sum _{k=0}^{i}p^{-(k+1)/4}a_k(x){\rm d}x. \end{eqnarray*}$$

The relative number of nuclei n_i and mean particle size u_i shown in figure 3 are

$$\begin{equation*} n_i(x(t))=\frac{N_i(x(t))}{N_3(x_*)},\quad u_i(x(t))=\frac{{\bar{L}}_i(x(t))}{{\bar{L}}_3(x_*)}. \end{equation*}$$

Here parameter x_* corresponds to the process time t_* when the system is completely desupercooled/desupersaturated (w = 0). This parametric dependence is determined by expression (32), where t_* = t(x_*).