## Abstract

This article is concerned with a new analytical description of nucleation and growth of crystals in a metastable mushy layer (supercooled liquid or supersaturated solution) at the intermediate stage of phase transition. The model under consideration consisting of the non-stationary integro-differential system of governing equations for the distribution function and metastability level is analytically solved by means of the saddle-point technique for the Laplace-type integral in the case of arbitrary nucleation kinetics and time-dependent heat or mass sources in the balance equation. We demonstrate that the time-dependent distribution function approaches the stationary profile in course of time.

This article is part of the theme issue ‘From atomistic interfaces to dendritic patterns’.

## 1. Introduction

One of the well-known mechanisms of phase transitions in supercooled melts and supersaturated solutions is the nucleation and growth of newly born crystals [1–3]. So, the transition of a metastable phase to a thermodynamically stable phase occurs as a result of the growth of particles of a new phase in nucleating centres that appear due to some fluctuations or on heterogeneous centres of crystallization or condensation. If at the initial stage of the process individual particles can be considered as independent, then, as it develops, the physical nonlinearity caused by the influence of the collective of growing particles on the degree of metastability (i.e. on the magnitude of supercooling or supersaturation) of the surrounding matrix phase (mushy layer) becomes significant.

Another possible mechanism of the solid-phase evolution in a mushy layer occurs in directional solidification processes when the planar solid/liquid interface becomes unstable [4–7]. So, a small perturbation of the phase interface will lead to favourable growth conditions for a small bump whose tip evolves faster than neighbouring interfacial regions [8]. This eventually will lead to the formation of dendritic structures evolving in a mushy layer and compensating the thermal or constitutional supercooling [9–11]. In those parts of the mushy layer where the metastability level is small enough so that nucleation does not happen, the processes of particle coalescence (Ostwald ripening) and agglomeration may occur [12–16]. In view of the fact that it is difficult to study all possible crystal growth processes simultaneously, the mushy layer can conditionally be divided into different regions with a predominant role of one or another of the aforementioned growth mechanisms. A schematic representation of the mushy layer structure is given in figure 1. Note that throughout this paper, we theoretically study the metastable region III (the phase transition process at its intermediate stage) where nucleation and growth of crystals are of primary importance.

The available theoretical results on the evolution of the polydisperse phase in pure systems (i.e. without heterogeneous nuclei) refer either to the intermediate or to the most final stage of the phase transition [17–22]. In studies of the first stage, the assumption of the independence of individual particles is usually used. The analysis of the second stage (coalescence and coagulation) is usually based on neglecting the appearance of new nuclei [23–28]. Note that the stage in which the processes of crystal growth and the continuing nucleation play an important role at the same time have been recently studied in detail in [29–31].

The dynamics of changes in the properties of the evolving metastable system at this stage is important not only in general theoretical, but also in applied terms. This is especially true for crystallization processes from supercooled melts or supersaturated solutions, when quite often most nuclei appear precisely due to the fluctuations, and the role of heterogeneous crystallization centres is relatively small [32]. Thus, this stage determines the dispersity of products in some types of crystallizers and granulators [33].

A distinctive feature of the processes of nucleation and growth of particles in crystallizers consists in the dependence of balance equations (for heat or mass) on the intensity of external sources and the dependence of the kinetic equation for the distribution function on the removal rate of crystals [34–36]. Note that some attempts to solve the integro-differential model taking these processes into account were previously undertaken in [35–37]. However, the model equations of these papers did not take into account the ‘diffusion’ term in the Fokker–Planck kinetic equation, which plays an important role at the initial stages of particle growth, and the removal of crystals of a given size. In this paper, we eliminate these model deficiencies and construct a complete analytical solution of the integro-differential model describing the nucleation and growth of particles in a crystallizer.

This article is organized as follows. Section 2 is devoted to a statement of the problem whose complete analytical solution is presented in §3, discussions of our results as well as numerical examples are given in §4 and, finally, §5 summarizes our conclusions.

## 2. Governing equations

Let us consider the process of nucleation and growth of particles when supercooling (supersaturation) of the system is macroscopically homogeneous throughout the volume under consideration. The physical properties of solid and liquid phases will be assumed to be independent of supercooling, spatial coordinates and time *τ*. We will also neglect the processes of coagulation and breaking of growing particles.

Under these assumptions, the phase transition process is described by the following balance equations [38]:
2.1and
2.2where equation (2.1) describes the supercooled melts and equation (2.2) is valid for the supersaturated solutions. Here, *θ* and *C* represent the temperature of a supercooled melt and the solute concentration of a supersaturated solution, *ρ*_{m} and *C*_{m} are the density and specific heat of the melt, *L*_{V} and *C*_{p} designate the latent heat and concentration at saturation, *Q*_{θ}<0 and *Q*_{C}>0 represent the time-dependent external heat and mass fluxes, *f* is the distribution function, d*r*/d*τ* is the rate of particle growth, *r*_{*} is the radius of critical nuclei and *D* is a function representing the rate of particle fluctuations. Note that the exact form of coefficient *D* is a hard task of statistical physics [17,30]. For the sake of simplicity, one can assume that *D* is proportional to the growth rate [39–42]: *D*=*d*_{1} d*r*/d*τ*, where *d*_{1} is a pertinent factor.

For definiteness, we assume that the crystals evolve in a metastable system according to the kinetic growth regime. In this case, we get (see, among others, [19,22,43])
2.3where *β*_{*} stands for the kinetic parameter, Δ*θ* and Δ*C* represent the system supercooling and supersaturation, and subscripts SM and SS designate the supercooled melts and supersaturated solutions, respectively.

The distribution function *f* satisfies the following Fokker–Planck kinetic equation [17]:
2.4The dependence *h*(*r*) can be connected with the classification of crystals inside a crystallizer. For the sake of simplicity, we assume that *h*=*q*/*V* is constant (*q* and *V* are, respectively, the feed rate and total volume of a crystallizer) [44,45].

The model equations (2.1)–(2.4) should be supplemented by the corresponding initial and boundary conditions 2.5and 2.6

where Δ*θ* and Δ*θ*_{0} should be replaced by Δ*C* and Δ*C*_{0}, respectively, in the case of supersaturated solutions. Here, Δ*θ*_{0} represents the initial supercooling and the nucleation rate *I* has the form [22,29]
where the first line in curly brackets describes the Weber–Volmer–Frenkel–Zeldovich kinetics and the second line determines the Meirs kinetics. In addition, *I*_{*} and *p* are empirical constants, and *r*_{p} is the radius of withdrawn crystals where *f*=0. Note that the first boundary condition (2.6) determines the number flux of newly formed nuclei that overcame the critical barrier.

For convenience, we introduce dimensionless variables as follows:
2.7where, in the case of supersaturated solutions, Δ*θ*, Δ*θ*_{0}, *L*_{V}/(*ρ*_{m}*C*_{m}) and −*Q*_{θ}/(*ρ*_{m}*C*_{m}) should be replaced by Δ*C*, Δ*C*_{0}, *C*_{p} and *Q*_{C}, respectively.

Substituting variables (2.7) into expressions (2.1)–(2.6), we come to the dimensionless integro-differential model
2.8
2.9
2.10where and
Here, subscripts, as before, describe the supercooled melts and supersaturated solutions, and *w*_{p}=*C*_{p}/Δ*C*_{0} [30].

## 3. Analytical solution

In this section, we present the exact analytical solution of the model (2.8)–(2.10) in the steady-state solidification conditions and construct its complete solution for the transient nucleation processes.

### (a) Steady-state solutions

In the case under consideration, equation (2.9) becomes an ordinary differential equation of the second order whose solution can be written down in the form (here subscript s designates the steady-state distribution function)
3.1where and constants *C*_{1} and *C*_{2} are determined from the corresponding boundary conditions (2.10):
Here, *w*_{s} can be easily found from equation (2.8) and represents the root of equation
Below we compare the steady-state distribution function (3.1) with the unsteady-state solution constructed in the next subsection.

### (b) Unsteady-state solutions

To solve equation (2.9), we apply the method of separation of variables. For this we introduce the following auxiliary function:
3.2where *J*=1 at *t*=0 (at *w*=1).

Now combining equations (2.9), (2.10) and (3.2), we obtain 3.3and 3.4

Next substituting
3.5into (3.3) and (3.4), keeping in mind that *X*(0)−*u*_{0}*X*^{′}(0)=0 and *X*(*x*_{0})=0, we arrive at the following eigenfunctions *X*_{k}(*x*) and eigenvalues *n*_{k}:
3.6and
3.7Expansion of functions *ν*(*x*,*t*) and *F*_{1}(*x*,0) in series in *X*_{k}(*x*) leads to
3.8
3.9
3.10
3.11

Now we have from expression (3.5)
3.12Then substituting (3.5) into (3.3) and taking into account (3.8)–(3.11), we get a differential equation of the first order for *T*_{k}(*t*). Its integration after substitution into (3.12) gives
3.13where
Now expression (3.2) determines the distribution function *F* as a function of the system supercooling (supersaturation), which in turn can be found from equation (2.8). To do this, let us estimate the integral term
3.14in the case *u*_{0}≪1. Note that the derivative *S*^{′}(*t*) is positive (i.e. *S*(*t*) increases and attains the maximum value at the upper limit of integration). Therefore, the Laplace-type integral (3.14) can be estimated on the basis of the saddle-point method [46]. Taking into account only the fundamental term, we get
3.15Next substitution (3.15) into (3.14) leads to
3.16

Finally, combining (2.8) and (3.16), we come to the Cauchy problem
3.17where function *M*_{0} is determined in appendix A.

In concluding this section, we note that the Cauchy problem (3.17) determines the system supercooling (supersaturation) *w*(*t*)=*U*^{′}(*t*) and, consequently, the distribution function *F*(*x*,*t*) in accordance with expressions (3.2) and (3.16).

## 4. Discussion

Figure 2 shows the dimensionless metastability level *w*=*U*^{′}(*t*) calculated according to expressions (3.17) (when we take into account only the fundamental contribution in the saddle-point method, formula (3.15), solid curves) and expression (B3) (when we take into account the fundamental contribution and its first correction, formula (B1), circles and triangles). It is easily seen that the curves and symbols practically coincide. By this is meant that the first correction (appendix B) to the fundamental term of the saddle-point method does not influence the system dynamics substantially. In other words, one can use more simple expressions (3.15)–(3.17) to calculate the system metastability *w* and distribution function *F*. What is more, *w*(*t*) decreases with time in an oscillatory manner. This behaviour reflects the oscillatory law for thermal flux , *ω*=*π*/15, which cools the system. As one would expect, the magnitude of supercooling *w* increases with increasing the thermal flux amplitude *A*.

The distribution function dynamics for crystals of a given size is demonstrated in figure 3. As one might expect, the metastable system contains more crystals of a smaller size (smaller *x*). The distribution function profiles shown in figure 4 demonstrate that *F* approaches its steady-state solution when time increases. In other words, the idealized crystallizer under consideration tends to the stationary operating regime in the course of time.

## 5. Conclusion

In summary, a method of analytical solutions for the integro-differential Fokker–Planck and balance equations that describe the processes of particle nucleation and growth in supercooled melts or supersaturated solutions is detailed on the basis of the saddle-point method for the Laplace-type integral. The model under consideration includes a sink term in the kinetic equation, which describes the removal rate of crystals from a metastable system, and external sources of heat or mass in the balance equations. The theory is developed for arbitrary nucleation kinetics. In addition, we consider two special regimes of the Meirs and Weber–Volmer–Frenkel–Zeldovich kinetics in more detail. We show that the fundamental contribution of the saddle-point expansion gives a physically reasonable solution so that its first correction does not change the obtained solutions. Also we demonstrate that the unsteady-state distribution functions approach the steady-state solutions in course of time. The developed analytical technique can be used to describe different phase transition processes encountered in materials science and geophysics. Hence, directional solidification with a two-phase (mushy) layer [47–50], evolution of crystals in terrestrial lava lakes, magma oceans and under-ice melt ponds [51–54], coarsening of particles in colloids and magnetic fluids [55,56] may be mentioned in this regard.

## Data accessibility

This article has no additional data.

## Authors' contributions

All authors contributed equally to this article.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by the Russian Foundation for Basic Research (grant no. 16-08-00932).

## Appendix A

The right-hand side of equation (3.17) takes the form where the following designations are introduced:

## Appendix B

Let us now determine the next contribution entering in the asymptotic expansion of the Laplace-type integral (3.15). Taking into account the general expressions found in [46], one can get B 1where Note that the first term in parenthesis in the right-hand side of equation (B.1) coincides with the fundamental contribution used in expansion (3.15).

Taking this into account, we have instead of expression (3.16)
B 2where
Now combining (2.8) and (B.2), we arrive at
B 3Here, *U*_{sd} can be easily found from the balance condition (2.8) after substitution *F*(*x*,0)=*F*_{0}(*x*) and *w*=*U*^{′}, and
where

## Footnotes

One contribution of 16 to a theme issue ‘From atomistic interfaces to dendritic patterns’.

- Accepted September 21, 2017.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.