Evolution of a Polydisperse Ensemble of Spherical Particles in a Metastable Medium with Allowance for Heat and Mass Exchange with the Environment

: Motivated by a wide range of applications in various ﬁelds of physics and materials science, we consider a generalized approach to the evolution of a polydisperse ensemble of spherical particles in metastable media. An integrodifferential system of governing equations, consisting of a kinetic equation for the particle-size distribution function (Fokker–Planck type equation) and a balance equation for the temperature (concentration) of a metastable medium, is formulated. The kinetic equation takes into account ﬂuctuations in the growth/reduction rates of individual particles, the velocity of particles in a spatial direction, the withdrawal of particles of a given size from the metastable medium, and their source/sink term. The heat (mass) balance equation takes into account the growth/reduction of particles in a metastable system as well as heat (mass) exchange with the environment. A generalized system of equations describes various physical and chemical processes of phase transformations, such as the growth and dissolution of crystals, the evaporation of droplets, the boiling of liquids and the combustion of a polydisperse fuel. The ways of analytical solution of the formulated integrodifferential system of equations based on the saddle-point technique and the separation of variables method are considered. The theory can be applied when describing the evolution of an ensemble of particles at the initial and intermediate stages of phase transformation when the distances between the particles are large enough, and interactions between them can be neglected.

A common important feature of all of the above processes is the highly non-linear relationship between the metastability degree of the medium (e.g., supercooling, superheating, supersaturation, undersaturation) and the evolution of crystal ensemble. For example, the degree of metastability decreases with time, which affects the dynamics of particle appearance/disappearance and the critical nucleus size of the new phase.
The evolution of particle ensemble in all the above-mentioned processes is based on similar physical laws governing the dynamics of a particular process or phenomenon in nature. Some unifying ideas about these processes can be found in Refs. [31][32][33][34][35]. In general, the evolution of a polydisperse particle system is described using a kinetic equation for the particle-size distribution function, which takes the form of the Fokker-Planck equation. In general, this equation contains sources and sinks of particles, their motion in a continuous medium as well as possible fluctuations of their growth/reduction rates. As this takes place, the evolution of metastability degree of the medium is described by a balance equation of heat or mass, which takes into account heat or mass transfer in space, release/absorption of heat or mass due to growth or degradation of particles, as well as thermal-mass exchange with the environment. Thus, the mathematical model of the process is an integrodifferential system with appropriate boundary and initial conditions. This model is complicated by the fact that particles grow/decay according to a predetermined law. Unfortunately, such a law cannot be derived analytically since the evolution of each particle is described by a moving curvilinear phase transformation boundary problem (Stefan problem) [36][37][38]. Summing up the above, we note that there are no general methods for solving the problem concerning the evolution of a system of particles (even of a simple spherical shape) in a metastable medium. Therefore the solution to all such problems represents a unique study allowing one to determine the dynamic characteristics of the system at each moment in time.
In this paper, we present the general theoretical approach based on kinetic and balance equations that describe all of the aforementioned phenomena using the same model. In addition, we discuss the main methods that can be used to solve this model containing heat and mass exchange with the environment: the saddle-point and separation of variables methods. These methods can be used to solve the problems under consideration with both the first-and second-order kinetic equations. The theoretical methods considered for solving non-linear models taking into account the heat and mass exchange with the environment can be applied to describe the dynamics of a crystal ensemble in supercooled melts and supersaturated solutions as well as the dissolution of a crystal system in undersaturated liquids.
In the case of bulk crystallization, the number of solid-phase nuclei and their dispersion determines the final structure and properties of the crystallized material. Each crystal grows from a single nucleus and therefore the number of nuclei arising in the melt/solution determines the size of the resulting crystallization grain. To achieve high mechanical and, in particular, strength properties it is desirable to obtain a fine grain structure. For this purpose, it is necessary that as many crystallization centers as possible arise in the supercooled melt (supersaturated solution) and the cooling rate must be selected in such a way that the resulting crystallization centers have the opportunity to grow. The cooling rate is directly determined by the heat exchange with the environment, the influence of which on the bulk crystallization process is investigated in this paper. The description of grain structure formation depending on the control parameters (e.g., cooling rate, supercooling, supersaturation) using modern theoretical methods is an important task that has been actively investigated in recent years to determine the formation of nano-and microstructures of alloys. One can distinguish several theoretical and practical aspects of the theory in question, which bring it to the forefront of contemporary materials science issues: (i) determination of the dispersivity of structure formation by theoretical methods, which makes cheaper and faster the procedure itself, depending on the control parameters; (ii) finding microstructure-property relations (for example, in high entropy alloys or alloys-glassformers Al50Ni50, Cu50Zr50 and Ni50Zr50).

The Generalised Model and Methods for Its Solution
Let us formulate a generalized model governing the evolution of a particulate assemblage of spherical particles in a metastable medium (e.g., liquid or gas phase). This model comprises various aforementioned phase transformation phenomena (such as bulk crystallization, dissolution, evaporation, boiling, and combustion) from the viewpoint of kinetic and balance equations for the particle-radius distribution function f and metastability degree w, which is frequently replaced by system temperature T or concentration C.

The Model
So, let us write down the kinetic equation as [39][40][41][42] ∂ f ∂t where r and t are the radial coordinate and time, x is the vector of spatial variables, f = f (r, x, t), V = dr/dt is the growth/reduction rate of particles, v is the velocity of particle motion, the function h(r, x, t) expresses the inverse of the mean residence time of particles of size r in a metastable medium, D is the coefficient of mutual Brownian diffusion of particles in space of their radii. Note that this coefficient may be defined by the "Einstein relation" in r-space or may not have an Einstein-like form (see, among others, [40][41][42][43]). Fluctuations in particle growth/reduction rates play a decisive role at the initial stage of phase transformation (when metastability degree is large enough). With allowance for such fluctuations, the nucleus radius r(t) represents a random quantity satisfying the following stochastic Equation [44] where W stands for the Wiener process. The function g(r, x, t) represents the source term of particles entering the system. So, for example, this function has a form of Dirac delta function when considering dissolution kinetics of particulate ensembles in channels [16,17]. A representative of heat/mass balance equation reads as ∂T ∂t where T(x, t) is the temperature or concentration, Q is the heat or mass sink/source term, r * and r * are the minimal and maximal sizes of particles existing in a metastable system (r * and r * are frequently replaced by 0 and ∞, respectively), and η i is constant defined by the physical meaning of the process under study. Let us especially note that Equation (3) can include additional terms when studying more complex phase transformation phenomena. For instance, dealing with combined bulk crystallization and polymerization, we include the polymerization rate on the right-hand side of Equation (3), and add an equation for polymerization degree (see, for details, [6,45]). The metastability degree w = w(x, t), which describes the ability of a system to undergo a phase transformation, is usually defined as a dimensionless positive value that varies from zero to unity. For example, considering a supercooled melt, we can define w = ∆T/∆T 0 as the ratio of the current melt supercooling ∆T to its initial supercooling ∆T 0 . By analogy, when dealing with the crystallization of a supersaturated solution, we have w = ∆C/∆C 0 , where ∆C and ∆C 0 stand for the current and initial supersaturations. Note that in the case of boiling of a superheated liquid and combustion of a polydispersed fuel, w can be defined as the dimensionless superheating of the medium [25,26,30,46,47], and in the case of dissolution of dispersed solids, w is described by the dimensionless undersaturation of the liquid [31][32][33].
The growth/reduction rate of particles V entering the governing Equations (1) and (3) substantially influences the phase transformation phenomenon. So, in the simplest case, it is proportional to the metastability degree and depends only on time t in a homogeneous system as a complex function of w, i.e., V = V(w(t)) [48][49][50]. However, it is often necessary to consider particle growth/reduction laws where V also depends on the particle radius r, i.e., V = V(r, w(t)) (see, among others, [51][52][53][54][55]). If the system is spatially inhomogeneous, V = V(r, w(x, t)).
Equations (1) and (3) should be supplemented by initial and boundary conditions that reflect the physical meaning of the process under consideration. So, for example, initial conditions can be written in the form where n 0 (x) is the initial number concentration of particles. Here we assume that the initial particle-size distribution function is f 0 (r, x), the system has an initial temperature T 0 (x), and its metastability degree is w 0 (x).
Note that the flux of particles of a given size (V f and V f − D∂ f /∂r in cases of the first-and second-order kinetic Equation (1), respectively) defines the rate I(w) of phase transformation process under consideration. This physical law determines the boundary condition of the form where r 1 is the boundary point of system transition from one phase state to another (e.g., r 1 can be equal to 0 or r * ). In the case of particle nucleation, Equation (5) determines the flux of particles crossing the nucleation barrier. So, dealing with crystal nucleation, we obtain r 1 = r * with I(w) being the rate (frequency) of nucleation [9]. Considering evaporation kinetics of a polydisperse ensemble of drops, we get r 1 = 0 with I(w) being the rate of evaporation [21]. When neglecting "diffusion term" D∂ f /∂r (neglecting fluctuations in the growth/reduction rates of particles), the kinetic Equation (1) and boundary condition (5) become simpler. This, for example, can occur when modeling bulk crystal growth in a supercooled (supersaturated) liquid [56] and intense boiling in a superheated liquid [25,26] (in the latter case, I(w) has the meaning of the rate of bubble appearance). If the kinetic Equation (1) takes "diffusion term" into account (the second-order equation with respect to the variable r), we need one more boundary condition. Such a condition can express the fact that there are no crystals of a certain size r 2 in the metastable system, i.e., So, for instance, r 2 = r p when modeling continuous operation of the crystalliser, taking into account the withdrawal process of product crystals of a given size r p [57]. When considering crystal growth without withdrawal mechanism [58] or dissolution of crystals [16,17], r 2 → ∞.
Thus, equations and boundary conditions (1)-(6) represent the closed model describing a phase transformation phenomenon accompanying the evolution of a particulate assemblage.

The Methods
Let us highlight below the main methods of theoretical modeling which allow obtaining analytical solutions to the aforementioned model equations. The first of them consists of applying the separation of variables method to the first-order kinetic Equation (1) in the absence of "diffusion term" (see, among others, [18,[31][32][33]49,50]). This method requires additional requirements to obtain the constant of variables separation (e.g., considering asymptotic behavior or comparing with experiments). Another technique in solving the first-order kinetic equation is based on the Laplace integral transform or the method of characteristics [48,59,60]. The second-order kinetic Equation (1) can be solved using the classical scheme of the separation of variables method [57] or the scheme suggested in Ref. [42] for obtaining a parametric solution. As this takes place, the distribution function found from (1) using one of these methods will be dependent on unknown metastability degree w. Therefore, we should substitute it into the balance Equation (3) and obtain a single integrodifferential equation for w. Such an equation can be solved using the saddle-point technique [33,61]. In the next section, we demonstrate the method of an analytical solution to the spatially homogeneous problem of particle nucleation and growth when the sink of product crystals plays a substantial role.

Bulk Crystal Growth
Below we consider the process of nucleation and growth of a polydisperse ensemble of spherical crystals in a metastable liquid (supercooled melt or supersaturated solution) with allowance for heat (mass) exchange with the environment and withdrawal mechanism of crystals of a given size. At first, we pay our attention to a simpler model containing the first-order kinetic equation. Then the model is extended to a more general case of the second-order kinetic equation that takes fluctuations in crystal growth rates into account.

Kinetic Equation of the First Order with a Sink of Crystals
Let us consider the process of nucleation and growth of crystals in a single-component metastable melt or solution neglecting fluctuations in the growth rates of spherical crystals (D = 0) and their external sources (g = 0). Furthermore, assuming that h = h(r), we simplify Equation (1) For simplicity, we suppose that r * = 0, r * → ∞, and rewrite Equation (3) if crystals evolve in a single-component supercooled melt as where Q T = Qρc < 0, w(t) = ∆T(t)/∆T 0 , η i = 4πL V /(ρc), L V represents the latent heat of crystallization, ρ is the melt density, and c is the specific heat. If crystals evolve in a metastable solution we should use the mass balance equation instead of heat balance, which reads as where Q = Q C > 0, w(t) = ∆C(t)/∆C 0 , η i = −4πC p , and C p stands for the saturation concentration.
The boundary condition (5) becomes where I(w) denotes the rate of nucleation and V f is the flux of crystals overcoming the nucleation barrier. The growth rate V entering this condition will be considered accordingly to the theory of steady-state crystal growth [51] and reads as for one-component supercooled melts (subscript "sm"), and for supersaturated solutions (subscript "ss"). Here k l is the coefficient of thermal conduction, and D l is the diffusion coefficient. The nucleation rate I(w) depends on the nucleation mechanism. Dealing with the Weber-Volmer-Frenkel-Zel'dovich (WVFZ) kinetics, we obtain [51] where I * and p are considered as constants for simplicity. Dealing with the Meirs nucleation kinetics, we have In this case, we also assume that I * and p are constant. Consider the simplest situation when the metastable system did not contain any crystals at the initial moment of time. In this case, the initial conditions (4) read as For the sake of convenience, we will solve the problem (7)-(15) using the dimensionless functions and variables when considering a supercooled liquid, and when considering a supersaturated liquid andĨ = I(1). Substituting (16) into (7)- (15), we rewrite the model in dimensionless form as where ψ(w) is presented in Table 1 and H(ζ) = t * h(ζr 0 ). Table 1. Analytical functions for various kinetic mechanisms, R 0 = ln −2 (1 + w −1 p ), w p = C p /∆C 0 , w 0p = 1 + w 0 /w p and C = ∆C + C p .

WVFZ Mechanism Meirs Mechanism WVFZ Mechanism Meirs Mechanism
Integrating Equation (17) and taking the initial condition (19) into account, we come to where and Heav(·) stands for the Heaviside function. Let us emphasise that Equation (20) depends on the unknown metastability degree w, which can be obtained from Equation (18).
Let us now integrate the crystal growth rate Taking the initial condition ζ = 0 at τ = µ into account, we have Here we assume that a growing crystal appears at time τ = µ. Note that the crystal of maximum radius ζ m (τ) appears at initial time µ = 0 Now we change the integration variable ζ in (18) by the new variable µ accordingly to ξ(µ) = ξ(τ) − η(ζ). Taking (20) and wdµ = −(1 + α k ζ)dζ into account, we obtain where Here the function ψ(w) is presented in Table 1. Further, for reasons of simplicity, we assume that the dimensionless crystal removal rate from the crystallizer is constant, i.e., H = const.
To find w, we evaluate the integral in (25) using the saddle-point technique [33,61] Here Γ(·) represents the gamma function and µ 0 is a high point of ψ(w(µ 0 )). Note that µ 0 = 0 if Q = 0 [51] (no heat/mass exchange with the environment). Dealing with the case Q = 0, we expand the integral term in (25) in series near the point τ = 0 and obtain Taking dψ/dµ ∼ dw/dµ into account, we conclude that the point µ 0 can be found from the equation dw/dµ = 0 at µ = µ 0 . As a result, we have from (29) This expression represents the first condition connecting w 0 and µ 0 . The second condition follows from the Cauchy problem (29) and w = 1 at µ = 0. This means that w 0 and µ 0 in expressions (27) and (28) are found. The coefficients of analytical solutions (27) and (28) are given in the Appendix A. Now keeping in mind the main contribution in the integral term (27), we arrive at where M(µ 0 ) is also written out in the appendix A. Let us especially emphasize that u in (31) depends on ξ(τ), which is a function of w(τ). Therefore, Equation (31) is integrodifferential.
The analytical solution based on expressions (20) and (32) is illustrated in Figures 1 and 2 for the WVFZ nucleation mechanism. The dynamics of desupercooling shown in Figure 1 for the WVFZ mechanism of crystal nucleation in the presence of heat sink (Q = 0.01 and Q = 0.05 in Figure 1) essentially differs from the case Q = 0 when heat/mass exchange with the environment is negligible. Such a theory with Q = 0 was developed in Ref. [51]. Indeed, the presence of a high point of w caused by Q(w) in (32) slightly increases w(τ) at the initial stage. After that w attains the maximum at a certain point τ m and then begins to decrease as a result of crystal growth in a supercooled melt (crystals release the latent heat and partially compensate w). Note that the point τ m → 0 in the limiting case Q → 0 (when the metastable liquid is not being cooled). In this case, w has no maximum point and is a monotonically decreasing function for all τ (the green dash-dotted curve shown in Figure 1) [51]. It is significant that w reduces slower as Q increases (e.g., compare the blue and red solid curves in Figure 1)). This is explained by the stronger cooling of the metastable liquid as Q increases. A similar effect is observed when product crystals are withdrawn from the crystallizer (e.g., compare the solid and dashed curves of the same color). In this case, the more H the slower w decreases.  The particle-radius distribution function Φ increases and slightly decreases up to the maximum crystal size ζ m , which is shown by the vertical line in Figure 2 at a fixed point in time τ. Namely, Φ contains a point of maximum in the presence of heat sink (when Q > 0). If Q = 0 the distribution function is a monotonically increasing function until ζ reaches a maximum value of ζ m (see also Ref. [51]). As this takes place, the distribution function is narrower and reaches large values at small times (compare two solid curves at τ = 0.2 and τ = 0.4). This is due to the fact that supercooling is greater at lower times and therefore more intense nucleation and crystal growth occur. As time passes, crystals grow in the supercooled system and the distribution function becomes wider and lower. If product crystals are withdrawn from the crystallizer (H = 3), the distribution function is smaller as compared to no withdrawal (H = 0) (compare the solid and dashed curves at the same τ). Furthermore, note that crystals of maximum radius ζ m become larger in the case of more intense cooling (as Q grows) (compare the green and red curves at the same τ).

Kinetic Equation of the Second Order with a Sink of Crystals
Let us now demonstrate the method of solution when fluctuations in crystal growth rates are taken into consideration. In this case, the kinetic equation has the form (as before, we neglect external sources of particles, i.e., g = 0) For the sake of simplicity, we consider Equation (33) with constant h = F R /V cr [62,63]. Here F R stands for the feed rate, and V cr is the volume of a metastable liquid where crystals evolve.
As before, we consider here two possible cases of a crystallizer filled with (i) a singlecomponent metastable melt and (ii) a metastable solution. In the first case, the heat balance law looks like In the second case, we have the mass balance condition Here we will show how to solve the problem dealing with the critical radius r * of new phase nuclei. Note that η i = 4πL V /(ρc) and η i = −4πC p for supercooled melts and supersaturated solutions, respectively. As before, Q T < 0 and Q C > 0. These heat/mass balances transform to corresponding balances (8) and (9) when neglecting fluctuations in particle growth rates (D = 0). To simplify the matter we consider the simplest case of the "diffusion" coefficient D = d 1 V in the space of particle radii [64][65][66], where d 1 is constant. In addition, we consider the simplest case for crystal growth rate: V = β k ∆T and V = β k ∆C when dealing with metastable melts and solutions, respectively.
The boundary conditions defining the flux of nucleating crystals and removal of product crystals with radius r p take the form Here the rate I of nucleation is defined accordingly to the WVFZ and Meirs kinetic mechanisms described in Section 3.1. The first boundary condition in (36) follows from (5) while the second condition (36) is written with allowance for total withdrawal of crystals of radius r p from the metastable liquid of the crystalliser.
Let us chose the initial conditions as Here we assume that the initial particle-size distribution function f 0 (r) is known. Thus, the governing equations, boundary, and initial conditions (33)-(37) determine the evolution of a polydisperse ensemble of spherical crystals in a crystallizer filled with a supercooled melt or supersaturated solution. This model includes heat/mass exchange with the environment through Q T /Q C and removal of product crystals of size r p from a metastable liquid through h.
Let us now use the following rescaled variables and parameters when considering the phase transition in a supercooled melt If crystals grow in a supersaturated solution we must take the substitutions ∆T → ∆C, ∆T 0 → ∆C 0 , L V /(ρc) → C p , and −Q T /(ρc) → Q C . Now rewriting the model (33)-(37) using (38), we arrive at where P(w) = w −1 exp[pψ(w)] and ψ(w) is defined in Table 1.
First we find the stationary solution when nothing depends on τ. In this case, Equation (39) leads to where subscript "s" designates the steady-state solutions. The boundary conditions (41) enable us to determine A 1 and A 2 as As this takes place, the steady-state metastability degree w s satisfies the equation following from (40) Thus, the steady-state solutions, Φ s (z) and w s , are defined by expressions (42) and (43). Equation (39) in unsteady-state case can be solved using the technique of variables separation. To make the boundary condition (41) at z = 0 homogeneous, let us make the following substitution Note that P = 1 for τ = 0 (when the metastability degree w = 1).
Substituting (44) into (39) and (41), we come to Now we separate the variables in (45) as Then the boundary conditions (46) become Z(0) − v 0 Z (0) = 0 and Z(z 0 ) = 0. Taking this into account and using the method of variables separation, we get the eigenfunctions Z j (z) and equation for eigenvalues m j as follows Next let us expand σ(z, τ) and Φ 1 (z, 0) in series of eigenfunctions as Keeping this in mind, we have from (48) Now combining (45), (50) and (51) we obtain an equation for Θ j (τ). Integration of this equation leads to Now the dimensionless distribution function Φ(z, τ) can be found from expression (44). An important point is that Φ depends on w. To find w we must use the balance Equation (40).
For this purpose, we evaluate the integral in (52) applying the saddle-point technique [33,61] for v 0 1. Introducing We evaluate the integral term in (52) keeping in mind that Λ j (τ) > 0. So, since Λ j (τ) grows and reaches a maximum at the upper limit, we obtain the main integral contribution [33,61]: Now combining (52) and (54), we obtain Substitution of (55) into (40) gives the following two-point Cauchy problem for the determination of Ω(τ): where N 0 is given in the Appendix B. Thus, we have constructed a complete analytical solution. Namely, the solution to the Cauchy problem (56) is a standard procedure programmed into various modern software packages. Therefore Ω(τ) is considered to be known. Then substituting this function into w(τ) = Ω (τ) and expressions (44) and (55), we find the metastability degree w(τ) and distribution function Φ(z, τ). An important point is that this solution is constructed taking the main contribution of the saddle-point method into account (see expansion (54)). A more precise solution to the problem, taking two terms in the asymptotic expansion into account, is given in Appendix C.
Our analytical solutions are illustrated in  In Figure 3, we compare the analytical solutions (56) found with allowance for the main contribution (54) and (C3) found with allowance for the main contribution and the first correction (C1). The main conclusion is that these solutions coincide. This indicates that the first correction does not influence the behavior of solutions essentially and we can use simpler formula (56) to plot the curves. Figure 3 also shows that the metastability degree w reduces with time and fluctuates due to the fluctuations of thermal outflux Q T (τ) = Q 0 (1 − sin( τ)/2), where = π/15. Furthermore, note that the more heat escapes the system (more Q 0 ), the more the melt cools and w increases.    (44) and (55). Here Q 0 = 0.005 and w s = 0.349. The steady-state distribution function Φ s (z) (red squares) is plotted accordingly to the analytical solution (42). Figure 4 illustrates the particle-radius distribution function versus time at fixed values of particle radii. One can easily see that the melt contains more small crystals (z = 0.01) and lesser larger particles (z = 0.05). Figure 5 shows the distribution of crystals of various sizes at fixed times τ. We can see that the particle-radius distribution function evolves with time to the steady-state distribution Φ s (z). That means, in particular, that a stationary regime establishes at large times.

Conclusions
In summary, we have formulated a generalized mathematical model of particle ensemble evolution suitable for describing various physical processes and phenomena of nature. For example, the aforementioned model can be used when describing (i) nucleation and growth of crystals in the metastable liquid or gas phase, (ii) dissolution of crystals in an under-saturated liquid, (iii) evaporation of liquid droplets in metastable gas-vapor systems, (iv) intense liquid boiling, and (v) dispersed fuel combustion.
In this study, we discussed the theoretical methods that enable constructing a complete analytical solution to the problem of nucleation and growth of spherical particles in metastable media with allowance for external heat sink (mass source) and outward crystal removal. These methods are based on the saddle-point technique to calculate the Laplace-type integral and the separation of variables method to find the eigenvalues and eigenfunctions of the corresponding boundary-value problem. As this takes place, the application of these methods depends on the order of the kinetic equation. Namely, in dealing with the first-order equation, we should use the first of these methods (Section 3.1). If fluctuations in growth rates of crystals are taken into account, we must use both of these methods (Section 3.2). The theory was developed for particular laws of crystal growth rates V(t) and nucleation kinetics I(w(t)). However, the aforementioned mathematical technique can be used when considering more complicated laws for V(t) and I(w(t)) (see, among others, [9,67]). What is more, we can restrict ourselves to the first contributions of the asymptotic solution of the saddle point technique [33,61] (as a rule, we need the main contribution and some first corrections to it).
The analytical technique described above can be used when studying similar problems about the evolution of particulate assemblages with allowance for heat/mass exchange with the environment and possible removal of particles of a given size (for example, nucleation and evolution of crystals with polymerization in a monomer metastable liquid [6,45], bulk solidification with a two-phase region [68,69], magma, lava and ice crystallization [70,71], phase transformations when producing medicines, food additives and various chemical compounds [5,49,50,72,73]). Funding: This paper comprises two parts of research studies, including (i) the model generalization and discussion of governing equations, analytical solutions, analysis of the results obtained, and (ii) checking the mathematics, numerical calculations, visualization, and discussion of the functions obtained. Two parts of this review article were supported by two financial sources. The first part was supported by the Russian Foundation for Basic Research (grant no. 20-08-00199). At the same time, the authors are grateful to the Ministry of Science and Higher Education of the Russian Federation (grant no. FEUZ-2020-0057) for the support of the second part of the research studies.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable. Data Availability Statement: All data generated or analysed during this study are included in this published article.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The metastability degree w near the maximum point µ = µ 0 is defined by The derivatives of this function near µ = µ 0 (w 0 = w(µ 0 )) take the form The coefficients of the analytical solution (28) at µ → µ 0 are given by To solve the problem more accurately, it may be necessary to use more coefficients of the asymptotic saddle point solution. In this more general case, the metastability degree w(τ) = Σ (τ) is defined by (ũ(µ 0 , µ, Σ)M(µ 0 , µ)) µ=µ 0 ,