Morphological/Dynamic Instability of Directional Crystallization in a Finite Domain with Intense Convection

: This study is devoted to the morphological/dynamic instability analysis of directional crystallization processes in ﬁnite domains with allowance for melt convection. At ﬁrst, a linear instability theory for steady-state crystallization with a planar solid/liquid interface in the presence of convection was developed. We derived and analyzed a dispersion relation showing the existence of morphological instability over a wide range of wavenumbers. This instability results from perturbations arriving at the solid/liquid interface from the cooled wall through the solid phase. Also, we showed that a planar solid/liquid interface can be unstable when it comes to dynamic perturbations with a zero wavenumber (perturbations in its steady-state velocity). A branch of stable solutions for dynamic perturbations is available too. The crystallizing system can choose one of these branches (unstable or stable) depending of the action of convection. The result of morphological and dynamic instabilities is the appearance of a two-phase (mushy) layer ahead of the planar solid/liquid interface. Therefore, our next step was to analyze the dynamic instability of steady-state crystallization with a mushy layer, which was replaced by a discontinuity interface between the purely solid and liquid phases. This analysis showed the existence of dynamic instability over a wide range of crystallization velocities. This instability appears in the solid material at the cooled wall and propagates to the discontinuity interface, mimicking the properties of a mushy layer. As this takes place, at a certain crystallization velocity, a bifurcation of solutions occurs, leading to the existence of unstable and stable crystallization branches simultaneously. In this case, the system chooses one of them depending of the effect of the convection as before. In general, the crystallizing system may be morphologically/dynamically unstable when it comes to small perturbations arriving at the phase interface due to ﬂuctuations in the heat and mass exchange equipment (e.g., ﬂuctuations in the freezer temperature).


Introduction
The dynamics of an interfacial boundary (crystallization front) separating the solid and liquid phases is responsible for the formation and growth of solid phase protrusions, dendritic crystals, absorption and redistribution of impurities and finally for the development of different microstructures (for example, cellular, banded, mixed-type and irregular structures) in the crystallized material [1][2][3][4][5][6][7]. In order to control the process of microstructure development by means of the physical and operating parameters of the crystallization process, the laws controlling the formation of one or another type of structure as well as the transitions between various types of structural formations must be fully established. One of the processes influencing the structuring of materials is the fluctuations in various physical quantities (e.g., temperature fluctuations always existing in nature; hydrodynamic fluctuations in the melt flow rate caused by convection, natural causes or the type of laboratory setup; fluctuations in the impurity concentration in liquid caused by hydrodynamic fluctuations or impurity supply from the outside; mechanical fluctuations in the entire crystallizing system caused by seismic processes in nature or the natural background of human activity in a laboratory setup; etc.). These fluctuations result in perturbations in the temperature and concentration fields as well as hydrodynamic perturbations in the fluid velocity, leading to perturbations in the solid/liquid interphase boundary (the crystallization front). The development of such morphological or dynamic perturbations is facilitated by the inhomogeneities in the temperature and concentration fields, convection and thermal or constitutional supercooling. As a result, such perturbations completely modify the crystallization process and lead to the formation of different microstructures (for example, the morphological instability of the solid/liquid interface is responsible for the development of cellular structures; its dynamic instability leads to the formation of banded structures, and the evolution of the dendritic forest in general generates a twophase (mushy) layer) [8][9][10][11][12][13]. Note that the linear analysis of morphological stability with applications for crystallization problems was first used by Mullins and Sekerka [14][15][16]. Then, their technique was extended to describe the different features of crystallization phenomena in refs. [17][18][19][20][21][22][23][24]. Namely, the effects of (1) anisotropic liquid-solid interfaces [17], (2) fluctuating and periodic growth rates [18,19], (3) transient solidification [20], (4) selfsimilar growth [21], (5) rapid crystallization [22], (6) a mushy layer with a changeover from instability [23] and (7) convective instability with a forced flow [24] were investigated.
In addition, on the one hand, convection can equalize temperature and concentration distributions in a liquid, and on the other hand, convective flows or cells can create cold or hot spots in certain areas of the solid/liquid interface [25,26]. As this takes place, in colder spots, the conditions of the preferential growth of solid-phase protrusions are created; i.e., a morphological instability can be developed. At present, the effect of convection on morphological/dynamic instability has been weakly investigated due to the strong nonlinearity of mathematical models, the variety of convective flows and the different boundary conditions at the phase interface in the presence of convection. For instance, a simplified mode of the heat and mass transfer laws was applied to investigate the effect of a plane-parallel fluid current on the interfacial boundary stability in ref. [27]. The authors of refs. [28][29][30] applied the boundary layer techniques to omit small nonlinear terms and simplify the model equations. The spatially periodic fluid currents were considered in ref. [31] to model the localized morphological patterns. A mathematical model in an unbounded area of space relying on conductive heat and mass transfer laws was applied to study the effects of fluid currents in refs. [32,33]. This theory showed that the dispersion relation and the neutral stability curve in the presence of convection depend significantly on the extension rate at the solid/liquid interface.
In this study, we substantially develop the morphological/dynamic stability theory in the presence of convective flows. First of all, our theory is built in a limited region of space, which allows us to model perturbations coming from the solid and liquid phases (e.g., from the solid boundaries of laboratory facilities or natural processes). Secondly, two types of instabilities are investigated below: morphological and dynamic. Initially, a quasistationary crystallization process with a planar solid/liquid interface separating the purely solid and liquid phases is considered. As a result of the development of the morphological instability, the planar crystallization front breaks up and a two-phase (mushy) layer between the purely solid and liquid phases is formed. This conclusion is confirmed by experiments in which aqueous solutions of isopropanol were cooled and crystallized from above in the presence of convection [34]. The mushy layer may be unstable when it comes to dynamic perturbations in the constant crystallization rate. Therefore, we further analyzed the dynamic instability of the quasistationary crystallization process with a two-phase layer. This analysis showed that dynamic perturbations are unstable and lead to fluctuations in the rate of directional crystallization.

Directional Crystallization with a Planar Solid/Liquid Interface
Let us first consider the steady-state directional crystallization of a binary liquid (melt or solution) with a planar solid/liquid interface along the spatial coordinate ζ caused by temperature gradients in the solid and liquid phases. The process under consideration is triggered by the cooling of a liquid from above and controlled by an intense convection in the liquid (Figure 1a). Using the coordinate system moving with a constant velocity u s together with the solid/liquid interface, we have the following heat transfer equations in both the phases: Here, θ l and θ s are the temperatures in the liquid and solid phases, respectively; τ is the time; v is the vector of the fluid velocity; D l and D s are the thermal diffusivities in the liquid and solid phases, respectively; u s is the crystallization velocity; and Σ (τ) is the solid/liquid interface coordinate (Σ (τ) = 0 when considering the unperturbed steady-state crystallization scenario). Given sufficiently intense convection in the liquid to equalize the concentration distribution therein, we do not consider the process of impurity diffusion [34]. At the phase interface, the temperatures from both sides of the solid/liquid interface are equal to the sum of the phase transition temperature θ * of the pure substance and the interface curvature term ΓH. In addition, the difference in thermal fluxes defines the crystallization heat L V released at the interface [34,35], i.e., where Γ = θ * γ/L V is the Gibbs coefficient, H is the interface curvature, γ is the solid/liquid interfacial energy, L V is the latent heat parameter, n is the normal vector of the interface, c l is the specific heat, θ ∞ is the temperature of the liquid far from the crystallization front, k s is the thermal conductivity of the solid and j(θ l ) is the convective heat flux transferred from the liquid to the solid phases. Note that H = 0 in the case of a planar solid/liquid interface and H ≈ ∇ 2 Σ in the case of small morphological perturbations (linear theory). The convective heat flux j(θ l ) (J m −2 s −1 ) is given by [34,35] where λ is a dimensionless constant, is the coefficient of thermal expansion and g (m s −2 ) is the acceleration due to gravity. Let us especially note that α can be a function of θ l , i.e., α(θ l ) = 10 −4 (2.25 + 0.15θ l ) for the isopropanol solution [34]. Note that Equation (5) was justified in ref. [35] based on the general principles of convective heat transfer, and the parameter λ was found for the isopropanol solution by a particular experiment in ref. [34].

Morphological Instability of a Planar Solid/Liquid Interface with Intense Convection in Liquid
Let us assume that heat transfer predominates in the direction of the spatial axis ζ. In this case, the steady-state temperature profiles in the liquid and solid phases are given by θ lo = θ lo (ζ) and θ so = θ so (ζ) (the subscript "o" denotes the steady-state solutions).
The boundary condition (4) enables us to find the steady-state temperature gradient in a solid at ζ = 0: where θ int = θ * for the planar solid/liquid phase interface. The heat transfer Equations (1) and (2) allow us to obtain the second derivatives at ζ = 0 (here we use the no-slip condition v = 0 at ζ = 0) in the form of where G l = dθ lo /dζ at ζ = 0. Now we morphologically perturbed the planar solid/liquid interface as Here, ξ and η are the Cartesian coordinates directed perpendicular to the solidification axis ζ; k ξ and k η are the perturbation wavenumbers along these directions, respectively; i is the imaginary unit; and ω is the amplification rate (frequency) of the perturbations. The phase interface perturbation Σ is in accordance with the temperature perturbations in the solid and liquid, i.e., θ s = θ sA E exp(β s ζ) and θ l = θ l A E exp(β l ζ), where Σ A , θ sA and θ l A designate the perturbation amplitudes and β s and β l are the perturbation amplification/damping coefficients found below. The linear theory under consideration implies that |θ s | θ so and |θ l | θ lo . The next step is to expand the boundary conditions (3) and (4) in a Taylor series in the vicinity of the unperturbed solid/liquid interface ζ = 0. Keeping in mind only the linear terms of the perturbations, we arrive at the following expressions at ζ = 0 (for details, see Appendix A): where k 2 h = k 2 ξ + k 2 η and j = dj/dθ l at θ l = θ int , i.e., The heat and mass transfer Equations (1) and (2) lead to the following equations at ζ = 0: where V ζ = −∂v ζo /∂ζ at ζ = 0 represents the extension rate (v ζo denotes the steady-state ζ component of the fluid velocity). Combining expressions (9) and (12), we come to β l in the form of An important point is that the signs ± in the expressions for β s and β l define the direction of the perturbation amplification/damping. So, for example, if a perturbation appears in the solid phase (ζ < 0) at some distance from the phase interface ζ = 0, it grows/decays in cases of positive/negative β s . When dealing with the liquid phase (ζ > 0), we have a similar behavior. Namely, a perturbation originating at some distance ζ > 0 in the liquid propagates to the phase interface in the −ζ direction and decays/grows if there is a positive/negative β l . In real laboratory setups or natural processes, this corresponds to a limited region of the crystallization process: the solid walls are located at certain distances ζ < 0 in the solid material and ζ > 0 in the liquid phase. Now combining Equations (8) and (9), we find the dispersion relation (ω(k h )) as follows: where Figure 2a shows the dispersion curves calculated according to the dispersion relation (15) for the isopropanol solution experimentally studied in ref. [34]. Analyzing Equation (15), we found only real solutions. As is easily seen, the directional crystallization process in the presence of intense convection is morphologically unstable (ω increases with an increasing wavenumber k h and decreasing steady-state velocity u s ). Such a behavior follows from ω > 0 calculated in a wide range of wavenumbers k h . As this takes place, all points shown in Figure 2a are described by β s < 0. Physically, it means that a perturbation in the solid phase temperature appearing at some distance from the crystallization front at ζ < 0 (e.g., a temperature fluctuation arising on the cooled wall) decreases with an increasing ζ and perturbs the solid/liquid interface at ζ = 0. Note that β s < 0 and β l < 0 describe the damped perturbations propagating from the cooled boundary at ζ < 0 into the liquid phase at ζ > 0. It is important to emphasize that the revealed morphological instability is a consequence of the perturbations appearing on the cooled wall at ζ < 0. In the case of an infinite crystallization domain, such a solution does not exist due to the requirement of damping perturbations at ζ → −∞. As a result of this mechanism, morphological perturbations develop in the liquid at the phase interface, and a mushy layer filled with dendrite-like structures appear, as was demonstrated in experiments [34]. In addition to the morphological perturbations where k h = 0, dynamic perturbations where k h = 0 may exist in the crystallizing system, which represent the perturbations in the steady-state crystallization velocity u s . Figure 2b shows the amplification rate versus u s at different temperatures θ ∞ . We found two branches of the solution, one branch of dynamic instability (ω > 0, the solid line in Figure 2b) and the other one of stability (ω = 0, the dashed line in Figure 2b). Note that the root ω = 0 of Equation (15) can be easily found analytically for dynamic perturbations where k h = 0. In addition, Σ = Σ A and dΣ /dτ = 0, i.e., the system has zero velocity perturbation and crystallizes with an unperturbed steady-state rate u s . Which of these two branches of solution will be realized depends on whether convection amplifies (by creating irregularities in the temperature and concentration fields) or attenuates (by equalizing the temperature and concentration) the dynamic perturbations near the solid/liquid interface. The answer to this question depends on the nature of the convective currents and the arrangement of the crystallizing facility. Note that the mushy layer often occurs in various geophysical phenomena of ice freezing and magma solidification, as well as in metallurgical and chemical processes of the equilibrium and nonequilibrium crystallization of melts and solutions [36][37][38][39][40][41]. This explains the need to develop a theory of dynamic stability of such processes under the influence of convection. When a mushy layer has appeared between purely solid and liquid phases, small temperature perturbations can produce a new oscillatory crystallization scenario when the mushy layer is dynamically unstable and oscillates near its steady-state crystallization velocity. Such a process substantially changes the impurity distribution in solid material and leads to the phenomenon of layered impurity liquation appearing as a result of the dynamic oscillations of a mushy layer. To describe this effect within the framework of intense convection in liquid, we develop the aforementioned model and analyze it against small dynamic perturbations. To simplify the matter, we replace a real mushy layer with the discontinuity interface between purely solid and liquid phases that reflect the properties of a real mushy layer by means of a new boundary condition. This condition represents the equality of the phase-transition temperature gradient and the liquid-phase temperature gradient at the discontinuity interface. In addition, this condition defines the fact that no supercooling exists ahead of the discontinuity surface (the mushy layer) in the liquid. Such an analysis of dynamic stability is carried out below in the spirit of previously developed theories without convection [42,43].

Dynamic Instability of a Discontinuity Solid/Liquid Interface Mimicking the Properties of a Mushy Layer
First of all, we consider a narrow quasiequilibrium mushy layer appearing ahead of the planar solid/liquid interface. The latent heat of solidification completely compensates for the supercooling in a mush, which is replaced by a discontinuity interface between the purely solid and liquid phases [44,45] (Figure 1b). Since impurities accumulate in the interdendritic spacing of a mushy layer, where convection is retarded and cannot equalize the impurities, it is necessary to include the impurity concentration C l in the mathematical model of directional crystallization. Thus, the process model consists of the heat-conduction Equations (1) and (2) and an impurity-diffusion equation in the liquid phase (diffusion in the solid material is traditionally neglected): where D C is the diffusion coefficient and Σ (τ) represents the dynamic perturbations in the discontinuity interface that replaces a mushy layer (Σ (τ) = 0 when dealing with the steady-state crystallization scenario). Since the phase transition temperature is dependent on the impurity concentration, the boundary condition (3) should be rewritten as follows: where the function f (C l ) is determined from the phase diagram ( f (C l ) = mC l for the linear phase diagram, with m being the equilibrium liquidus slope). The heat balance condition (4) is satisfied at the discontinuity interface Σ (τ). To close the mathematical model, we have the following quasiequilibrium condition of the mushy layer, which is replaced by a discontinuity interface [44,45]: where f (C l ) = d f /dC l and f = m in the case of the linear liquidus equation.
Thus, the mathematical model of directional crystallization with a mushy layer that is replaced by a discontinuity interface consists of the heat and mass transfer Equations (1), (2) and (18) as well as the boundary conditions (4), (19) and (20).
In the case of steady-state solidification, Σ (τ) = 0 and the expressions (6) and (7) take place. What is more, the impurity diffusion Equation (18) gives the following at ζ = 0: The subscript "o", as before, designates the steady-state solutions.
Since the mushy layer can be dynamically unstable as a whole object when its steadystate velocity oscillates around u s , let us analyze the dynamic stability of crystallization with a discontinuity surface below (Figure 2b). In this case, the perturbations have the same form with E = exp(ωτ) and C l = C l A E exp(βζ), where β is the perturbation amplification/damping coefficient for the impurity concentration and C l A stands for the amplitude of its perturbations. Now substituting perturbations into the boundary conditions (4), (19) and (20), we obtain four equations at ζ = 0 for the perturbation amplitudes. As before, Equations (8) and (10) take place. The other two equations are as follows: where f = 0 in the case of the linear liquidus equation.
Equations (1), (2) and (18) enable us to express β l (ω), β s (ω) and β(ω) as follows: As before, the signs of β l , β s and β determine the direction of the perturbations propagating to the discontinuity interface from the solid or liquid phases. Now eliminating the perturbation amplitudes from Equations (8), (10), (22) and (23), we arrive at the following expression for ω: The sign of ω (or the real part of ω) in Equation (27) defines the process stability/instability when it comes to small dynamic perturbations (stability occurs at ω < 0 and instability occurs at ω > 0).
Analyzing expression (27) for the isopropanol solution, we found only the real roots of this equation. They are shown in Figure 3 in the plane ω(u s ) for different temperatures θ ∞ . As is easily seen, the amplification rate can be either positive or negative for any given value of θ ∞ . Moving from left to right along the horizontal axis u s , we see that at first, ω is positive (solid line) up to the bifurcation point of the solution, which is marked by the symbol "black square, " (Figure 3a). Then, at the bifurcation point, the solution splits into two branches marked with dashed lines (the unstable branch continues in the region ω > 0 and the stable branch jumps into the region ω < 0). Then, both of these branches remain in their own semiplanes. As this takes place, the stable branch (ω < 0) asymptotically approaches zero from below (Figure 3b). An important point is that β s < 0, β l < 0 and β < 0 for all points in Figure 3. As before, it means that a perturbation appearing at the cooled wall at ζ < 0 propagates across the discontinuity interface to the liquid phase and leads to dynamic perturbations for a single solution to the problem (the solid line in Figure 3a). When the crystallization velocity u s is to the right of the bifurcation point, there are two possibilities: (i) the process jumps into the stable region ω < 0 or (ii) the process stays in the unstable region ω > 0. Which path the system chooses depends on convection, which can either (i) equalize the temperature and concentration fields in the liquid and lead to stability or (ii) create local irregularities in these fields near the discontinuity surface and lead to instability. In other words, the perturbation coming from the cooled wall can either attenuate (ω < 0, stability) or enhance (ω > 0, instability) due to convection. In the latter case, the discontinuity interface becomes dynamically unstable, and the quasistationary crystallization process with constant velocity u s is destroyed.

Conclusions
In summary, we study how small morphological/dynamic perturbations influence the directional crystallization process in a finite domain with allowance for vigorous convection in liquid. First of all, we investigate whether the planar solid/liquid interface moving in a steady-state manner is stable when it comes to small morphological perturbations in the temperature and concentration fields, as well as to perturbations in the fluid flow velocity. To perform this, we derive the dispersion relation defining the amplification rate (the frequency of perturbations) as a function of the perturbation wavenumber and other system parameters. Analyzing this relation for the isopropanol solution, we conclude that ω is positive in a broad range of wavenumbers. This solution corresponds to the negative values of the perturbation amplification/damping coefficients β s and β l . The negative signs of these coefficients show that the morphological perturbations propagate along the crystallization direction ζ. Namely, when the temperature fluctuates on a cooled wall (at ζ < 0) at some distance from the phase interface, this perturbation propagates to the liquid phase and makes the solid/liquid interface morphologically unstable. In addition, we show that the planar solid/liquid interface can be unstable when it comes to dynamic perturbations with a zero wavenumber (perturbations in the steady-state crystallization velocity u s ) in a broad range of u s . However, the unstable branch of solutions simultaneously coexists with the stable one. As this takes place, the crystallizing melt/solution chooses one of them depending on the influence of convection, which can either strengthen or weaken a dynamic perturbation coming from the solid phase (cooled wall). This morphological/dynamic instability evolves with time and leads to the formation of a two-phase (mushy) layer between the purely solid material and liquid (see the experiments in ref. [34]). The crystallization process then takes place in the presence of this mushy layer.
To study the stability of such a process when it comes to dynamic perturbations (perturbations in the steady-state crystallization velocity with a mush), we carry out a linear dynamic stability analysis with convection. The result of this theory, where a mushy layer is replaced by a discontinuity interface reflecting its properties, is the equation for the amplification rate as a function of the process and physical parameters. Analyzing this equation, we see that the dynamic perturbations can evolve from solid to liquid material. Namely, if a perturbation appears on the cooled wall at ζ < 0, it decreases and propagates to the liquid phase at ζ > 0 with β s < 0, β l < 0 and β < 0. This perturbation can evolve with time and lead to dynamic instability in a wide range of crystallization velocities u s . In addition, a bifurcation of solutions occurs at a certain velocity u s , and the system has two branches of solutions (unstable and stable) that coexist simultaneously. The system chooses one of them depending on convection, which strengthens or weakens the perturbations in cases of instability and stability, respectively. As a result, the discontinuity interface mimicking a mushy layer can be unstable when it comes to small dynamic perturbations in the crystallization velocity; i.e., the steady-state process with mush becomes broken and a more complex crystallization scenario with an unsteady velocity and variable mushy layer thickness can occur.
Let us summarize the main assumptions used in the theory under development below. (1) The steady-state distributions of the temperature and impurity concentration depend on only one spatial coordinate ζ directed along the crystallization process. (2) The convective heat flux is defined by using the experimentally valid expression (5). (3) Mass transfer in the liquid phase is described by a convective-type equation of impurity diffusion (18), which is valid for slow local-equilibrium crystallization processes (at u s V D , where V D is the solute diffusion velocity in bulk liquid [47,48]). (4) The temperature of the phase transformation at the solid/liquid interface does not depend on atomic kinetics, which is the case for nonhigh crystallization velocities (u s V D ). (5) All the thermophysical parameters of the solid and liquid phases are assumed to be constant near the phase transformation temperature. (6) The morphological/dynamic perturbations are sufficiently small compared to the characteristic steady-state solutions, which allow us to expand the required functions in the Taylor series near the solid/liquid interface.
As a result of this work, the main conclusion can be summarized as follows. A planar solid/liquid interface (mushy layer) can be unstable when it comes to small morphological/dynamic perturbations in the presence of vigorous convection in liquids when considering a finite crystallization domain. The evolution of instability is schematically illustrated in Figure 4. As this takes place, the key factor here is the finite process domain with solid walls, which is characteristic of any real crystallization phenomenon occurring in nature, laboratory or industrial facilities. Namely, the origin of perturbations occurs on solid walls due to fluctuations in the heat and mass transfer facility, which lead to morphological/dynamic instability. If we were to carry out this study in an infinite domain, we would be forced to discard the solutions found from the conditions of the bounded temperature and concentration perturbations at an infinite distance from the solid/liquid interface in both the phases. This allows us to conclude that the finiteness of the crystallization domain is one of the main factors influencing the morphological/dynamic instability.  This means that a number of theories of stable/unstable crystallization with convection (and perhaps even in its absence) need to be revised to take a limited crystallization domain into account, such as the theory of the stable growth mode of a dendritic crystal tip, which allows for the selection of the tip velocity depending on the tip curvature and supercooling [49][50][51][52]. Another example is the morphological stability of the ice/ocean interface when it comes to morphological perturbations considering a finite depth of liquid [53][54][55]. v ζ = v ζo + ∂v ζo ∂ζ Σ + v ζ = 0, at ζ = 0, or v ζ = V ζ Σ , at ζ = 0.
By substituting the perturbations here, we come to expression (14) at ζ = 0.