## Abstract

The processes of particle nucleation and their evolution in a moving metastable layer of phase transition (supercooled liquid or supersaturated solution) are studied analytically. The transient integro-differential model for the density distribution function and metastability level is solved for the kinetic and diffusionally controlled regimes of crystal growth. The Weber–Volmer–Frenkel–Zel’dovich and Meirs mechanisms for nucleation kinetics are used. We demonstrate that the phase transition boundary lying between the mushy and pure liquid layers evolves with time according to the following power dynamic law: , where *Z*_{1}(*t*)=*βt*^{7/2} and *Z*_{1}(*t*)=*βt*^{2} in cases of kinetic and diffusionally controlled scenarios. The growth rate parameters *α*, *β* and *ε* are determined analytically. We show that the phase transition interface in the presence of crystal nucleation and evolution propagates slower than in the absence of their nucleation.

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

## 1. Introduction

The processes of particle nucleation, evolution, coarsening and agglomeration entirely describe the nonlinear dynamics of metastable phase transition regions in supercooled melts and supersaturated solutions. Such dynamic processes play a very important role in numerous fundamental and applied problems ranging from condensed matter physics and materials science to production of some kinds of food and pharmaceuticals in the life science, chemical and medical industries [1–8]. A transient behaviour of crystal nucleation and growth in a metastable system depends on many physical peculiarities of the phase transition process and represents a very difficult task for mathematical modelling. This explains why the nucleation process is frequently divided into several stages described by special simplified physical models. So, for instance, the initial, intermediate and concluding stages of the phase transition process may be mentioned. The initial stage usually characterizes a metastable system filled with nucleating and growing crystallites when the level of system metastability (supercooling or supersaturation) remains nearly unchanged in the beginning. The intermediate stage describes the crystal evolution and continuing nucleation of newly born particles which decrease the level of system supercooling or supersaturation. At the concluding stage when the Ostwald ripening and agglomeration processes predominate, the growing crystals are capable of coagulating and evolving at the expense of dissolution of smaller crystals. Nowadays, there are a number of theories devoted to mathematical models of crystal evolution at these stages in a stationary metastable region (see, among others, [9–18]).

However, if the bulk and directional phase transitions occur in a metastable layer synchronously in the presence of temperature and/or solute concentration gradients, the mathematical description becomes more complex. Such phase transitions take place when the metastable liquid layer filled with nucleating and evolving crystals moves as a whole due to the presence of driving force in the form of temperature or concentration gradient (see, for details, [19–21], where the steady-state solidification regime was considered). If this is really the case, the time-dependent distributions of temperature and concentration become functions of nucleation kinetics and particle growth rates. In addition, the phase transition interface (dividing regions occupied by pure liquid and metastable mushy layer) propagates in unsteady-state manner with unknown velocity and interface coordinate (see, among others, [22–26]).

The main objective of the present study is to consider a combined unsteady-state bulk and directional phase transition process when the phase interface propagates in a binary supercooled/supersaturated mixture. The theory under consideration can also be used for the theoretical description of phase transition phenomena met in such systems as heterogeneous materials [27], ferroelectrics [28], terrestrial magma oceans and lava lakes [29,30], metastable colloids and magnetic fluids [31,32] and so on.

## 2. Mathematical model

We consider a directional solidification process of a binary system, which occupies the region *ξ*>0 (figure 1). Let the initial temperature *θ*_{l} in liquid (at time *τ*=0) be greater than the phase transition temperature *θ*_{p}−*mσ*_{l} (*σ*_{l} is the solute concentration in liquid). Then the temperature at *ξ*=0 changes abruptly to the boundary temperature *θ*_{b}<*θ*_{p}−*mσ*_{l}. Furthermore, a transient layer 0<*ξ*<*Σ*(*τ*) of supercooling Δ*θ*=*θ*_{p}−*θ*−*mσ* appears and propagates with velocity d*Σ*/d*τ* into the liquid region *ξ*>*Σ*(*τ*) (*σ* is the solute concentration at 0<*ξ*<*Σ*(*τ*)). In addition, the temperature continuity at the phase interface *Σ*(*τ*) holds true, i.e. *θ*=*θ*_{l} (*θ* and *θ*_{l} represent the temperature fields in regions 0<*ξ*<*Σ*(*τ*) and *ξ*>*Σ*(*τ*), respectively). The newly born particles nucleate and evolve in the supercooled layer 0<*ξ*<*Σ*(*τ*) so that Δ*θ*=0 at *ξ*=*Σ*(*τ*). Note that the latent heat of phase transition only decreases the system supercooling Δ*θ* and does not cancel it completely in the mushy layer 0<*ξ*<*Σ*(*τ*) (figure 1).

Let us consider the growth of spherical crystallites by means of the distribution function over crystal radii *r* normed to their denumerable concentration. The kinetic, heat and mass transport equations for the distribution function *f*(*τ*,*ξ*,*r*), temperature and solute concentration fields in the mushy layer (0<*ξ*<*Σ*(*τ*)) and liquid region (*ξ*>*Σ*(*τ*)) have the form [8,15,33]
2.1
2.2
2.3where *ρ* and *ρ*_{l} are the densities of the mushy and liquid phases, *c* and *c*_{l} are their thermal capacities, λ_{l} is the thermal conductivity of the liquid phase, *κ*_{1}=4*πL*_{V}/3, *b*_{C}=4*π*(1−*k*_{0})/3, *L*_{V} is the latent heat of phase transition, *D* and *D*_{l} are the diffusion coefficients in the mushy layer and liquid, and *k*_{0} is the partition coefficient. Let us especially note that equation (2.1) does not include the ‘diffusion’ term in the space of particle radii.

Equations (2.1)–(2.3) should be supplemented by the corresponding boundary and initial conditions 2.4 2.5 2.6Note that the first condition (2.4) defines the flux of nuclei that have overcome the critical barrier.

The growth rate d*r*/d*τ* of a single spherical particle is expressed as [7,16]
2.7where *β*_{*} is the kinetic coefficient and *q*=*L*_{V}/λ_{l}. If the crystals are small enough (*r*≪(*β*_{*}*q*)^{−1}) their growth velocity is independent of their size *r* and the crystal growth scenario is termed ‘kinetic’. The opposite scenario (*r*≫(*β*_{*}*q*)^{−1}), when the growth velocity is controlled by the heat removal rate is called ‘diffusionally controlled growth’ [8].

The rate (frequency) of particle nucleation represents an exponential function of the energy barrier height. In the case of Weber–Volmer–Frenkel–Zel’dovich (WVFZ) kinetics, the nucleation rate is given by [7]
2.8where *I*_{*} is the pre-exponential factor and is the dimensionless Gibbs number (see, for details, ref. [7]). Here *γ*_{i} is the surface tension, is the initial supercooling, and *k*_{B} is the Boltzmann constant.

Also, let us write down another empirical nucleation rate (Meirs kinetics) frequently used in many industrial and engineering applications [34,35]
2.9where *I*_{*} and *p* are empirical constants. An important point is that parameters *I*_{*} and *p* in expressions (2.8) and (2.9) are different.

## 3. Analytical solutions

The nonlinear integro-differential model (2.1)–(2.9) represents a boundary-value problem with unknown mushy layer/liquid phase moving boundary *Σ*(*τ*). In order to construct its analytical solution, we use the previously developed theoretical approach [7,16] based on the well-known Laplace method. In the first instance, we evaluate the integral contributions entering in equations (2.2) using the exact analytical solution of kinetic equation (2.1) and the saddle-point method for the Laplace integral [36]. The next step is to perturb the moving boundary *Σ*(*τ*) in order to find the main contributions of the temperature, solute concentration, supercooling and *Σ*(*τ*) by means of the small parameter method.

Within this framework, we define the dimensionless variables and parameters as follows
3.1where and Δ*θ*_{l}=*θ*_{p}−*θ*_{l}−*mσ*_{l} represent the initial and liquid phase supercoolings and

The kinetic equation (2.1) and conditions (2.4), (2.7)–(2.9) in dimensionless variables (3.1) take the form 3.2and 3.3where 3.4

The solution of the problem (3.2), (3.3) can be written down as [7] 3.5

where
3.6and *η* is the Heaviside function.

Next, rewriting equations (2.2), (2.3) and boundary conditions (2.5) and (2.6) in dimensionless variables (3.1), we obtain
3.7
3.8
3.9
3.10
3.11where, for the sake of simplicity, we consider the case λ=λ_{l}, *ρ*=*ρ*_{l} and *c*=*c*_{l}. Let us especially emphasize that the integral contributions in the right-hand sides of equations (3.7) are given accordingly to the previous approach [7], where *ν* satisfies the expression *x*(*ν*,*z*)=*x*(*t*,*z*)−*y*(*s*), and

In order to evaluate the integral terms in expressions (3.7), let us use the saddle-point method for the Laplace integral [36]. For clarity’s sake, we note that ∂*g*/∂*ν*=(d*g*/d*w*)∂*w*/∂*ν*<0 for the WVFZ and Meirs kinetic mechanisms. Therefore, we conclude from expression (3.4) that *g* attains the maximum point at its boundary *ν*=0. Using equations (3.7) and to evaluate the derivatives of *w* with respect to *ν*, we conclude that the first three of them are zero at *ν*=0 and the fourth derivative at *ν*=0 is −12*χb* and −6*χb* for the WVFZ and Meirs kinetics (). Let us next retain only the main contributions of the asymptotic expansions for integrals in equations (3.7). In this case, we get [7,36]
3.12where
3.13*Γ* is the Euler gamma function, *b*_{0}=2^{7/4} and *b*_{0}=4^{3/4} for the WVFZ and Meirs kinetics, respectively.

For the sake of simplicity, we study below the kinetic (KG, *α*_{*}≪1) and diffusionally controlled (DCG, *α*_{*}≫1) growth regimes. In these cases, we obtain from (3.13)
3.14Next, expressions (3.12)–(3.14) give
3.15where
3.16

Let us seek for a solution to the problem (3.8)–(3.11), (3.15), (3.16) as a power expansion in small parameter *ε*, i.e.
3.17Substituting expansions (3.17) into equations (3.8), (3.15), taking into account that , expanding the boundary conditions at *z*=*Z*(*t*) in series and equating the terms at the same powers of small parameter *ε*, we come to the following solution in self-similar variables
3.18and
3.19where *α* and *β* are the constant coefficients describing a nonlinear behavior of the moving boundary *Z*(*t*). In addition, the temperature and concentration contributions satisfy the linear differential equations
3.20Here *Ψ*(*ζ*) is given by

The analytical solution of equations (3.20) supplemented by the boundary conditions (3.9)–(3.11) can be written down as
3.21where
Here, the growth rate constants *α* and *β* are defined by the following expressions
3.22Note that the first correction *w*_{1} to the system supercooling follows from (3.19) and (3.21) and takes the form
3.23Thus, expressions (3.18), (3.19) and (3.21)–(3.23) completely determine the analytical solution of mushy layer equations in the presence of simultaneous directional and bulk phase transitions with a moving boundary.

## 4. Numerical examples and conclusion

Figures 2 and 3 demonstrate the analytical solutions (3.17)–(3.19) and (3.21)–(3.23) for the kinetic regime of crystal growth. The dimensionless supercooling profiles plotted in figure 2 show that the processes of nucleation and evolution of crystallites in a metastable region significantly change the growth dynamics. As may be seen from the figure, the real supercooling taking into account the crystal growth process in a mushy layer (time-dependent functions *w*) is always below than the main contribution *w*_{0}, which characterizes the phase evolution in the absence of particle nucleation in a mush.

The evolutionary behaviour of the phase boundary plotted in figure 3 shows that the nucleation and growth processes in a mushy region substantially constrict its width, namely the bell-shaped dashed line taking into account particle nucleation and growth lies below the corresponding curve without nucleation (solid line which behaves as the square root of time self-similar law [38,39]). This very different transient behaviour arose from the fact that the growing nuclei release the latent heat of phase transition and compensate the mushy layer supercooling.

In summary, the theory under consideration is concerned with approximate solution of nonlinear integro-differential model that describes the transient phase transition dynamics complicated by crystal nucleation and growth processes. Taking into account the previously developed theory of particle nucleation in a motionless metastable liquid [7,16], we analytically describe the mushy layer evolution in unsteady-state manner which is caused by the crystal nucleation phenomenon. Let us especially underline here that the phase transition boundary evolves as , where *Z*_{1}(*t*)=*βt*^{7/2} and *Z*_{1}(*t*)=*βt*^{2} for the kinetic and diffusionally controlled scenarios of particle evolution. In addition, the growth rate parameters *α*>0 and *β*<0 are determined analytically too.

The present analytical approach can be easily extended to more complex phase transition regimes in the presence of moving phase transition regions (e.g. ternary and multicomponent systems [40]). For instance, it can be applied to ternary melts with allowance for the fact that the temperature relaxation time is much smaller than the corresponding relaxation times of solute impurities. The theory under consideration can also be used to describe the non-equilibrium solidification processes with a mushy layer [19,20] taking into consideration the transient behaviour of the ‘solid phase—mushy layer’ and ‘mushy layer—liquid phase’ interfaces.

## Data accessibility

This article has no additional data.

## Authors' contributions

All authors contributed equally to the present research article.

## Competing interests

The authors declare that they have no competing interests.

## Funding

This work was supported by project no. 16-08-00932 from the Russian Foundation for Basic Research.

## Footnotes

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

- Accepted August 9, 2017.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.