1 Introduction
The solidification of multicomponent mixtures plays an important role in many natural and industrial systems. These solidification processes are often accompanied by the formation of so-called mushy layers, which are regions of an evolving dendritic solid matrix surrounded by the melt phase. Mushy layers influence solidification growth rates and affect the transport of heat, mass and solute. Further, processes that occur within the interior of these mushy layers can have a profound impact on features of interest in the adjacent liquid and solid phases. For example, the mushy layer that forms during the growth of sea ice plays an important role in the prediction of the salt flux into the ocean, an important component of accurate ocean and climate models (e.g. Notz & Worster Reference Notz and Worster2009; Rees Jones & Worster Reference Rees Jones and Worster2014). In the context of industrial materials processing, it is well known that the appearance of channel or ‘freckle’ defects in cast solids have origins tied to convective instabilities that occur within the mushy layer (e.g. McDonald & Hunt Reference McDonald and Hunt1969, Reference McDonald and Hunt1970; Copley et al. Reference Copley, Giamei, Johnson and Hornbecker1970; Giamei & Kear Reference Giamei and Kear1970; Sample & Hellawell Reference Sample and Hellawell1982, Reference Sample and Hellawell1984; Sarazin & Hellawell Reference Sarazin and Hellawell1988). Thus, the examination of mushy layer dynamics, even in isolation from the parent liquid phase and resultant solid phase, plays a useful role. From a fluid dynamics perspective, mushy layers in multicomponent systems constitute a rich area in which novel interactions between convection, heat and solute diffusion and phase change may occur. Recent work (Anderson et al. Reference Anderson, McFadden, Coriell and Murray2010; Guba & Anderson Reference Guba and Anderson2014) on linear stability analyses of an ideal, isolated mushy layer that forms during the solidification of ternary solutions has identified particularly novel convective instabilities that occur despite the lack of statically unstable density gradients with respect to temperature and composition. The goal in our present work is to further explore, quantify and understand these novel convective instabilities within the context of a weakly nonlinear stability analysis.
Convective instabilities that occur within mushy layers formed during the solidification of binary alloys have, and continue to receive, a significant amount of attention. Several excellent reviews have been written (e.g. Worster Reference Worster, Davis, Huppert, Muller and Worster1992, Reference Worster1997, Reference Worster, Batchelor, Moffatt and Worster2000; Davis Reference Davis2001). Convection motion that occurs in a binary alloy mushy layer, a reactive porous medium, can be compared to classical models of convection within inert porous layers with differentially heated upper and lower boundaries (e.g. Horton & Rogers Reference Horton and Rogers1945; Lapwood Reference Lapwood1948; Palm, Weber & Kvernvold Reference Palm, Weber and Kvernvold1972; Nield & Bejan Reference Nield and Bejan1998). One of the key differences between the reactive and inert cases is that when convection occurs in a reactive mushy layer, the local solid fraction evolves in response to heat and solute transport – a well-known consequence of this interaction is chimney convection (e.g. Worster Reference Worster1991; Schulze & Worster Reference Schulze and Worster1998; Katz & Worster Reference Katz and Worster2008; Wells, Wettlaufer & Orszag Reference Wells, Wettlaufer and Orszag2010, Reference Wells, Wettlaufer and Orszag2013). Another aspect of this coupling between phase transformation and the chemical thermodynamics within the mushy layer is that the liquid and solid are in local thermodynamic equilibrium and so the temperature and liquid composition are directly linked through the liquidus relationship of the binary phase diagram. This direct coupling of temperature and composition, which is not present during convection of a binary mixture in an inert porous medium, eliminates the possibility of double-diffusive convection (e.g. Turner Reference Turner1973) occurring within the binary alloy mushy layer.
Mushy layers that arise in multicomponent solidification have also been of interest, particularly due to their importance in metallurgy (e.g. Fujii, Poirier & Flemings Reference Fujii, Poirier and Flemings1979; Felicelli, Poirier & Heinrich Reference Felicelli, Poirier and Heinrich1997, Reference Felicelli, Poirier and Heinrich1998; Krane & Incropera Reference Krane and Incropera1997; Krane, Incropera & Gaskell Reference Krane, Incropera and Gaskell1997, Reference Krane, Incropera and Gaskell1998; Schneider et al.
Reference Schneider, Gu, Beckermann, Boettinger and Kattner1997). The solidification of aqueous ternary solutions, like their binary solution counterparts, have been the focus of numerous laboratory investigations (e.g. Aitta, Huppert & Worster Reference Aitta, Huppert and Worster2001a
,Reference Aitta, Huppert, Worster, Ehrhard, Riley and Steen
b
; Bloomfield & Huppert Reference Bloomfield and Huppert2003; Thompson et al.
Reference Thompson, Huppert, Worster and Aitta2003b
). One notable feature in the solidification of the ternary solution,
$\text{H}_{2}\text{O}$
–
$\text{KNO}_{3}$
–
$\text{NaNO}_{3}$
, considered by Aitta et al. (Reference Aitta, Huppert and Worster2001a
,Reference Aitta, Huppert, Worster, Ehrhard, Riley and Steen
b
) and Thompson et al. (Reference Thompson, Huppert, Worster and Aitta2003b
), is the appearance of two distinct mushy layers. The primary mushy layer is distinguished by the solid phase being composed of primarily a single component (e.g. ice) and the rejection of the other two components into the liquid. The secondary mushy layer is distinguished by the solid phase being composed of two components of the solution (e.g. ice and solid
$\text{KNO}_{3}$
) with the third component being rejected back into the liquid. In such a case, the temperature and two compositions within the primary mushy layer are coupled directly and lie along a liquidus surface (surface of onefold saturation) in the ternary phase diagram (e.g. Cocks & Brower Reference Cocks and Brower1974; Lupis Reference Lupis1993; Smallman & Bishop Reference Smallman and Bishop1999). In the secondary mushy layer the temperature and compositions lie along the intersection of two liquidus surfaces (cotectic boundary, or line of twofold saturation). Theoretical models for the growth of such mushy layers have been developed for diffusion-controlled growth (e.g. Anderson Reference Anderson2003; Thompson, Huppert & Worster Reference Thompson, Huppert and Worster2003a
; Alexandrov & Ivanov Reference Alexandrov and Ivanov2009). Scenarios for convective instabilities that may occur in the absence of solute diffusion in the primary and secondary mushy layers have also been explored (e.g. Anderson & Schulze Reference Anderson and Schulze2005; Riahi Reference Riahi2014).
While the ternary alloy double-mushy layer configuration gives rise to the possibility of convective interactions between two mushy layers as well as the adjacent liquid region, it has been shown recently that novel and unexpected instabilities originating within the primary mushy layer alone may occur (Anderson et al. Reference Anderson, McFadden, Coriell and Murray2010; Guba & Anderson Reference Guba and Anderson2014). Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010) explored the possibility of double-diffusive-type convection occurring within the primary mushy layer using linear stability theory. The temperature and compositions in the primary mushy layer are constrained to the liquidus surface (the analogue of the liquidus constraint present in binary mushy layer models). This extra degree of freedom was shown to lead to double-diffusive-type instabilities in the primary mushy layer in the sense that a statically stable configuration (motionless system with density increasing with depth in the ternary mushy layer) could be dynamically unstable if the slower diffusing agent (typically solute rather than heat) was unstably stratified. However, perhaps more surprisingly, Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010) also showed that convective instabilities could occur even in the event that all agencies (heat and solutes) were statically stably stratified so that no positive buoyancy was present in the system. Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010) gave physical explanations of the mechanisms behind these instabilities in the form of parcel arguments that required disparate solute diffusion coefficients associated with the two solutes that were both much smaller than the thermal diffusion coefficient. Their examples focused on the cases in which the buoyancy was associated with the faster diffusing solute, but in general the buoyancy is characterized by thermal and solute effects as measured by a combination of effective Rayleigh numbers. Guba & Anderson (Reference Guba and Anderson2014) further explored these instabilities and identified other solute ‘imbalance’ mechanisms – for example, with respect to solute rejection and/or the initial solute compositions. These ‘double-solute-diffusion’ or ‘double-solute-rejection’ effects could lead to stationary and oscillatory instabilities that occur despite all diffusing agencies (heat and solutes) being stably stratified from a static point of view. It is the nonlinear development of these novel convective instabilities that we wish to explore in the present paper.
The plan of the paper is as follows. In § 2, the problem is formulated and the approximations that are employed are made explicit. In § 3, we revisit the linear stability problem. In § 4, we perform a weakly nonlinear stability analysis and explain the physical mechanism of the pattern selection. The predictions in terms of the relative stability between different convection patterns and global stability criteria are established. In § 5, we discuss the results and compare some of our predictions with the relevant experimental observations. We identify a particular aqueous ternary system as an appropriate candidate for potential experimental confirmation of steady convection with two stabilizing solutal gradients. The supplementary material (see https://doi.org/10.1017/jfm.2017.382) contains the detailed derivation of the weakly nonlinear solutions.
2 Formulation
2.1 Model equations
Model equations describing convection and directional solidification in a single (primary) ternary mush have been given by Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010) and later used by Guba & Anderson (Reference Guba and Anderson2014). We take them in the form
$$\begin{eqnarray}c^{\ast }\left(\frac{\unicode[STIX]{x2202}T^{\ast }}{\unicode[STIX]{x2202}t^{\ast }}-V^{\ast }\frac{\unicode[STIX]{x2202}T^{\ast }}{\unicode[STIX]{x2202}z^{\ast }}\right)+c^{\ast }\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\ast }T^{\ast }=k^{\ast }\unicode[STIX]{x1D6FB}^{\ast 2}T^{\ast }+L_{v}^{\ast }\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t^{\ast }}-V^{\ast }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z^{\ast }}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle (1-\unicode[STIX]{x1D719})\left(\frac{\unicode[STIX]{x2202}C_{j}^{\ast }}{\unicode[STIX]{x2202}t^{\ast }}-V^{\ast }\frac{\unicode[STIX]{x2202}C_{j}^{\ast }}{\unicode[STIX]{x2202}z^{\ast }}\right)+\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\ast }C_{j}^{\ast }=D_{j}^{\ast }\unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }[(1-\unicode[STIX]{x1D719})\unicode[STIX]{x1D735}^{\ast }C_{j}^{\ast }]\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,(1-k_{j})C_{j}^{\ast }\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t^{\ast }}-V^{\ast }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z^{\ast }}\right),\quad \text{for}\quad j=1,2,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}^{\ast }=-\frac{\unicode[STIX]{x1D6F1}^{\ast }(\unicode[STIX]{x1D719})}{\unicode[STIX]{x1D707}^{\ast }}(\unicode[STIX]{x1D735}^{\ast }p^{\ast }+\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }g^{\ast }\hat{\boldsymbol{k}}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }\boldsymbol{u}^{\ast }=0, & \displaystyle\end{eqnarray}$$
$T^{\ast }$
is the temperature,
$C_{1}^{\ast }$
and
$C_{2}^{\ast }$
are the liquid compositions,
$\unicode[STIX]{x1D719}$
is the solid fraction,
$p^{\ast }$
is the reduced pressure (i.e. the dynamic pressure plus a hydrostatic part,
$p^{\ast }=p_{dyn}+\unicode[STIX]{x1D70C}_{top}^{\ast }g^{\ast }z^{\ast }$
) and
$\boldsymbol{u}^{\ast }$
is the Darcy velocity. The material parameters are
$c^{\ast }$
, the specific heat;
$k^{\ast }$
, the thermal conductivity;
$L_{v}^{\ast }$
, the latent heat per unit volume;
$D_{j}^{\ast }$
, the solutal diffusivity in the liquid for species
$j$
;
$k_{j}$
, the segregation coefficient for species
$j$
;
$\unicode[STIX]{x1D6F1}^{\ast }(\unicode[STIX]{x1D719})$
, the permeability of the mush;
$\unicode[STIX]{x1D707}^{\ast }$
, the fluid dynamic viscosity;
$g^{\ast }$
, the gravitational acceleration, and
$\hat{\boldsymbol{k}}$
is the upward unit vector. For simplicity we have assumed that, with the exception of solute diffusion which is taken to be negligible in the solid phase, material properties in the solid and liquid phases are equal. In writing (2.1), we have decomposed the fluid motion into a portion due to the uniform downward translation of the system at a prescribed constant rate
$V^{\ast }$
relative to the fixed laboratory frame (an effect we refer to as background solidification) and a buoyancy-driven portion
$\boldsymbol{u}^{\ast }$
. A word of caution is in order at the form (2.1b
) of the solute conservation equations. It assumes that the solid phase of the mush, a ternary solid solution, is of sufficient ternary composition so as the release of solutes on local dissolution keeps the solid at internal interfaces in chemical equilibrium with the liquid. This assumption, though somewhat arbitrary and not particularly realistic, avoids a difficulty of history-dependent behaviour in the case of local dissolution (see also Fowler Reference Fowler1997).
In addition, we assume that local thermodynamic equilibrium prevails so that the liquid compositions
$C_{1}^{\ast }$
and
$C_{2}^{\ast }$
and the temperature are coupled by the liquidus constraint defined by
$$\begin{eqnarray}T^{\ast }=T_{top}^{\ast }+m_{1}^{\ast }(C_{1}^{\ast }-C_{1top}^{\ast })+m_{2}^{\ast }(C_{2}^{\ast }-C_{2top}^{\ast }),\end{eqnarray}$$
where
$m_{j}^{\ast }$
are the liquidus surface slopes, and
$T_{top}^{\ast }$
and
$C_{jtop}^{\ast }$
(
$j=1,2$
) are the reference temperature and liquid compositions taken to correspond to the top of the mush at
$z^{\ast }=H^{\ast }$
.
Convective motions in the mush are driven by a density differential
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}^{\ast }-\unicode[STIX]{x1D70C}_{top}^{\ast }$
in (2.1c
), with the density assumed to be linearly related to the temperature and compositions,
$$\begin{eqnarray}\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}_{top}^{\ast }[1-\unicode[STIX]{x1D6FC}^{\ast }(T^{\ast }-T_{top}^{\ast })-\unicode[STIX]{x1D6FC}_{1}^{\ast }(C_{1}^{\ast }-C_{1top}^{\ast })-\unicode[STIX]{x1D6FC}_{2}^{\ast }(C_{2}^{\ast }-C_{2top}^{\ast })],\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FC}^{\ast }$
and
$\unicode[STIX]{x1D6FC}_{j}^{\ast }$
(
$j=1,2$
) are coefficients of thermal and compositional expansion, respectively.
The boundary conditions at the top and bottom of the mush are
$$\begin{eqnarray}\displaystyle C_{1}^{\ast }=C_{1top}^{\ast },\quad C_{2}^{\ast }=C_{2top}^{\ast },\quad \unicode[STIX]{x1D719}=\unicode[STIX]{x1D6F7},\quad \boldsymbol{u}^{\ast }\boldsymbol{\cdot }\hat{\boldsymbol{k}}=0\quad \text{at}~z^{\ast }=H^{\ast }, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle C_{1}^{\ast }=C_{1bot}^{\ast },\quad C_{2}^{\ast }=C_{2bot}^{\ast },\quad \boldsymbol{u}^{\ast }\boldsymbol{\cdot }\hat{\boldsymbol{k}}=0\quad \text{at}~z^{\ast }=0, & & \displaystyle\end{eqnarray}$$
where
$C_{jtop}^{\ast }$
,
$C_{jbot}^{\ast }$
(
$j=1,2$
) are prescribed constant values of the concentrations at the top and bottom of the mushy layer and
$\unicode[STIX]{x1D6F7}$
is a prescribed constant value for the solid fraction at the top of the mushy layer. The thickness of the mushy layer
$H^{\ast }$
is assumed to be an input parameter rather than an unknown quantity to be determined. Note that, in light of (2.2), the temperature is automatically fixed at
$T_{top}^{\ast }$
at the mush top, while the temperature at the mush bottom attains a dependent value
$T_{bot}^{\ast }=T_{top}^{\ast }+m_{1}^{\ast }(C_{1bot}^{\ast }-C_{1top}^{\ast })+m_{2}^{\ast }(C_{2bot}^{\ast }-C_{2top}^{\ast })$
. These conditions, which were used in Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010) and Guba & Anderson (Reference Guba and Anderson2014) to mimic conditions at the boundaries of the primary mushy layer in a realistic setting, permit the analysis of the dynamics within the primary mush, free from direct coupling with any adjacent layers of the ternary alloy system.
2.2 Non-dimensionalization
We shall employ dimensional units that follow closely those used by Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010), with the exception of the temperature and compositions, which we rescale as in Guba & Anderson (Reference Guba and Anderson2014). Thus, we define
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{x}=\frac{\boldsymbol{x}^{\ast }}{H^{\ast }},\quad t=\frac{t^{\ast }}{H^{\ast 2}/\unicode[STIX]{x1D705}^{\ast }},\\ \displaystyle T=\frac{T^{\ast }-T_{top}^{\ast }}{\unicode[STIX]{x0394}T^{\ast }},\quad C_{j}=\frac{C_{j}^{\ast }-C_{jtop}^{\ast }}{\unicode[STIX]{x0394}C_{j}^{\ast }},\quad p=\frac{p^{\ast }}{\unicode[STIX]{x1D707}^{\ast }\unicode[STIX]{x1D705}^{\ast }/\unicode[STIX]{x1D6F1}^{\ast }(\unicode[STIX]{x1D6F7})},\quad \boldsymbol{u}=\frac{\boldsymbol{u}^{\ast }}{\unicode[STIX]{x1D705}^{\ast }/H^{\ast }},\end{array}\right\}\end{eqnarray}$$
where
$\unicode[STIX]{x1D705}^{\ast }=k^{\ast }/c^{\ast }$
is the thermal diffusivity,
$\unicode[STIX]{x0394}T^{\ast }=T_{top}^{\ast }-T_{bot}^{\ast }$
,
$\unicode[STIX]{x0394}C_{j}^{\ast }=C_{jtop}^{\ast }-C_{jbot}^{\ast }$
(
$j=1,2$
), and
$\unicode[STIX]{x1D6F1}^{\ast }(\unicode[STIX]{x1D6F7})$
is a reference permeability.
The dimensionless governing equations then have the form
$$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}-V\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\right)T+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T=\unicode[STIX]{x1D6FB}^{2}T+S\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}-V\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\right)\unicode[STIX]{x1D719}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle (1-\unicode[STIX]{x1D719})\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}-V\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\right)C_{j}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}C_{j}=\frac{1}{\mathit{Le}_{j}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }[(1-\unicode[STIX]{x1D719})\unicode[STIX]{x1D735}C_{j}]\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,(1-k_{j})(C_{j}^{top}+C_{j})\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}-V\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\right)\unicode[STIX]{x1D719},\quad \text{for}~j=1,2,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle T=m_{1}C_{1}+m_{2}C_{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle K(\unicode[STIX]{x1D719})\boldsymbol{u}=-\unicode[STIX]{x1D735}p-\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,\hat{\boldsymbol{k}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle C_{1}=0,\quad C_{2}=0,\quad \unicode[STIX]{x1D719}=\unicode[STIX]{x1D6F7},\quad \boldsymbol{u}\boldsymbol{\cdot }\hat{\boldsymbol{k}}=0\quad \text{at}~z=1, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle C_{1}=-1,\quad C_{2}=-1,\quad \boldsymbol{u}\boldsymbol{\cdot }\hat{\boldsymbol{k}}=0\quad \text{at}~z=0. & & \displaystyle\end{eqnarray}$$
Here
$K(\unicode[STIX]{x1D719})=\unicode[STIX]{x1D6F1}^{\ast }(\unicode[STIX]{x1D6F7})/\unicode[STIX]{x1D6F1}^{\ast }(\unicode[STIX]{x1D719})$
is a resistance function, and the dimensionless density differential is
$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=-\mathit{Ra}T-\mathit{Ra}_{1}C_{1}-\mathit{Ra}_{2}C_{2}.\end{eqnarray}$$
The non-dimensional groups that are formed as a result of scalings adopted are
$$\begin{eqnarray}\displaystyle V=\frac{V^{\ast }}{\unicode[STIX]{x1D705}^{\ast }/H^{\ast }},\quad S=\frac{L_{v}^{\ast }}{c^{\ast }\unicode[STIX]{x0394}T^{\ast }},\quad \mathit{Le}_{j}=\frac{\unicode[STIX]{x1D705}^{\ast }}{D_{j}^{\ast }},\quad C_{j}^{top}=\frac{C_{jtop}^{\ast }}{\unicode[STIX]{x0394}C_{j}^{\ast }},\quad m_{j}=\frac{m_{j}^{\ast }\unicode[STIX]{x0394}C_{j}^{\ast }}{\unicode[STIX]{x0394}T^{\ast }},\qquad & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \mathit{Ra}=\frac{\unicode[STIX]{x1D70C}_{top}^{\ast }\unicode[STIX]{x1D6FC}^{\ast }\unicode[STIX]{x0394}T^{\ast }g^{\ast }\unicode[STIX]{x1D6F1}^{\ast }(\unicode[STIX]{x1D6F7})H^{\ast }}{\unicode[STIX]{x1D707}^{\ast }\unicode[STIX]{x1D705}^{\ast }},\quad \mathit{Ra}_{j}=\frac{\unicode[STIX]{x1D70C}_{top}^{\ast }\unicode[STIX]{x1D6FC}_{j}^{\ast }\unicode[STIX]{x0394}C_{j}^{\ast }g^{\ast }\unicode[STIX]{x1D6F1}^{\ast }(\unicode[STIX]{x1D6F7})H^{\ast }}{\unicode[STIX]{x1D707}^{\ast }\unicode[STIX]{x1D705}^{\ast }}. & & \displaystyle\end{eqnarray}$$
As mentioned earlier, the system described above has been analysed analytically and numerically in linear stability settings for a wide range of parameter regimes by Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010) and Guba & Anderson (Reference Guba and Anderson2014). Of particular note were convective instabilities that could arise when not only was the density stably stratified in an overall sense (e.g. when accounting for its net dependence on temperature and two solutes) but also when the density was stably stratified with respect to each agency (i.e. no individual temperature or solute field promoted an unstable density gradient). These novel instabilities were found to be a very robust feature of this system and could be identified even when many of the interactions represented in (2.7) were neglected. In particular, Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010) showed that the qualitative nature of these unexpected instabilities was still retained even with the simplifications of no latent-heat release (
$S=0$
), no solute rejection (
$k_{j}=1$
) and no background (or overall) solidification (
$V=0$
). With respect to the assumption of no solute rejection in the solute balance (
$k_{j}=1$
) it is important to recognize that a concentration gradient is still imposed across the mushy layer through conditions on
$C_{j}$
at the boundaries. In an actual solidification system it is the presence of solute rejection that would generate a solute gradient across a mushy layer, and solute gradients are certainly one of the requirements for the mushy layer instabilities of interest here. We recognize that while some quantitative features of the instabilities change under these simplifying assumptions – for example in the simplified case above (
$S=V=k_{j}-1=0$
) two different Lewis numbers (
$Le_{1}\neq Le_{2}$
) are required for the appearance of the unexpected instabilities while in the more general case this condition is not required – necessary ingredients for the existence of the instabilities are still present. Thus, we proceed in the simplest known context that admits these unexpected instabilities. As demonstrated by Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010) and outlined briefly below, in this simplified setting the linear stability problem can be characterized with closed-form analytical expressions. This closed-form nature of the linear stability results allow many of our new weakly nonlinear results to be expressed in closed form; this enhances our ability to physically interpret the results.
2.3 A base state
We consider the situation in which
$V=0$
so that the background solidification as the system is pulled is ignored;
$S=0$
so that the effects of latent-heat release are omitted; and
$k_{j}=1$
so that there is no influence of rejection of either solute upon solidification in the local solute balances. Recalling that we fix the concentrations at the top and bottom of the mushy layer, a steady state in our model is produced by extracting or inputting solutes at the top and bottom of a pre-existing mush, with the steady-state compositional profiles maintained through the action of solute diffusion. With these assumptions, the equations given in § 2.2 admit a steady one-dimensional solution
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{T}=-(1-z), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{C}_{j}=-(1-z), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D6F7}. & \displaystyle\end{eqnarray}$$
$m_{1}+m_{2}=1$
, a relation which follows from the liquidus constraint (2.7c
) considered at the bottom of the mush. Note that in this base-state solution the temperature and solutal fields are linear, and the solid fraction throughout the mushy layer is constant. This base state has no flow and a hydrostatic pressure distribution.
2.4 Nonlinear disturbance equations
We shall seek nonlinear solutions to the system (2.7) by perturbing the base state, denoted by a tilde, with a convective perturbation, denoted by a caret, as
$$\begin{eqnarray}\displaystyle & \displaystyle T=\tilde{T}(z)+\unicode[STIX]{x1D716}\hat{T}(x,y,z,t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle C_{j}=\tilde{C}_{j}(z)+\unicode[STIX]{x1D716}{\hat{C}}_{j}(x,y,z,t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}=\tilde{\unicode[STIX]{x1D719}}(z)+\unicode[STIX]{x1D716}\hat{\unicode[STIX]{x1D719}}(x,y,z,t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}=\mathbf{0}+\unicode[STIX]{x1D716}\hat{\boldsymbol{u}}(x,y,z,t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p=\tilde{p}(z)+\unicode[STIX]{x1D716}\hat{p}(x,y,z,t), & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D716}$
is a small-but-finite amplitude,
$0<\unicode[STIX]{x1D716}\ll 1$
.
The disturbance equations, after eliminating the variable
$\hat{p}$
, are
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\hat{T}}{\unicode[STIX]{x2202}t}+{\hat{w}}-\unicode[STIX]{x1D6FB}^{2}\hat{T}=-\unicode[STIX]{x1D716}\hat{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\hat{T}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle (1-\unicode[STIX]{x1D6F7})\frac{\unicode[STIX]{x2202}{\hat{C}}_{j}}{\unicode[STIX]{x2202}t}+{\hat{w}}-\frac{1}{\mathit{Le}_{j}}\left[\left(1-\unicode[STIX]{x1D6F7}\right)\unicode[STIX]{x1D6FB}^{2}{\hat{C}}_{j}-\frac{\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x2202}z}\right]\nonumber\\ \displaystyle & & \displaystyle \qquad \quad =\unicode[STIX]{x1D716}\hat{\unicode[STIX]{x1D719}}\frac{\unicode[STIX]{x2202}{\hat{C}}_{j}}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D716}\hat{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\hat{C}}_{j}-\unicode[STIX]{x1D716}\frac{1}{\mathit{Le}_{j}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\hat{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D735}{\hat{C}}_{j}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{T}-m_{1}{\hat{C}}_{1}-m_{2}{\hat{C}}_{2}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\hat{u} +\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}(\mathit{Ra}\hat{T}+\mathit{Ra}_{1}{\hat{C}}_{1}+\mathit{Ra}_{2}{\hat{C}}_{2})=\unicode[STIX]{x1D716}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}[\hat{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{K}(\hat{\unicode[STIX]{x1D719}})]-\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FB}^{2}[\bar{K}(\hat{\unicode[STIX]{x1D719}})\hat{u} ], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\hat{v}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y\unicode[STIX]{x2202}z}(\mathit{Ra}\hat{T}+\mathit{Ra}_{1}{\hat{C}}_{1}+\mathit{Ra}_{2}{\hat{C}}_{2})=\unicode[STIX]{x1D716}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}[\hat{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{K}(\hat{\unicode[STIX]{x1D719}})]-\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FB}^{2}[\bar{K}(\hat{\unicode[STIX]{x1D719}})\hat{v}], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}{\hat{w}}-\unicode[STIX]{x1D6FB}_{2}^{2}(\mathit{Ra}\hat{T}+\mathit{Ra}_{1}{\hat{C}}_{1}+\mathit{Ra}_{2}{\hat{C}}_{2})=\unicode[STIX]{x1D716}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}[\hat{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{K}(\hat{\unicode[STIX]{x1D719}})]-\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FB}^{2}[\bar{K}(\hat{\unicode[STIX]{x1D719}}){\hat{w}}], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\hat{C}}_{1}=0,\quad {\hat{C}}_{2}=0,\quad \hat{\unicode[STIX]{x1D719}}=0,\quad {\hat{w}}=0\quad \text{at}~z=1, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\hat{C}}_{1}=0,\quad {\hat{C}}_{2}=0,\quad {\hat{w}}=0\quad \text{at}~z=0. & & \displaystyle\end{eqnarray}$$
The operator
$\unicode[STIX]{x1D6FB}_{2}^{2}$
denotes the horizontal Laplacian. In writing (2.15d
–
f
), we have assumed that the function
$K(\unicode[STIX]{x1D719})$
has a Taylor series about the (constant) base-state solid fraction
$\unicode[STIX]{x1D6F7}$
for small
$\unicode[STIX]{x1D716}$
,
$$\begin{eqnarray}K(\unicode[STIX]{x1D719})=1+\unicode[STIX]{x1D716}\bar{K}(\hat{\unicode[STIX]{x1D719}})\equiv 1+\unicode[STIX]{x1D716}K_{1}\hat{\unicode[STIX]{x1D719}}+\unicode[STIX]{x1D716}^{2}K_{2}\hat{\unicode[STIX]{x1D719}}^{2}+\cdots \,,\end{eqnarray}$$
where the coefficients
$K_{1},K_{2},\ldots$
represent the sensitivity of the permeability to changes in solid fraction. A particular form of the permeability function is then given by the choice of these coefficients.
An equivalent form of the disturbance equations (2.15), which is more convenient for the stability analysis, is
$$\begin{eqnarray}\displaystyle (1-\unicode[STIX]{x1D6F7})\left[m_{1}(\mathit{Le}_{1}-1)\frac{\unicode[STIX]{x2202}{\hat{C}}_{1}}{\unicode[STIX]{x2202}t}+m_{2}(\mathit{Le}_{2}-1)\frac{\unicode[STIX]{x2202}{\hat{C}}_{2}}{\unicode[STIX]{x2202}t}\right]+\frac{\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x2202}z}+\frac{1}{\unicode[STIX]{x1D6E4}}{\hat{w}}=\hat{\unicode[STIX]{x1D6F7}}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1-\unicode[STIX]{x1D6F7}}{\mathit{Le}_{1}}\left[[m_{1}(\mathit{Le}_{1}-1)-\mathit{Le}_{1}]\frac{\unicode[STIX]{x2202}{\hat{C}}_{1}}{\unicode[STIX]{x2202}t}+m_{2}(\mathit{Le}_{2}-1)\frac{\unicode[STIX]{x2202}{\hat{C}}_{2}}{\unicode[STIX]{x2202}t}\right]\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,\frac{1-\unicode[STIX]{x1D6F7}}{\mathit{Le}_{1}}\unicode[STIX]{x1D6FB}^{2}{\hat{C}}_{1}-\frac{(1-\unicode[STIX]{x1D6F7})\unicode[STIX]{x1D6FA}_{1}}{\mathit{Le}_{1}}{\hat{w}}=\hat{\unicode[STIX]{x1D6E9}}_{1},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1-\unicode[STIX]{x1D6F7}}{\mathit{Le}_{2}}\left[[m_{2}(\mathit{Le}_{2}-1)-\mathit{Le}_{2}]\frac{\unicode[STIX]{x2202}{\hat{C}}_{2}}{\unicode[STIX]{x2202}t}+m_{1}(\mathit{Le}_{1}-1)\frac{\unicode[STIX]{x2202}{\hat{C}}_{1}}{\unicode[STIX]{x2202}t}\right]\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,\frac{1-\unicode[STIX]{x1D6F7}}{\mathit{Le}_{2}}\unicode[STIX]{x1D6FB}^{2}{\hat{C}}_{2}-\frac{(1-\unicode[STIX]{x1D6F7})\unicode[STIX]{x1D6FA}_{2}}{\mathit{Le}_{2}}{\hat{w}}=\hat{\unicode[STIX]{x1D6E9}}_{2},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\hat{u} +\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}(\mathit{Ra}_{C1}{\hat{C}}_{1}+\mathit{Ra}_{C2}{\hat{C}}_{2})=\hat{U} , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\hat{v}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y\unicode[STIX]{x2202}z}(\mathit{Ra}_{C1}{\hat{C}}_{1}+\mathit{Ra}_{C2}{\hat{C}}_{2})=\hat{V}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}{\hat{w}}-\unicode[STIX]{x1D6FB}_{2}^{2}(\mathit{Ra}_{C1}{\hat{C}}_{1}+\mathit{Ra}_{C2}{\hat{C}}_{2})={\hat{W}}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\mathit{Ra}_{Cj}=m_{j}\mathit{Ra}+\mathit{Ra}_{j},\quad \text{for}~j=1,2.\end{eqnarray}$$
Note that
$\mathit{Ra}_{C1}+\mathit{Ra}_{C2}=\mathit{Ra}+\mathit{Ra}_{1}+\mathit{Ra}_{2}$
. To condense the notation a little, we have also introduced the parameter groups
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}_{1}=1+\frac{m_{2}\left(\mathit{Le}_{1}-\mathit{Le}_{2}\right)}{1-\unicode[STIX]{x1D6F7}},\quad \unicode[STIX]{x1D6FA}_{2}=1+\frac{m_{1}\left(\mathit{Le}_{2}-\mathit{Le}_{1}\right)}{1-\unicode[STIX]{x1D6F7}}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\unicode[STIX]{x1D6E4}=1/\left[m_{1}\mathit{Le}_{1}+m_{2}\mathit{Le}_{2}-(1-\unicode[STIX]{x1D6F7})\right].\end{eqnarray}$$
Equation (2.19a
) for the solid-fraction perturbation
$\hat{\unicode[STIX]{x1D719}}$
has been obtained by first combining the two solutal balances (2.15b
) (namely, multiplying (2.15b
) for
$j=1$
by
$m_{1}\mathit{Le}_{1}/(1-\unicode[STIX]{x1D6F7})$
, equation (2.15b
) for
$j=2$
by
$m_{2}\mathit{Le}_{2}/(1-\unicode[STIX]{x1D6F7})$
, and adding the two), and then subtracting the result from the heat balance (2.15a
) with
$\hat{T}$
eliminated by the liquidus constraint (2.15c
). A modified solute balance (2.19b
) for
${\hat{C}}_{1}$
, and similarly (2.19c
) for
${\hat{C}}_{2}$
, have been obtained by first eliminating
$\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D719}}/\unicode[STIX]{x2202}z$
between the two solutal balances (2.15b
) (namely, multiplying (2.15b
) for
$j=1$
by
$\mathit{Le}_{1}$
, equation (2.15b
) for
$j=2$
by
$\mathit{Le}_{2}$
, and subtracting the two), and then solving the resulting equation and the heat balance (2.15a
), with
$\hat{T}$
eliminated using (2.15c
), for the Laplacians of
${\hat{C}}_{1}$
and
${\hat{C}}_{2}$
. Equations (2.19d
–
f
) for the flow-field components
$\hat{u}$
,
$\hat{v}$
and
${\hat{w}}$
follow from (2.15d
–
f
) on elimination of
$\hat{T}$
using (2.15c
).
The parameter group
$\unicode[STIX]{x1D6FA}_{j}$
measures a coupling between the modified heat equation for
$\hat{\unicode[STIX]{x1D719}}$
and the modified solute balance for
${\hat{C}}_{j}$
. For the typical conditions that require
$\mathit{Le}_{j}>1$
, noting that
$m_{1}+m_{2}=1$
and
$0<\unicode[STIX]{x1D6F7}<1$
, we find that
$\unicode[STIX]{x1D6FA}_{j}$
can take either sign, while
$\unicode[STIX]{x1D6E4}$
is always positive. Note that
$\unicode[STIX]{x1D6FA}_{1}$
and
$\unicode[STIX]{x1D6FA}_{2}$
are related by
$m_{1}\unicode[STIX]{x1D6FA}_{1}+m_{2}\unicode[STIX]{x1D6FA}_{2}=1$
, implying that the two cannot be simultaneously negative. The groups
$\unicode[STIX]{x1D6FA}_{1}$
and
$\unicode[STIX]{x1D6FA}_{2}$
emerge as the key parameters determining the type of convective instability of the static diffusive state, as highlighted by the linear stability results of § 3, and embody in a combined sense the influence of solute and heat diffusion, the liquidus slopes, concentration and temperature differences across the mushy layer and the base-state solid phase volume fraction
$\unicode[STIX]{x1D6F7}$
. One advantage of (2.19) is that the symmetry between indices 1 and 2 carries on to the solvability conditions for the nonlinear field contributions to exist, providing some check on the weakly nonlinear calculations.
We expand the perturbation quantities in powers of
$\unicode[STIX]{x1D716}$
as
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{T}=T^{0}+\unicode[STIX]{x1D716}T^{1}+\unicode[STIX]{x1D716}^{2}T^{2}+\cdots \,, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\hat{C}}_{j}=C_{j}^{0}+\unicode[STIX]{x1D716}C_{j}^{1}+\unicode[STIX]{x1D716}^{2}C_{j}^{2}+\cdots \,, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D719}^{0}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D719}^{1}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D719}^{2}+\cdots \,, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{\boldsymbol{u}}=\boldsymbol{u}^{0}+\unicode[STIX]{x1D716}\boldsymbol{u}^{1}+\unicode[STIX]{x1D716}^{2}\boldsymbol{u}^{2}+\cdots \,, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}_{Cj}=\mathit{Ra}_{Cj}^{0}+\unicode[STIX]{x1D716}\mathit{Ra}_{Cj}^{1}+\unicode[STIX]{x1D716}^{2}\mathit{Ra}_{Cj}^{2}+\cdots \,. & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D716}$
, we obtain a sequence of linear equations to be solved in series. The details of these equations are given in the supplementary material (appendix B). The expansion of both
$\mathit{Ra}_{C1}$
and
$\mathit{Ra}_{C2}$
in (2.23e
) rationalizes a freedom in the choice between the two as a bifurcation parameter.
3 Linear stability theory
At
$O(\unicode[STIX]{x1D716}^{0})$
we recover the linear stability problem studied by Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010). The eigensolution in terms of a normal mode takes the form
$$\begin{eqnarray}\left[\begin{array}{@{}c@{}}C_{j}^{0}\\ w^{0}\end{array}\right]=\left[\begin{array}{@{}c@{}}\tilde{C}_{j}^{0}\\ \tilde{w}^{0}\end{array}\right]\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}}\text{e}^{\unicode[STIX]{x1D70E}t}\sin n\unicode[STIX]{x03C0}z+\text{c.c.},\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{C}_{1}^{0}=\frac{J\unicode[STIX]{x1D6FA}_{1}+\unicode[STIX]{x1D70E}[\mathit{Le}_{2}+m_{2}(\unicode[STIX]{x1D6FA}_{1}-\unicode[STIX]{x1D6FA}_{2})]}{J+\unicode[STIX]{x1D70E}(m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1})}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{C}_{2}^{0}=\frac{J\unicode[STIX]{x1D6FA}_{2}+\unicode[STIX]{x1D70E}[\mathit{Le}_{1}+m_{1}(\unicode[STIX]{x1D6FA}_{2}-\unicode[STIX]{x1D6FA}_{1})]}{J+\unicode[STIX]{x1D70E}(m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1})}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{w}^{0}=-(J+\unicode[STIX]{x1D70E}). & \displaystyle\end{eqnarray}$$
$\boldsymbol{r}$
is the horizontal position vector,
$\boldsymbol{k}$
is a single horizontal wavevector of arbitrary orientation,
$n\unicode[STIX]{x03C0}$
(
$n$
is an integer) is the vertical wavenumber and
$J=n^{2}\unicode[STIX]{x03C0}^{2}+k^{2}$
with
$k=|\boldsymbol{k}|$
. The associated dispersion relation for the complex growth rate
$\unicode[STIX]{x1D70E}$
is 
$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle [J+\unicode[STIX]{x1D70E}(m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1})]J(J+\unicode[STIX]{x1D70E})\nonumber\\ \displaystyle & & \displaystyle +\,k^{2}\mathit{Ra}_{C1}^{0}\{J\unicode[STIX]{x1D6FA}_{1}+\unicode[STIX]{x1D70E}[\mathit{Le}_{2}+m_{2}(\unicode[STIX]{x1D6FA}_{1}-\unicode[STIX]{x1D6FA}_{2})]\}\nonumber\\ \displaystyle & & \displaystyle +\,k^{2}\mathit{Ra}_{C2}^{0}\left\{J\unicode[STIX]{x1D6FA}_{2}+\unicode[STIX]{x1D70E}[\mathit{Le}_{1}+m_{1}(\unicode[STIX]{x1D6FA}_{2}-\unicode[STIX]{x1D6FA}_{1})]\right\}.\end{eqnarray}$$
The remaining equations are then solved to yield
$$\begin{eqnarray}\displaystyle & \displaystyle u^{0}=-\frac{n\unicode[STIX]{x03C0}\left(J+\unicode[STIX]{x1D70E}\right)}{k^{2}}\frac{\unicode[STIX]{x2202}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}}}{\unicode[STIX]{x2202}x}\text{e}^{\unicode[STIX]{x1D70E}t}\cos n\unicode[STIX]{x03C0}z+\text{c.c.}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle v^{0}=-\frac{n\unicode[STIX]{x03C0}\left(J+\unicode[STIX]{x1D70E}\right)}{k^{2}}\frac{\unicode[STIX]{x2202}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}}}{\unicode[STIX]{x2202}y}\text{e}^{\unicode[STIX]{x1D70E}t}\cos n\unicode[STIX]{x03C0}z+\text{c.c.}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}^{0} & = & \displaystyle \left[\frac{J}{\unicode[STIX]{x1D6E4}}+\unicode[STIX]{x1D70E}\right.\nonumber\\ \displaystyle & & \displaystyle \left.\times \,\left(\frac{1}{\unicode[STIX]{x1D6E4}}-\frac{(1-\unicode[STIX]{x1D6F7})[m_{1}\unicode[STIX]{x1D6FA}_{1}(\mathit{Le}_{1}-1)(J+\unicode[STIX]{x1D70E}\mathit{Le}_{2})+m_{2}\unicode[STIX]{x1D6FA}_{2}(\mathit{Le}_{2}-1)(J+\unicode[STIX]{x1D70E}\mathit{Le}_{1})]}{J+\unicode[STIX]{x1D70E}(m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1})}\right)\right]\nonumber\\ \displaystyle & & \displaystyle \times \,\frac{(-1)^{n}-\cos n\unicode[STIX]{x03C0}z}{n\unicode[STIX]{x03C0}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}}\text{e}^{\unicode[STIX]{x1D70E}t}+\text{c.c}.\end{eqnarray}$$
For neutrally stable real modes,
$\unicode[STIX]{x1D70E}=0$
, equation (3.3) gives
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{1}\mathit{Ra}_{C1}^{0}+\unicode[STIX]{x1D6FA}_{2}\mathit{Ra}_{C2}^{0}=-\frac{J^{2}}{k^{2}}=-\frac{(n^{2}\unicode[STIX]{x03C0}^{2}+k^{2})^{2}}{k^{2}}.\end{eqnarray}$$
For neutrally stable oscillatory modes,
$\unicode[STIX]{x1D70E}=0+\text{i}\unicode[STIX]{x1D714}$
, the real and imaginary parts of (3.3) provide the stability boundary
$$\begin{eqnarray}\left[1+\frac{m_{2}\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D6FA}_{1}-\unicode[STIX]{x1D6FA}_{2})}{m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1}}\right]\mathit{Ra}_{C1}^{0}+\left[1+\frac{m_{1}\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D6FA}_{2}-\unicode[STIX]{x1D6FA}_{1})}{m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1}}\right]\mathit{Ra}_{C2}^{0}=-\frac{J^{2}}{k^{2}}\left[1+\frac{1}{m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1}}\right]\end{eqnarray}$$
and a dispersion relation for the frequency of the oscillation given by
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}=\frac{k^{2}}{m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1}}\left[\frac{J^{2}}{k^{2}}+\unicode[STIX]{x1D6FA}_{1}\mathit{Ra}_{C1}^{0}+\unicode[STIX]{x1D6FA}_{2}\mathit{Ra}_{C2}^{0}\right],\end{eqnarray}$$
provided
$\unicode[STIX]{x1D714}^{2}>0$
is insured.
Figure 1 shows a set of neutral stability curves in the
$(k,\mathit{Ra}_{C1})$
-plane for representative values (a)
$\mathit{Ra}_{C2}^{0}<\mathit{Ra}_{C2}^{\dagger }$
, (b)
$\mathit{Ra}_{C2}^{0}=\mathit{Ra}_{C2}^{\dagger }$
, (c)
$\mathit{Ra}_{C2}^{\dagger }<\mathit{Ra}_{C2}^{0}<\mathit{Ra}_{C2}^{\ast }$
, (d)
$\mathit{Ra}_{C2}^{0}=\mathit{Ra}_{C2}^{\ast }$
and (e)
$\mathit{Ra}_{C2}^{0}>\mathit{Ra}_{C2}^{\ast }$
, where
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}_{C2}^{\dagger }=-\frac{4\unicode[STIX]{x03C0}^{2}m_{2}(1-\unicode[STIX]{x1D6F7})(\mathit{Le}_{2}-1)}{m_{2}(\mathit{Le}_{1}-\mathit{Le}_{2})+(1-\unicode[STIX]{x1D6F7})\mathit{Le}_{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}_{C2}^{\ast }=-\frac{4\unicode[STIX]{x03C0}^{2}\left[1-\unicode[STIX]{x1D6F7}+m_{2}(\mathit{Le}_{1}-\mathit{Le}_{2})(m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1}+1-\unicode[STIX]{x1D6F7})\right]}{(\mathit{Le}_{1}-\mathit{Le}_{2})(m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1}-\unicode[STIX]{x1D6F7})}. & \displaystyle\end{eqnarray}$$
$n=1$
and the other parameters remain fixed so as this figure applies to the case
$(1-\unicode[STIX]{x1D6F7})/m_{2}<\mathit{Le}_{2}-\mathit{Le}_{1}$
(cf. figure 2
a where
$\unicode[STIX]{x1D6FA}_{1}<0$
). In figure 1, each of the five sketches shows the neutral stability curves for the direct (solid) and oscillatory (dashed) modes, and a transition boundary separating the growing direct modes from the growing oscillatory modes (dotted). The regions where the system is dynamically unstable to linear perturbations are shaded. In figure 1(a) the oscillatory and transition boundaries each have two disjoint branches, with terminations on the direct branch where the frequency vanishes. There is a window of wavenumbers in which the neutral oscillatory mode does not exist and where the system is unstable for any value of
$\mathit{Ra}_{C1}$
. In figure 1(b), the two branches of the transition boundary intersect and then for yet larger
$\mathit{Ra}_{C2}$
the curves are transformed as shown in figure 1(c). In figure 1(d), the oscillatory mode meets the direct mode at a single wavenumber
$k=\unicode[STIX]{x03C0}$
, while the upper branch of the transition boundary ceases to exist. In figure 1(e), the neutral stability and transition boundaries are all disjoint, with the regions of instability above the direct branch and below the oscillatory branch separated by a region of linear stability. A feature illustrated in figure 1(e) is that an instability is possible for either sign of
$\mathit{Ra}_{C1}$
. The physical significance of this result is discussed below.
Figure 1. Schematic diagrams of the neutral stability curves in the (
$k$
,
$\mathit{Ra}_{C1}$
)-plane for (a)
$\mathit{Ra}_{C2}<\mathit{Ra}_{C2}^{\dagger }$
, (b)
$\mathit{Ra}_{C2}=\mathit{Ra}_{C2}^{\dagger }$
, (c)
$\mathit{Ra}_{C2}^{\dagger }<\mathit{Ra}_{C2}<\mathit{Ra}_{C2}^{\ast }$
, (d)
$\mathit{Ra}_{C2}=\mathit{Ra}_{C2}^{\ast }$
and (e)
$\mathit{Ra}_{C2}>\mathit{Ra}_{C2}^{\ast }$
, where
$\mathit{Ra}_{C2}^{\dagger }$
and
$\mathit{Ra}_{C2}^{\ast }$
are defined by (3.8). In each plot, the solid (dashed) curve corresponds to the direct (oscillatory) modes, while the dotted curve indicates the transition boundary separating the growing direct modes from the growing oscillatory modes. The regions of linear instability are shaded. The insets in (a) show the locations of the two normal-mode growth rates in the complex plane. Diagram (b) represents the transition between (a) and (c), while diagram (d) represents the transition between (c) and (e). Note that the extrema along the real branch and the disconnected oscillatory branch occur always for
$n=1$
and
$k=\unicode[STIX]{x03C0}$
.
Figure 2. Schematic diagrams showing sketches of the linear stability regimes in the
$(\mathit{Ra}_{C1},\mathit{Ra}_{C2})$
-plane for representative cases (a)
$\unicode[STIX]{x1D6FA}_{1}<0<\unicode[STIX]{x1D6FA}_{2}$
, (b)
$0<\unicode[STIX]{x1D6FA}_{1}<\unicode[STIX]{x1D6FA}_{2}$
, (c)
$\unicode[STIX]{x1D6FA}_{1}=\unicode[STIX]{x1D6FA}_{2}=1$
, (d)
$0<\unicode[STIX]{x1D6FA}_{2}<\unicode[STIX]{x1D6FA}_{1}$
and (e)
$\unicode[STIX]{x1D6FA}_{2}<0<\unicode[STIX]{x1D6FA}_{1}$
. The red solid curve indicates the critical linear stability boundary for the neutrally stable direct modes, the red dashed curve indicates the critical linear stability boundary for the neutrally stable oscillatory modes, and the red dotted curve represents the critical transition boundary between the growing direct and oscillatory modes. The yellow-coloured region indicates where the base-state solution is net statically stable. The cyan-coloured region indicates where the base-state solution is individually statically stable. The system is linearly unstable in the grey-shaded region. The insets in (a) show the locations of the growth rates in the complex plane. Horizontal sections at the points marked ● and ○ in (a) correspond to vertical sections at
$k=\unicode[STIX]{x03C0}$
in figure 1(b,d) respectively.
A set of qualitatively different representations of the linear stability boundaries in the
$(\mathit{Ra}_{C1},\mathit{Ra}_{C2})$
-plane are sketched in figure 2. These are categorized as follows: (a)
$(1-\unicode[STIX]{x1D6F7})/m_{2}<\mathit{Le}_{2}-\mathit{Le}_{1}$
, (b)
$0<\mathit{Le}_{2}-\mathit{Le}_{1}<(1-\unicode[STIX]{x1D6F7})/m_{2}$
, (c)
$\mathit{Le}_{1}=\mathit{Le}_{2}$
, (d)
$0<\mathit{Le}_{1}-\mathit{Le}_{2}<(1-\unicode[STIX]{x1D6F7})/m_{1}$
and (e)
$(1-\unicode[STIX]{x1D6F7})/m_{1}<\mathit{Le}_{1}-\mathit{Le}_{2}$
. In light of the significance of parameter groups
$\unicode[STIX]{x1D6FA}_{1}$
and
$\unicode[STIX]{x1D6FA}_{2}$
, these cases can be equivalently represented as (a)
$\unicode[STIX]{x1D6FA}_{1}<0<\unicode[STIX]{x1D6FA}_{2}$
, (b)
$0<\unicode[STIX]{x1D6FA}_{1}<\unicode[STIX]{x1D6FA}_{2}$
, (c)
$\unicode[STIX]{x1D6FA}_{1}=\unicode[STIX]{x1D6FA}_{2}=1$
, (d)
$0<\unicode[STIX]{x1D6FA}_{2}<\unicode[STIX]{x1D6FA}_{1}$
and (e)
$\unicode[STIX]{x1D6FA}_{2}<0<\unicode[STIX]{x1D6FA}_{1}$
. The results shown correspond to the most dangerous modes that always have
$n=1$
and
$k=\unicode[STIX]{x03C0}$
at the critical points of the neutral stability curves. As a reference to interpret the nonlinear results below, we identify the regions where the system is stable from a static viewpoint: the yellow-coloured region indicates the static stability when the net base-state density gradient is stabilizing; the cyan-coloured region indicates the static stability when all the individual base-state density-gradient contributions (thermal and the two solutal) are stabilizing. We pay particular attention to the prediction of instabilities that occur in the cyan-coloured region since, in that region, no convective instability would be expected from a static density stratification point of view. Instabilities in these regions are what we refer to as ‘novel’ or ‘unexpected’. Examples of such instabilities can be seen in figure 2(a,e).
While the linear stability results are most conveniently represented in terms of the effective Rayleigh numbers
$\mathit{Ra}_{C1}$
and
$\mathit{Ra}_{C2}$
, the physical understanding of convection driven by double-diffusive phase-change effects can only be gained through the natural Rayleigh numbers
$\mathit{Ra}$
,
$\mathit{Ra}_{1}$
and
$\mathit{Ra}_{2}$
. For the purposes of interpretation and for definiteness in the remainder of this section we assume that the interdendritic liquid has the solute-
$j$
composition at the mush top less than at the mush bottom (
$\unicode[STIX]{x0394}C_{j}^{\ast }=C_{jtop}^{\ast }-C_{jbot}^{\ast }<0$
). With the convention that the density increases (decreases) with concentrations when
$\unicode[STIX]{x1D6FC}_{j}^{\ast }<0$
(
$\unicode[STIX]{x1D6FC}_{j}^{\ast }>0$
),
$\mathit{Ra}_{j}$
positive (negative) denotes a solutally stable (unstable) base-state density gradient.
The diagrams shown in figure 2(a–e) are consistent with those in Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010). Here we have drawn these in terms of the two effective Rayleigh numbers
$\mathit{Ra}_{C1}$
and
$\mathit{Ra}_{C2}$
and, as noted above, have indicated the associated growth rates as well as the transition curve separating direct and oscillatory instabilities. In what follows, we revisit the parcel arguments put forward in Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010) with an additional connection to the reduction in the solute diffusivity associated with freezing. Specifically, while
$D_{j}^{\ast }$
are the actual diffusion coefficients for species
$j$
in the liquid phase,
$(1-\unicode[STIX]{x1D719})D_{j}^{\ast }$
is the effective diffusion coefficient in the porous mush which can be increased by melting or decreased by freezing (e.g. by phase fraction perturbation variations).
We follow simplifications outlined in Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010) and consider a system in which buoyancy may be present only with respect to composition 1. That is, we shall assume
$\mathit{Ra}=\mathit{Ra}_{2}=0$
while
$\mathit{Ra}_{1}\neq 0$
. In terms of the diagrams in figure 2 this situation corresponds to the horizontal line
$\mathit{Ra}_{C2}=0$
. In the following three paragraphs we shall outline three parcel arguments that correspond to figure 2(a) (
$\mathit{Le}_{2}>\mathit{Le}_{1}$
), figure 2(c) (in which
$\mathit{Le}_{1}=\mathit{Le}_{2}$
) and figure 2(e) (
$\mathit{Le}_{1}>\mathit{Le}_{2}$
).
As noted earlier, a novel feature of figure 2(a) is the appearance of a direct mode of instability that occurs in the region where the base-state temperature and composition fields are each individually stably stratified (in the cyan-coloured region). In terms of the special case with
$\mathit{Ra}=\mathit{Ra}_{2}=0$
this instability occurs at a positive value of
$\mathit{Ra}_{1}$
. A parcel argument in this case, adopted from Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010), can be outlined as follows. Consider the situation in which, in terms of diffusion rates, solute 2 is the slowest and heat is the fastest, with solute 1 diffusing at an intermediate rate. Additionally, the layer is cold and rich with respect to both solutes at the bottom and warmer and fresher at the top. When a parcel is displaced upwards, the parcel warms relatively rapidly while remaining at a nearly unchanged concentration of solute 2. In order to maintain thermodynamic equilibrium – that is, for the temperature and compositions to remain on the liquidus surface – the parcel’s composition 1 must overcompensate and hence take on a composition value that is less than its surroundings. When the density increases with increasing composition 1 (i.e.
$\unicode[STIX]{x1D6FC}_{1}^{\ast }<0$
) this means that the parcel will be lighter than its surroundings and will continue to rise. This physical description is consistent with a direct mode of instability with
$\mathit{Ra}_{1}>0$
, as is indicated in figure 2(a). As noted above, this is the mechanism reported and described in Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010) and requires that the diffusion rates of solute are different from the thermal diffusion rate and, in this instance, the solute 2 diffusion rate is smaller than that for solute 1. We note that under the set of simplifying assumptions made here, namely
$V=0$
and
$k_{j}=1$
, the solute balance equations (2.7b
) show that the compositional adjustments, particularly associated with solute 1, are achieved through a disparity in the rates of solute diffusion tied to the assumption
$Le_{2}>Le_{1}$
. We make the additional note here that solidification plays a role in amplifying the difference between solute and thermal diffusion. That is, local freezing reduces the effective solute diffusion coefficients (in an equal way with respect to both solute 1 and solute 2) by a factor of the liquid fraction,
$1-\unicode[STIX]{x1D719}$
. This follows from an underlying assumption in the model that species diffusion occurs only in the liquid phase and not in the solid phase.
It is instructive at this point to now consider the situation described in figure 2(c) in which the two solute diffusion rates are the same (i.e.
$\mathit{Le}_{1}=\mathit{Le}_{2}$
). Here if we again focus on the case with
$\mathit{Ra}=\mathit{Ra}_{2}=0$
and
$\mathit{Ra}_{1}\neq 0$
we observe that direct modes of instability are predicted to occur only for the case with
$\mathit{Ra}_{1}<0$
. In the case with cold and rich fluid at the bottom of the layer and warmer and fresher fluid at the top of the layer, we see that only when
$\unicode[STIX]{x1D6FC}_{1}^{\ast }>0$
(i.e. density decreases with increasing composition 1) does an instability arise. This would be a statically expected mode of instability in that the layer is top-heavy. In this case, a parcel displaced upwards equilibrates relatively rapidly with respect to temperature and relatively slowly with respect to composition (and in an unbiased way as
$\mathit{Le}_{1}=\mathit{Le}_{2}$
). A parcel that starts out relatively rich in terms of its composition 1 remains relatively rich. Then, in the case
$\unicode[STIX]{x1D6FC}_{1}^{\ast }>0$
(i.e.
$\mathit{Ra}_{1}<0$
since
$\unicode[STIX]{x0394}C_{1}^{\ast }<0$
) the parcel would find itself lighter than its surroundings and would continue to rise. Again, local solidification would reduce the effective solute diffusivities
$(1-\unicode[STIX]{x1D719})D_{j}^{\ast }$
relative to the thermal diffusivity and further promote the direct instability.
The situation in figure 2(e) with
$\mathit{Ra}=\mathit{Ra}_{2}=0$
and
$\mathit{Ra}_{1}\neq 0$
parallels the case described in the previous paragraph for figure 2(c) in the sense that a statically expected mode of instability would be predicted to occur with
$\mathit{Ra}_{1}<0$
. Unlike figure 2(c), this figure now captures the situation when
$\mathit{Le}_{1}>\mathit{Le}_{2}$
. Therefore, the sluggishness of solute 1 that was a key element of the parcel argument in the previous paragraph, is now even more pronounced. It is worth pointing out here that the sluggishness of solute 1 would suggest an overcompensation of solute 2, and indeed we see a direct mode of instability that can occur in the region of individual static stability (cyan-coloured region) for the case
$\mathit{Ra}=\mathit{Ra}_{1}$
and
$\mathit{Ra}_{2}>0$
. This latter mode (with
$\mathit{Ra}_{2}>0$
) is physically equivalent to the one discussed above for figure 2(a); only the indices 1 and 2 are reversed.
The two convective regimes of interest are distinguished also by the structure of the corresponding eigenfunctions. For the special case
$\mathit{Ra}=\mathit{Ra}_{2}=0$
while
$\mathit{Ra}_{1}\neq 0$
, figure 3 illustrates the spatial structure of the neutrally stable direct modes with
$\boldsymbol{k}=(\unicode[STIX]{x03C0},0)$
for figure 3(a)
$\unicode[STIX]{x1D6FA}_{1}>0$
(a mode with
$\mathit{Ra}_{1}<0$
on the horizontal axis in figure 2
e) and figure 3(b)
$\unicode[STIX]{x1D6FA}_{1}<0$
(a mode with
$\mathit{Ra}_{1}>0$
on the horizontal axis in figure 2
a). In both cases, the flow field takes the form of a well-known roll-like pattern. In figure 3(a), the upward vertical flow is associated with positive perturbations to the solute-1 composition (solute-1 enriched), which is the expected correlation for the standard compositional buoyancy argument (recall that only solute 1 contributes to buoyancy in this simple case). In contrast, in figure 3(b), the phase relation is just the opposite: the rising fluid is associated with negative compositional anomalies (solute-1 depleted). Note that this convective regime is consistent with our parcel argument for figure 2(a) discussed above. Curiously, in both figures 3(a) and 3(b), there is a positive correlation between the vertical flow and the solid-fraction perturbation, leading to vertically aligned channels of reduced solid fraction where the flow is downwards. If the permeability is allowed to vary with solid fraction, such an unusual phase relation is suggestive of the formation of ‘inverse’ chimneys hidden within the mush. This unusual phase relation is apparently a consequence of the zero-pulling-speed approximation (
$V=0$
), which excludes the effects of heat and frame advection in the present ternary mush model, and perhaps warrants further investigation in an analysis where these assumptions are not applied. This should be contrasted with the results obtained for the binary (Anderson & Worster Reference Anderson and Worster1995, Reference Anderson and Worster1996) and ternary (Anderson et al.
Reference Anderson, McFadden, Coriell and Murray2010, full numerical model; Guba & Anderson Reference Guba and Anderson2014) mushy layer systems incorporating the solidification caused by overall cooling.
Figure 3. Structure of the neutrally stable steady convective modes for
$\mathit{Ra}=0$
,
$\mathit{Ra}_{1}\neq 0$
,
$\mathit{Ra}_{2}=0$
and (a)
$\unicode[STIX]{x1D6FA}_{1}>0$
and (b)
$\unicode[STIX]{x1D6FA}_{1}<0$
. The onset of linear instability is statically expected in (a) and statically unexpected in (b). Each plot illustrates the eigensolution in terms of the convective streamlines (solid) and the contours of the total solute-1 composition (dashed). Also indicated are the vertical channels of minimum solid fraction in the mush and the compositional stripes in the resulting ternary eutectic solid (dotted), assuming no secondary mush formation. A notable difference between cases (a) and (b) is the phase relation between the solute-1 composition and the flow field. Rather surprisingly, the solid-fraction channels occur where the flow is downwards in both cases.
It is the nonlinear development of the direct instability which we shall elucidate here. We shall restrict the analysis to the parameter regime in which the oscillatory instability does not arise so as to avoid its possible interaction with the direct mode near onset.
4 Weakly nonlinear analysis
In this section we report on the results of a weakly nonlinear analysis extending the above-described linear stability results. The calculations were performed symbolically using Mathematica. The presentation of this section is given in terms of parameter groups that arise naturally in the calculations (e.g. parameter groups
$\unicode[STIX]{x1D6FA}_{j}$
and effective Rayleigh numbers
$\mathit{Ra}_{Cj}$
). In the section that follows we revisit these results with the goal of providing physically based interpretations using more physically tangible parameters.
We restrict our investigation to two cases based on two and three interacting convection rolls. In the first case, the motion is the superposition of two perpendicular rolls with
$$\begin{eqnarray}\boldsymbol{k}_{1}=(\unicode[STIX]{x03C0},0),\quad \boldsymbol{k}_{2}=(0,\unicode[STIX]{x03C0}),\end{eqnarray}$$
resulting in the formation of square cells. In the second case, the motion is the superposition of three rolls whose axes intersect at an angle
$(1/3)\unicode[STIX]{x03C0}$
with
$$\begin{eqnarray}\boldsymbol{k}_{1}=(0,\unicode[STIX]{x03C0}),\quad \boldsymbol{k}_{2}=\left(\frac{\sqrt{3}}{2}\unicode[STIX]{x03C0},\frac{1}{2}\unicode[STIX]{x03C0}\right),\quad \boldsymbol{k}_{3}=\left(\frac{\sqrt{3}}{2}\unicode[STIX]{x03C0},-\frac{1}{2}\unicode[STIX]{x03C0}\right),\end{eqnarray}$$
resulting in the formation of hexagonal cells.
4.1 Roll/square interaction
For this case, we take a solution to the
$O(\unicode[STIX]{x1D716}^{0})$
problem in the form
$$\begin{eqnarray}\displaystyle & \displaystyle C_{j}^{0}=\unicode[STIX]{x1D6FA}_{j}\sin (\unicode[STIX]{x03C0}z)\unicode[STIX]{x1D702}^{0}(x,y,\unicode[STIX]{x1D70F}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}^{0}=-\frac{2\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D6E4}}\left[1+\cos (\unicode[STIX]{x03C0}z)\right]\unicode[STIX]{x1D702}^{0}(x,y,\unicode[STIX]{x1D70F}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u^{0}=-2\unicode[STIX]{x03C0}\cos (\unicode[STIX]{x03C0}z)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{0}(x,y,\unicode[STIX]{x1D70F})}{\unicode[STIX]{x2202}x}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle v^{0}=-2\unicode[STIX]{x03C0}\cos (\unicode[STIX]{x03C0}z)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{0}(x,y,\unicode[STIX]{x1D70F})}{\unicode[STIX]{x2202}y}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle w^{0}=-2\unicode[STIX]{x03C0}^{2}\sin (\unicode[STIX]{x03C0}z)\unicode[STIX]{x1D702}^{0}(x,y,\unicode[STIX]{x1D70F}) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\unicode[STIX]{x1D702}^{0}(x,y,\unicode[STIX]{x1D70F})=\mathop{\sum }_{i=1}^{2}A_{i}(\unicode[STIX]{x1D70F})\text{e}^{\text{i}\boldsymbol{k}_{i}\boldsymbol{\cdot }\boldsymbol{ r}}+\text{c.c.},\end{eqnarray}$$
with
$\boldsymbol{k}_{i}$
given in (4.1). The modal amplitudes
$A_{i}(\unicode[STIX]{x1D70F})$
evolve on a slow time
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D716}^{2}t$
as a consequence of the nonlinear interactions in the system, and are to be determined via the solvability conditions. When
$A_{1}=A_{2}\neq 0$
the motion corresponds to three-dimensional squares and when
$A_{i}\neq 0$
and
$A_{j}=0$
for
$i\neq j$
(
$i,j=1,2$
) the motion corresponds to two-dimensional rolls.
The
$O(\unicode[STIX]{x1D716}^{1})$
problem gives rise to a solvability condition
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{1}\mathit{Ra}_{C1}^{1}+\unicode[STIX]{x1D6FA}_{2}\mathit{Ra}_{C2}^{1}=0.\end{eqnarray}$$
The procedure for obtaining this result involves identifying the solution to an adjoint problem, as outlined in the supplementary material (appendix C). Note that
$\mathit{Ra}_{Cj}^{k}=m_{j}\mathit{Ra}^{k}+\mathit{Ra}_{j}^{k}$
(
$j=1,2$
;
$k=0,1,\ldots$
). When
$\mathit{Ra}$
is considered to be the independent bifurcation parameter, we hold
$\mathit{Ra}_{1}$
and
$\mathit{Ra}_{2}$
constant by setting
$\mathit{Ra}_{1}^{1}=\mathit{Ra}_{2}^{1}=0$
, so that (4.5) reduces to
$\mathit{Ra}^{1}=0$
. In the event that
$\mathit{Ra}_{1}$
(or
$\mathit{Ra}_{2}$
) is viewed as the bifurcation parameter, we set
$\mathit{Ra}^{1}=\mathit{Ra}_{2}^{1}=0$
(or
$\mathit{Ra}^{1}=\mathit{Ra}_{1}^{1}=0$
), so that (4.5) reduces to
$\mathit{Ra}_{1}^{1}=0$
(or
$\mathit{Ra}_{2}^{1}=0$
). We conclude that the solvability condition (4.5) is satisfied by
$\mathit{Ra}_{C1}^{1}=\mathit{Ra}_{C2}^{1}=0$
. This is consistent with the time scale
$O(\unicode[STIX]{x1D716}^{2})$
adopted. The solution at this order is detailed in the supplementary material (appendix D).
At
$O(\unicode[STIX]{x1D716}^{2})$
, we apply the solvability conditions to obtain the amplitude equations
$$\begin{eqnarray}\displaystyle & \displaystyle a\frac{\text{d}A_{1}}{\text{d}\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D707}A_{1}-\left(c|A_{1}|^{2}+d|A_{2}|^{2}\right)A_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle a\frac{\text{d}A_{2}}{\text{d}\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D707}A_{2}-\left(c|A_{2}|^{2}+d|A_{1}|^{2}\right)A_{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle a=\left[1+m_{2}\left(\mathit{Le}_{1}\unicode[STIX]{x1D6FA}_{1}-\mathit{Le}_{2}\unicode[STIX]{x1D6FA}_{2}\right)\right]\mathit{Ra}_{C1}^{0}+\left[1+m_{1}\left(\mathit{Le}_{2}\unicode[STIX]{x1D6FA}_{2}-\mathit{Le}_{1}\unicode[STIX]{x1D6FA}_{1}\right)\right]\mathit{Ra}_{C2}^{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}=2\unicode[STIX]{x03C0}^{2}\left(\unicode[STIX]{x1D6FA}_{1}\mathit{Ra}_{C1}^{2}+\unicode[STIX]{x1D6FA}_{2}\mathit{Ra}_{C2}^{2}\right) & \displaystyle\end{eqnarray}$$
$c$
and
$d$
of the cubic terms are given in the supplementary material (appendix E). The form of these equations is identical to those arising in a variety of other physical systems including convection (Jenkins & Proctor Reference Jenkins and Proctor1984; Jenkins Reference Jenkins1987), combustion (Kuske & Matkowsky Reference Kuske and Matkowsky1994) and vegetative morphogenesis (Kealy & Wollkind Reference Kealy and Wollkind2012). Our goal in the next subsection is to extract some physical interpretation of these amplitude equations and their coefficients.
4.1.1 Results
To begin our study of pattern selection described by equations (4.6), we first note that, from (3.7) and (3.6), the oscillation frequency
$\unicode[STIX]{x1D714}^{2}$
can be expressed as
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}=k^{2}a/[(m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1}+1)(m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1})].\end{eqnarray}$$
It follows that the possibility of oscillatory convection is characterized by the condition
$a>0$
. In such a case, the present analysis breaks down owing to the interaction between direct and oscillatory modes. Consequently we focus on the case
$a<0$
in the analysis which follows.
The steady-state solutions of (4.6) are the conduction state (
$A_{1}=A_{2}=0$
), the two-dimensional rolls (one of
$A_{1}$
or
$A_{2}$
zero) and the three-dimensional squares (
$A_{1}=A_{2}\neq 0$
). For an idealized ternary system with a single stratifying solutal agency, i.e.
$\mathit{Ra}=0$
,
$\mathit{Ra}_{1}\neq 0$
and
$\mathit{Ra}_{2}=0$
, recall that both signs of the parameter group
$\unicode[STIX]{x1D6FA}_{1}$
occur, and the results divide naturally into two groups according to
$\unicode[STIX]{x1D6FA}_{1}>0$
(convection under statically unstable conditions) and
$\unicode[STIX]{x1D6FA}_{1}<0$
(convection under statically stable conditions). By analysing the structure and stability of these solutions, we find that the rolls are stable if
$c<0$
and
$d-c<0$
, and the squares are stable if
$c+d<0$
and
$d-c>0$
, irrespective of the sign of
$\unicode[STIX]{x1D6FA}_{1}$
. The results of the analysis are summarized in figure 4(a) for
$\unicode[STIX]{x1D6FA}_{1}>0$
and figure 4(b) for
$\unicode[STIX]{x1D6FA}_{1}<0$
. The solid curves indicate stable solute branches, while the dashed curves unstable ones. We see that stable solutions are present provided both roll and square branches bifurcate supercritically. In other words, it is not possible for a pattern to be stable when another one is subcritical. Further, if both branches bifurcate supercritically, the stable solution is the one with the larger amplitude.
Figure 4. Bifurcation diagrams in the (c,d)-plane, showing the steady-state amplitude of the rolls (R) and squares (S),
$\unicode[STIX]{x1D716}(A_{1}^{2}+A_{2}^{2})^{1/2}$
, as a function of
$\mathit{Ra}_{1}$
in the two cases: (a)
$\unicode[STIX]{x1D6FA}_{1}>0$
and (b)
$\unicode[STIX]{x1D6FA}_{1}<0$
. Solid curves indicate stable finite-amplitude states, while dashed curves correspond to unstable states. The (c,d)-plane divides into six regions, labelled
$\text{I}_{-}$
–
$\text{VI}_{-}$
in (a) and
$\text{I}_{+}$
–
$\text{VI}_{+}$
in (b), each characterized by distinct bifurcation diagrams and stability assignments of roll and square patterns. The subscripts
$\pm$
indicate the sign of the linear critical Rayleigh number. Note that stable solutions can be identified in the shaded regions.
To establish specific predictions for the ternary mush system, we investigate the dependence of the stability assignments of the steady patterns on the control parameters
$m_{j}$
,
$\mathit{Le}_{j}$
,
$\unicode[STIX]{x1D6F7}$
,
$K_{1}$
and
$K_{2}$
. In fact, there are three relevant combinations of the coefficients
$c$
and
$d$
that determine the behaviour:
$c$
,
$d-c$
and
$c+d$
. The two regimes of interest (
$\unicode[STIX]{x1D6FA}_{1}>0$
,
$\unicode[STIX]{x1D6FA}_{1}<0$
) are most clearly distinguished in the two limiting cases,
$\mathit{Le}_{1}\gg 1$
and
$\mathit{Le}_{2}\gg 1$
, which force a clear separation of orders between the neutral stability bound for ‘statically expected’ direct instability (when
$\unicode[STIX]{x1D6FA}_{1}>0$
) and its ‘statically unexpected’ counterpart (when
$\unicode[STIX]{x1D6FA}_{1}<0$
). As an aside, they provide a simple form for the pattern selection criteria. For as
$\mathit{Le}_{1}\rightarrow \infty$
,
$\unicode[STIX]{x1D6FA}_{1}\sim m_{2}\mathit{Le}_{1}/(1-\unicode[STIX]{x1D6F7})$
and (3.5) produces
$\mathit{Ra}_{1}^{0}\sim -4\unicode[STIX]{x03C0}^{2}(1-\unicode[STIX]{x1D6F7})/(m_{2}\mathit{Le}_{1})$
, corresponding to convection onset under statically unstable conditions. As
$\mathit{Le}_{2}\rightarrow \infty$
,
$\unicode[STIX]{x1D6FA}_{1}\sim -m_{2}\mathit{Le}_{2}/(1-\unicode[STIX]{x1D6F7})$
and (3.5) produces
$\mathit{Ra}_{1}^{0}\sim 4\unicode[STIX]{x03C0}^{2}(1-\unicode[STIX]{x1D6F7})/(m_{2}\mathit{Le}_{2})$
, corresponding to onset under statically stable conditions.
A closer inspection of (E.1a,b) in the supplementary material reveals that the coefficients
$c$
and
$d$
are dominated by the terms
$$\begin{eqnarray}c=-\frac{98\unicode[STIX]{x03C0}^{4}}{3}\frac{O{\mathcal{R}}}{\unicode[STIX]{x1D6E4}^{2}(1-\unicode[STIX]{x1D6F7})^{2}}+\frac{\unicode[STIX]{x03C0}^{2}}{6}\frac{O^{2}{\mathcal{R}}^{2}}{\unicode[STIX]{x1D6E4}^{2}(1-\unicode[STIX]{x1D6F7})^{2}}-4\unicode[STIX]{x03C0}^{4}\frac{O{\mathcal{M}}{\mathcal{R}}}{1-\unicode[STIX]{x1D6F7}}+2\unicode[STIX]{x03C0}^{4}\frac{(m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1}){\mathcal{L}}{\mathcal{R}}}{(1-\unicode[STIX]{x1D6F7})^{2}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle d & = & \displaystyle -\frac{1468\unicode[STIX]{x03C0}^{4}}{7}\frac{O{\mathcal{R}}}{\unicode[STIX]{x1D6E4}^{2}(1-\unicode[STIX]{x1D6F7})^{2}}+\frac{445\unicode[STIX]{x03C0}^{2}}{14}\frac{O^{2}{\mathcal{R}}^{2}}{\unicode[STIX]{x1D6E4}^{2}(1-\unicode[STIX]{x1D6F7})^{2}}-5\unicode[STIX]{x03C0}^{4}\frac{O{\mathcal{R}}{\mathcal{M}}}{1-\unicode[STIX]{x1D6F7}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{20\unicode[STIX]{x03C0}^{4}}{7}\frac{(m_{1}\mathit{Le}_{2}+m_{2}\mathit{Le}_{1}){\mathcal{L}}{\mathcal{R}}}{(1-\unicode[STIX]{x1D6F7})^{2}}\end{eqnarray}$$
$\mathit{Le}_{j}$
limits; the omitted terms are higher-order effects. Here
$O=\unicode[STIX]{x1D6FA}_{1}-\unicode[STIX]{x1D6FA}_{2}$
,
${\mathcal{L}}=\mathit{Le}_{1}\unicode[STIX]{x1D6FA}_{1}-\mathit{Le}_{2}\unicode[STIX]{x1D6FA}_{2}$
,
${\mathcal{M}}=m_{1}\mathit{Le}_{1}\unicode[STIX]{x1D6FA}_{1}+m_{2}\mathit{Le}_{2}\unicode[STIX]{x1D6FA}_{2}$
and
${\mathcal{R}}=m_{2}\mathit{Ra}_{C1}^{0}-m_{1}\mathit{Ra}_{C2}^{0}$
. Expressions (4.9) embody a variety of physical effects which play a role in the pattern selection process. The first two terms in each equation represent double-diffusive effects induced by the
$O(\unicode[STIX]{x1D716}^{0})$
and
$O(\unicode[STIX]{x1D716}^{1})$
solid-fraction perturbations. The third term represents a combined effect of double diffusion induced by the
$O(\unicode[STIX]{x1D716}^{0})$
solid-fraction perturbation and double advection, the transport of the two solutes at differing rates, induced by the
$O(\unicode[STIX]{x1D716}^{0})$
flow. The final term arises from double-advective effects associated with the
$O(\unicode[STIX]{x1D716}^{0})$
and
$O(\unicode[STIX]{x1D716}^{1})$
disturbance flows. These particular physical effects control the presence of stable finite-amplitude convective states.
As
$\mathit{Le}_{1}\rightarrow \infty$
, the limiting form of (4.9) is
$$\begin{eqnarray}\displaystyle & \displaystyle c\sim \frac{400\unicode[STIX]{x03C0}^{6}}{3}\frac{m_{1}^{2}}{\left(1-\unicode[STIX]{x1D6F7}\right)^{2}}\mathit{Le}_{1}^{2}+16\unicode[STIX]{x03C0}^{6}\frac{m_{1}m_{2}}{\left(1-\unicode[STIX]{x1D6F7}\right)^{2}}\mathit{Le}_{1}^{2}-8\unicode[STIX]{x03C0}^{6}\frac{m_{2}^{2}}{\left(1-\unicode[STIX]{x1D6F7}\right)^{2}}\mathit{Le}_{1}^{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle d\sim \frac{9432\unicode[STIX]{x03C0}^{6}}{7}\frac{m_{1}^{2}}{\left(1-\unicode[STIX]{x1D6F7}\right)^{2}}\mathit{Le}_{1}^{2}+\frac{148\unicode[STIX]{x03C0}^{6}}{7}\frac{m_{1}m_{2}}{\left(1-\unicode[STIX]{x1D6F7}\right)^{2}}\mathit{Le}_{1}^{2}-\frac{80\unicode[STIX]{x03C0}^{6}}{7}\frac{m_{2}^{2}}{\left(1-\unicode[STIX]{x1D6F7}\right)^{2}}\mathit{Le}_{1}^{2} & \displaystyle\end{eqnarray}$$
$\mathit{Le}_{2}$
. We note also that
$a\sim -4\unicode[STIX]{x03C0}^{2}m_{2}\mathit{Le}_{1}$
, which is negative. The results for the alternate limit
$\mathit{Le}_{2}\rightarrow \infty$
follow from (4.10) on exchanging indices
$1$
and
$2$
.
From these results, we deduce that for instability under statically unstable conditions (
$\mathit{Le}_{1}\gg 1$
), the stable finite-amplitude states occur when
$$\begin{eqnarray}\displaystyle & \displaystyle 0<m_{1}<\frac{3(\sqrt{51073}-21)}{12658}\approx 0.049\quad (\text{rolls of type I}_{-}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{3(\sqrt{51073}-21)}{12658}<m_{1}<\frac{13\sqrt{18993}-399}{14954}\approx 0.093\quad (\text{squares of type II}_{-}). & \displaystyle\end{eqnarray}$$
$\mathit{Le}_{2}\gg 1$
), the conditions for stability are the same as (4.11) except that
$m_{2}$
,
$\text{I}_{+}$
and
$\text{II}_{+}$
rather than
$m_{1}$
,
$\text{I}_{-}$
and
$\text{II}_{-}$
appear. Here the subscripts
$\pm$
indicate the sign of the linear critical Rayleigh number. In either case, these results suggest that, somewhat surprisingly, an essential ingredient in the pattern selection process is the asymmetry between the liquidus slopes
$m_{1}$
and
$m_{2}$
. Specifically, in a system with top-heavy stratification in the slower diffusing agency, the linear onset is reminiscent of the classical buoyancy-driven convection. The nonlinear development of this instability into stable finite-amplitude states of either planform (cf. region
$\text{I}_{-}$
or
$\text{II}_{-}$
of figure 4
a) requires that the slope
$m_{1}$
is sufficiently small, i.e. the scaled rate at which the phase-equilibrium temperature varies with solute-1 concentration is sufficiently small. In a system with bottom-heavy stratification in the faster diffusing agency, the linear onset is characterized by a phase-change-induced double-diffusive convection, despite the statically stable density stratification. Further from onset, this instability develops into stable steady patterns (cf. region
$\text{I}_{+}$
or
$\text{II}_{+}$
of figure 4
b) provided the slope
$m_{1}$
is sufficiently large. Without accounting for the asymmetry in the ternary phase diagram, i.e. when
$m_{1}=m_{2}=1/2$
, no stable convection patterns would be predicted in the limits considered.
Figure 5. Parameter regimes for roll/square interaction. (a)
$m_{1}=0.05$
,
$m_{2}=0.95$
, (b)
$m_{1}=m_{2}=1/2$
, (c)
$m_{1}=0.95$
,
$m_{2}=0.05$
. The remaining parameters are fixed at
$K_{1}=K_{2}=0$
and
$\unicode[STIX]{x1D6F7}=0.1$
. The different regions, labelled by the roman numerals, correspond to the distinct bifurcation diagrams as identified in figure 4. The dark-shaded region corresponds to where an oscillatory instability occurs. For visual purposes, the regions of stable patterns are also light-shaded.
A wider representation of the regions in the
$(\mathit{Le}_{1},\mathit{Le}_{2})$
-plane, where different bifurcation diagrams are located, is shown in figure 5 for (a)
$m_{1}=0.05$
,
$m_{2}=0.95$
, (b)
$m_{1}=m_{2}=1/2$
and (c)
$m_{1}=0.95$
,
$m_{2}=0.05$
. For simplicity, we take
$K_{1}=K_{2}=0$
, which removes the effects of solid-fraction dependence of the permeability. In the dark-shaded region an oscillatory instability precedes the steady-state bifurcation of interest here, and the present analysis breaks down (i.e.
$a\geqslant 0$
). The lower and upper boundaries of this region correspond to
$a=0$
and
$a=\infty$
, respectively. Note that along the latter also
$\unicode[STIX]{x1D6FA}_{1}=0$
, so that the oscillatory region separates a region of static instability below from a region of static stability above. Outside the oscillatory region, the different regions, labelled by
$\text{I}_{-}$
–
$\text{VI}_{-}$
in figure 5(a) and
$\text{I}_{+}$
–
$\text{VI}_{+}$
in figure 5(b), are associated with the distinct bifurcation diagrams located therein (cf. figure 4). The regions where stable finite-amplitude states are predicted,
$\text{I}_{-}$
,
$\text{II}_{-}$
in figure 5(a) and
$\text{I}_{+}$
,
$\text{II}_{+}$
in figure 5(b), are also light-shaded for convenience. For symmetric ternary phase diagrams (figure 5
b), the stable patterns are concentrated in a narrow region around the line
$\mathit{Le}_{1}=\mathit{Le}_{2}$
; this region gets gradually thinner as
$\mathit{Le}_{j}$
increase. Along the line
$\mathit{Le}_{1}=\mathit{Le}_{2}$
the ternary system reduces to an effective binary one, with rolls predicted as a stable pattern. This appears to be a new result in the context of steady convection in binary mushy layers not discussed before.
The nature of the stability characteristics changes dramatically for systems with phase diagrams departed from the symmetric ones. As
$m_{1}$
decreases (figure 5
a), we find that the oscillatory region decreases, and the regions
$\text{I}_{-}$
and
$\text{II}_{-}$
of stable patterns increase at the expense of a much smaller region
$\text{IV}_{-}$
of unstable states. We observe that squares (region
$\text{II}_{-}$
) are the stable pattern at large values of
$\mathit{Le}_{1}$
, while the stable squares are replaced by the stable rolls (region
$\text{I}_{-}$
) as
$m_{1}$
decreases through the value
$m_{1}\approx 0.049$
. The latter feature is illustrated by supplementary figure S1a and is consistent with the asymptotic result (4.11a
) obtained in the limit
$\mathit{Le}_{1}\rightarrow \infty$
.
In contrast, as
$m_{1}$
increases (figure 5
c), the oscillatory region increases, noting that the upper boundary of this region moves towards larger values of
$\mathit{Le}_{2}$
as
$(1-\unicode[STIX]{x1D6F7})/(1-m_{1})$
when
$m_{1}\rightarrow 1$
. Of particular importance is that the parameter plane contains a region
$\text{II}_{+}$
with stable squares. For yet larger values of
$m_{1}$
(see supplementary figure S1b), the stable squares give way to stable rolls of type
$\text{I}_{+}$
, particularly for large
$\mathit{Le}_{2}$
. The changeover occurs at
$m_{1}\approx 0.951$
, consistent with the asymptotic results analogous to (4.11) for the limit
$\mathit{Le}_{2}\rightarrow \infty$
.
4.2 Roll/hexagon interaction
In general, a multiple time-scale approach is required to consistently derive amplitude equations describing the roll/hexagon competition near threshold (see e.g. Fujimura Reference Fujimura1991; Kuske & Milewski Reference Kuske and Milewski1999; Skeldon & Guidoboni Reference Skeldon and Guidoboni2007; Golovin, Matkowsky & Volpert Reference Golovin, Matkowsky and Volpert2008). Here we wish to study a special case where the bifurcation to hexagons is nearly vertical so that a single appropriate time scale is
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D716}^{2}t$
as in the roll/square case; see the supplementary material (appendix B) for further details.
At
$O(\unicode[STIX]{x1D716}^{0})$
, we recover again the linear stability problem with the solution of the same form as (4.3), except for the planform function replaced by
$$\begin{eqnarray}\unicode[STIX]{x1D702}^{0}(x,y,\unicode[STIX]{x1D70F})=\mathop{\sum }_{i=1}^{3}A_{i}(\unicode[STIX]{x1D70F})\text{e}^{\text{i}\boldsymbol{k}_{i}\boldsymbol{\cdot }\boldsymbol{ r}}+\text{c.c.},\end{eqnarray}$$
with
$\boldsymbol{k}_{i}$
given by (4.2). When
$A_{1}=A_{2}=A_{3}\neq 0$
the motion corresponds to three-dimensional hexagons and when
$A_{i}\neq 0$
and
$A_{j}=0$
for
$i\neq j$
the motion corresponds to two-dimensional rolls.
We next consider the solvability condition for
$O(\unicode[STIX]{x1D716}^{2})$
. The particular distinguished scaling taken in the time variable is very important. For had we scaled
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D716}t$
, rather than
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D716}^{2}t$
, the existence of an
$O(\unicode[STIX]{x1D716}^{2})$
solution would require that
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{1}\mathit{Ra}_{C1}^{1}+\unicode[STIX]{x1D6FA}_{2}\mathit{Ra}_{C2}^{1}=12\unicode[STIX]{x03C0}^{3}\frac{K_{1}}{\unicode[STIX]{x1D6E4}}-3\unicode[STIX]{x03C0}\frac{\left(\unicode[STIX]{x1D6FA}_{1}-\unicode[STIX]{x1D6FA}_{2}\right)\left(m_{1}\mathit{Ra}_{C2}^{0}-m_{2}\mathit{Ra}_{C1}^{0}\right)}{\unicode[STIX]{x1D6E4}(1-\unicode[STIX]{x1D6F7})}.\end{eqnarray}$$
Now, in order to ensure that the hexagonal bifurcation is nearly vertical, i.e. the difference
$\unicode[STIX]{x1D6FA}_{1}(\mathit{Ra}_{C1}-\mathit{Ra}_{C1}^{0})+\unicode[STIX]{x1D6FA}_{2}(\mathit{Ra}_{C2}-\mathit{Ra}_{C2}^{0})$
is equal to zero correct to
$O(\unicode[STIX]{x1D716}^{1})$
, we make an assumption that
$$\begin{eqnarray}K_{1}=K_{1c}+\unicode[STIX]{x1D716}\bar{K}_{1},\end{eqnarray}$$
where the critical value
$K_{1c}$
is given by
$$\begin{eqnarray}\frac{K_{1c}}{\unicode[STIX]{x1D6E4}}=\frac{1}{4\unicode[STIX]{x03C0}^{2}}\frac{(\unicode[STIX]{x1D6FA}_{1}-\unicode[STIX]{x1D6FA}_{2})(m_{1}\mathit{Ra}_{C2}^{0}-m_{2}\mathit{Ra}_{C1}^{0})}{\unicode[STIX]{x1D6E4}(1-\unicode[STIX]{x1D6F7})}\end{eqnarray}$$
and
$\bar{K}_{1}$
is an
$O(1)$
detuning parameter that captures the small deviations of the transcritical bifurcation from vertical. A similar approach has been taken by Anderson & Worster (Reference Anderson and Worster1995) within the context of binary mushy layers. It is in the origin of the physical effects which give rise to hexagons that the ternary and binary systems differ. The first term on the right-hand side of (4.13) is the contribution arising from the nonlinear permeability variations due to
$O(\unicode[STIX]{x1D716}^{0})$
solid-fraction perturbations. The second term represents the nonlinear double-diffusive effects induced by
$O(\unicode[STIX]{x1D716}^{0})$
solid-fraction perturbations. It is a combination of these two effects which determines the nature of the transcritical bifurcation to hexagons in the ternary systems.
At
$O(\unicode[STIX]{x1D716}^{2})$
, we apply the solvability conditions to obtain a set of amplitude equations
$$\begin{eqnarray}\displaystyle & \displaystyle a\frac{\text{d}A_{1}}{\text{d}\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D707}A_{1}+bA_{2}A_{3}^{\ast }-cA_{1}|A_{1}|^{2}-d(|A_{2}|^{2}+|A_{3}|^{2})A_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle a\frac{\text{d}A_{2}}{\text{d}\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D707}A_{2}+bA_{1}A_{3}-cA_{2}|A_{2}|^{2}-d(|A_{1}|^{2}+|A_{3}|^{2})A_{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle a\frac{\text{d}A_{3}}{\text{d}\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D707}A_{3}+bA_{1}^{\ast }A_{2}-cA_{3}|A_{3}|^{2}-d(|A_{1}|^{2}+|A_{2}|^{2})A_{3}, & \displaystyle\end{eqnarray}$$
$a$
and
$\unicode[STIX]{x1D707}$
have the same form as (4.7), a coefficient of the quadratic terms 
$$\begin{eqnarray}b=-24\unicode[STIX]{x03C0}^{5}\frac{\bar{K}_{1}}{\unicode[STIX]{x1D6E4}},\end{eqnarray}$$
and the cubic interaction coefficients
$c$
and
$d$
are given explicitly in the supplementary material (appendix E). The form of amplitude equations (4.16) is generic and has been widely studied, particularly in the context of Rayleigh–Bénard convection models (e.g. Busse Reference Busse1978; Bodenschatz, Pesch & Ahlers Reference Bodenschatz, Pesch and Ahlers2000). The bifurcation properties of (4.16) have been discussed from the viewpoint of group theory in considerable detail by Buzano & Golubitsky (Reference Buzano and Golubitsky1983) and Golubitsky, Swift & Knobloch (Reference Golubitsky, Swift and Knobloch1984). As in the case of roll patterns, our goal in the next subsection is to extract some physical interpretation of these amplitude equations and their coefficients.
4.2.1 Results
We now examine the dynamics described by the amplitude equations (4.16). The equations admit steady convecting solutions in the form of two-dimensional rolls (
$A_{1}$
real,
$A_{2}=A_{3}=0$
), three-dimensional regular hexagons (
$A_{1}=A_{2}=A_{3}$
real) and three-dimensional mixed modes (
$A_{1}$
real,
$A_{2}=A_{3}$
real). A linear stability analysis indicates that the roll and hexagonal solutions can be stable, depending on the sign of
$\unicode[STIX]{x1D707}$
,
$b$
,
$c$
and
$d$
. Here we describe only cases which ensure the presence of stable solution branches and which can be invoked via our large-
$\mathit{Le}_{j}$
limits discussed below. It transpires that such situations occur for
$c<0$
,
$c+2d<0$
and
$d-c<0$
and either sign of
$\unicode[STIX]{x1D707}$
and
$b$
. This reduces our exposition to manageable size.
Figure 6. Bifurcation diagrams for the roll/hexagon competition. Sketch of the amplitude
$\unicode[STIX]{x1D716}A_{1}$
versus Rayleigh number
$\mathit{Ra}_{1}$
. Stable solutions are marked by solid curves, while dashed curves indicate unstable solutions. Panels (a,c) show the case of
$\unicode[STIX]{x1D6FA}_{1}>0$
(statically unstable regime), while (b,d) show the case of
$\unicode[STIX]{x1D6FA}_{1}<0$
(statically stable regime). Panels (a,b) show the case of
$\unicode[STIX]{x1D716}b>0$
, while (c,d) show the case of
$\unicode[STIX]{x1D716}b<0$
. Diagrams for
$\unicode[STIX]{x1D716}b<0$
are related to those for
$\unicode[STIX]{x1D716}b>0$
by the sign change
$A_{i}\mapsto -A_{i}$
, as can be seen from (4.16). Note that here
$c<0$
,
$c+2d<0$
and
$d-c<0$
is assumed. The Rayleigh numbers
$\mathit{Ra}_{1}^{0}$
,
$\mathit{Ra}_{1}^{(g)}$
,
$\mathit{Ra}_{1}^{(1)}$
and
$\mathit{Ra}_{1}^{(2)}$
mark correspondingly the linear critical point, a global stability limit point, and the two secondary bifurcation points associated with a stability exchange between rolls and hexagons.
The bifurcation diagrams shown in figure 6 represent the steady-state solution branches in the
$(\mathit{Ra}_{1},\unicode[STIX]{x1D716}A_{1})$
-plane; for simplicity, the
$A_{2}$
and
$A_{3}$
dependence is not exhibited explicitly. Parallel to the analysis of roll/square interaction, we consider the case where a single solute component contributes to buoyancy (
$\mathit{Ra}=0,\mathit{Ra}_{1}\neq 0,\mathit{Ra}_{2}=0$
) and divide the classification into two regimes of interest (
$\unicode[STIX]{x1D6FA}_{1}>0$
in figure 6
a,c;
$\unicode[STIX]{x1D6FA}_{1}<0$
in figure 6
b,d). Further, the sign of
$b$
controls the nature of the hexagonal pattern that is stable: the flow at the cell centres is either upwards when
$b<0$
(figure 6
a,b) or downwards when
$b>0$
(figure 6
c,d). Common to all four cases is the result that the roll solution branch bifurcates supercritically (since
$c<0$
) and the hexagonal solution branch opens in the direction of increasing
$|\mathit{Ra}_{1}|$
(since
$c+2d<0$
). The stable (unstable) sections of the branches are marked solid (dashed). The bifurcation points are indicated by the solid dots. The subcritical portion of the hexagonal branch acquires stability at a turning, or global stability limit, point
$(\mathit{Ra}_{1}^{(g)},\unicode[STIX]{x1D716}A_{1}^{(g)})$
. The hexagons transfer their stability to rolls via an unstable mixed mode solution, with stability exchanges at the secondary bifurcation points
$(\mathit{Ra}_{1}^{(1)},\unicode[STIX]{x1D716}A_{1}^{(1)})$
and
$(\mathit{Ra}_{1}^{(2)},\unicode[STIX]{x1D716}A_{1}^{(2)})$
. The associated Rayleigh numbers and amplitudes are given by
$$\begin{eqnarray}\displaystyle \mathit{Ra}_{1}^{(g)}=-\frac{1}{\unicode[STIX]{x1D6FA}_{1}}\left[4\unicode[STIX]{x03C0}^{2}+\frac{\unicode[STIX]{x1D716}^{2}b^{2}}{8\unicode[STIX]{x03C0}^{2}\left(c+2d\right)}\right],\quad \unicode[STIX]{x1D716}A_{1}^{(g)}=\frac{\unicode[STIX]{x1D716}b}{2\left(c+2d\right)}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \mathit{Ra}_{1}^{(1)}=-\frac{1}{\unicode[STIX]{x1D6FA}_{1}}\left[4\unicode[STIX]{x03C0}^{2}-\frac{\unicode[STIX]{x1D716}^{2}b^{2}c}{2\unicode[STIX]{x03C0}^{2}(d-c)^{2}}\right],\quad \unicode[STIX]{x1D716}A_{1}^{(1)}=\frac{\unicode[STIX]{x1D716}b}{d-c}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \mathit{Ra}_{1}^{(2)}=-\frac{1}{\unicode[STIX]{x1D6FA}_{1}}\left[4\unicode[STIX]{x03C0}^{2}-\frac{\unicode[STIX]{x1D716}^{2}b^{2}\left(2c+d\right)}{2\unicode[STIX]{x03C0}^{2}(d-c)^{2}}\right],\quad \unicode[STIX]{x1D716}A_{1}^{(2)}=\frac{\unicode[STIX]{x1D716}b}{d-c}. & & \displaystyle\end{eqnarray}$$
Referring to figure 6(a), note that there is a hysteresis between the purely diffusive solution and the hexagonal solution in the range
$\mathit{Ra}_{1}^{0}<\mathit{Ra}_{1}<\mathit{Ra}_{1}^{(g)}$
, and between the hexagonal solution and the roll solution in the range
$\mathit{Ra}_{1}^{(2)}<\mathit{Ra}_{1}<\mathit{Ra}_{1}^{(1)}$
.
We now proceed to make specific predictions for the ternary systems in the two large-Lewis-number limits as considered in the roll/square case. First, examining the asymptotic behaviour of (4.17) in light of (4.14), we deduce that
$$\begin{eqnarray}\unicode[STIX]{x1D716}b=-24\unicode[STIX]{x03C0}^{5}\left(K_{1}-\frac{1}{1-\unicode[STIX]{x1D6F7}}\right)m_{j}\mathit{Le}_{j}\end{eqnarray}$$
to leading order in either of the two large-
$\mathit{Le}_{j}$
limits. Observe that the coefficient
$b$
can change sign as a function of the permeability parameter
$K_{1}$
. This implies that the stable hexagons with upflow at the cell centres can be found for
$K_{1}>1/(1-\unicode[STIX]{x1D6F7})$
, while those with downflow at the cell centres can be found for
$K_{1}<1/(1-\unicode[STIX]{x1D6F7})$
. Note that this result is independent of whether the onset is statically expected (figure 6
a,c) or statically unexpected (figure 6
b,d).
The general results (4.18) can be used to give the ranges of
$\mathit{Ra}_{1}$
for stable hexagons and rolls. For the statically unstable case, we can write these as
$$\begin{eqnarray}\displaystyle & \displaystyle -\frac{259}{844}(1-\unicode[STIX]{x1D6F7})^{2}\bar{K}_{1}^{2}m_{1}^{2}\unicode[STIX]{x1D716}^{2}<\frac{\mathit{Ra}_{1}-\mathit{Ra}_{1}^{0}}{\mathit{Ra}_{1}^{0}}<\frac{346\,542}{34\,969}(1-\unicode[STIX]{x1D6F7})^{2}\bar{K}_{1}^{2}m_{1}^{2}\unicode[STIX]{x1D716}^{2}\quad (\text{hexagons}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\mathit{Ra}_{1}-\mathit{Ra}_{1}^{0}}{\mathit{Ra}_{1}^{0}}>\frac{67\,081}{34\,969}(1-\unicode[STIX]{x1D6F7})^{2}\bar{K}_{1}^{2}m_{1}^{2}\unicode[STIX]{x1D716}^{2}\quad (\text{rolls}), & \displaystyle\end{eqnarray}$$
$\mathit{Le}_{1}$
. Note that the range of predicted hexagons approaches zero with both
$m_{1}$
and
$\bar{K}_{1}$
in the limit of
$\mathit{Le}_{1}\rightarrow \infty$
. The corresponding results for the statically stable case can be obtained by taking the limit
$\mathit{Le}_{2}\rightarrow \infty$
,
$m_{2}\rightarrow 0$
to yield relations identical to (4.20) but with
$m_{1}$
replaced with
$m_{2}$
.
Figure 7 shows the computed bifurcation sets in the
$(\mathit{Le}_{j},\mathit{Ra}_{1})$
-plane, indicating the paths of (i) global stability limit points
$\mathit{Ra}_{1}^{(g)}$
, (ii) linear critical points
$\mathit{Ra}_{1}^{0}$
, (iii) bifurcation points
$\mathit{Ra}_{1}^{(1)}$
at which rolls stabilize and (iv) bifurcation points
$\mathit{Ra}_{1}^{(2)}$
at which hexagons destabilize. The value of
$\mathit{Ra}_{1}^{0}$
is rather sensitive to the magnitude of
$\mathit{Le}_{j}$
, so the vertical axis has been scaled to represent the relative distance of
$\mathit{Ra}_{1}$
from the critical Rayleigh number
$\mathit{Ra}_{1}^{0}$
. With this choice, the horizontal axis corresponds to the loci of
$\mathit{Ra}_{1}^{0}$
. The solid curves correspond to the full expressions (4.18), not restricted by the large-
$\mathit{Le}_{j}$
limits, while the dashed lines indicate the asymptotic results as given by (4.20). In the statically unstable regime (figure 7
a), we see that increasing
$\mathit{Le}_{1}$
diminishes the degree of subcriticality of the hexagonal bifurcation, thus rendering the system globally more stable. Also, the region of stable hexagons is narrowed, with rolls preferred over hexagons at sufficiently large
$\mathit{Le}_{1}$
. However, in the statically stable regime (figure 7
b), we see that increasing
$\mathit{Le}_{2}$
destabilizes the system, as indicated by the hexagonal branch becoming more subcritical. At
$\unicode[STIX]{x1D6FA}_{1}=0$
, corresponding to
$\mathit{Le}_{2}=10$
for the case shown, all the transition boundaries coalesce, as confirmed by a local analysis. Also note that the region of stable hexagons grows at the expense of rolls. This suggests that hexagons is a preferred pattern under statically stable conditions in the sense that the range of preferred
$\mathit{Ra}_{1}$
increases with
$\mathit{Le}_{2}$
. It is worth pointing out that the trend in the global stability limit
$\mathit{Ra}_{1}^{(g)}$
is dominated by that in the linear stability limit
$\mathit{Ra}_{1}^{0}$
in both cases figures 7(a) and 7(b). Recall that the two Lewis numbers
$\mathit{Le}_{j}$
influence the linear stability bound through
$\unicode[STIX]{x1D6FA}_{1}$
involving a difference
$\mathit{Le}_{1}-\mathit{Le}_{2}$
. Thus, increasing
$\mathit{Le}_{1}$
has a similar effect to decreasing
$\mathit{Le}_{2}$
. The results shown have
$\bar{K}_{1}=1$
, corresponding to stable up hexagons in both figures 7(a) and 7(b). Note also that the parameter ranges correspond to where
$a<0$
, i.e. away from the region of oscillatory instability. The dashed lines are the asymptotic results as calculated from (4.20), indicating good agreement with the results obtained from the full expressions.
Figure 7. Parameter regimes for roll/hexagon interaction. The paths of (i) global stability limit points, (ii) linear critical points, (iii) secondary bifurcation points along the roll branch and (iv) secondary bifurcation points along the hexagonal branch are shown in (a) the
$(\mathit{Ra}_{1}-\mathit{Ra}_{1}^{0})/\mathit{Ra}_{1}^{0}$
versus
$\mathit{Le}_{1}$
plane for
$\mathit{Le}_{2}=1$
,
$m_{1}=0.01$
and (b) the
$(\mathit{Ra}_{1}-\mathit{Ra}_{1}^{0})/\mathit{Ra}_{1}^{0}$
versus
$\mathit{Le}_{2}$
plane for
$\mathit{Le}_{1}=1$
,
$m_{2}=0.1$
. Here
$(\mathit{Ra}_{1}-\mathit{Ra}_{1}^{0})/\mathit{Ra}_{1}^{0}$
represents the Rayleigh number scaled relative to the linear critical value. The other parameter values are selected at
$\unicode[STIX]{x1D6F7}=0.1$
and
$\bar{K}_{1}=1$
, yielding only (stable) up hexagons in both (a) and (b). The predicted stable convection patterns identified in (a) are statically expected, while those in (b) are statically unexpected (see § 3). The dashed lines are the combined large-
$\mathit{Le}_{j}$
, small-
$m_{j}$
asymptotic results as given by (4.20).
5 Discussion
In a physically relevant situation, we may view the two solute compositions at the mush top,
$C_{jtop}^{\ast }$
, or equivalently the two compositional drops
$\unicode[STIX]{x0394}C_{j}^{\ast }$
(assuming
$C_{jbot}^{\ast }$
fixed), to play the role of the experimental control parameters. Of the eight dimensionless parameters identified in our model problem, though, these control parameters are involved in four (
$m_{j}$
and
$\mathit{Ra}_{j}$
). Thus, while the dimensionless parameters used in the theoretical development are mathematically convenient, they do not isolate the dependence on
$\unicode[STIX]{x0394}C_{j}^{\ast }$
. One way to overcome this difficulty is to define a new set of appropriate scalings that would separate
$\unicode[STIX]{x0394}C_{j}^{\ast }$
as independently controllable parameters. Alternately, one can revert back to dimensional variables. In view of directly addressing only a few material systems below, we follow the second avenue. In particular, we test our linear-theory predictions against the experiments on a ternary alloy
$\text{H}_{2}\text{O}$
–
$\text{KNO}_{3}$
–
$\text{NaNO}_{3}$
by Aitta et al. (Reference Aitta, Huppert and Worster2001a
) performed in the
$\text{H}_{2}\text{O}$
field of the ternary phase diagram, and by Thompson et al. (Reference Thompson, Huppert, Worster and Aitta2003b
) conducted in the
$\text{KNO}_{3}$
field. Perhaps most importantly, we identify a candidate ternary alloy system which allows convection under the statically stable operating conditions, for which there is no experimental confirmation at present.
Table 1 lists the values of the physical constants used for each system. We label the components 1, 2 and 3 in order of decreasing density. Since water as a solvent is the component in excess in each system (component 3), the expansion coefficients
$\unicode[STIX]{x1D6FC}_{j}^{\ast }$
of the two solutes (components 1 and 2) are negative. The transport coefficients
$D_{j}^{\ast }$
have been approximated by the values reported for the dilute aqueous binary systems. Also included are the values of the dimensionless parameters which are independent of the controls
$\unicode[STIX]{x0394}C_{j}^{\ast }$
. Further input parameters, not listed in table 1, include
$\unicode[STIX]{x1D6F1}^{\ast }$
,
$H^{\ast }$
and
$\unicode[STIX]{x1D708}^{\ast }$
. The permeability of the mush is known to be difficult to measure or predict. Here we adopt the relation between the permeability and the local solid fraction from Tait & Jaupart (Reference Tait and Jaupart1992) and evaluate it at the mush top with
$\unicode[STIX]{x1D6F7}=0.1$
, giving the characteristic permeability
$\unicode[STIX]{x1D6F1}^{\ast }=6.235\times 10^{-8}~\text{m}^{2}$
. In the typical experiments on the solidification from a fixed cooled base, such as those by Aitta et al. (Reference Aitta, Huppert and Worster2001a
) and Thompson et al. (Reference Thompson, Huppert, Worster and Aitta2003b
), the mush thickness
$H^{\ast }$
is time-dependent and not known in advance. Our simple model is not designed to address this issue, and in our calculations we use a constant value
$H^{\ast }=10^{-2}$
m, corresponding roughly to the reported experimental values at the commencement of convection. For the kinematic viscosity of the solution, we take a single value
$\unicode[STIX]{x1D708}^{\ast }=\unicode[STIX]{x1D707}^{\ast }/\unicode[STIX]{x1D70C}_{top}^{\ast }=9.3\times 10^{-7}~\text{m}^{2}~\text{s}^{-1}$
(Worster & Kerr Reference Worster and Kerr1994) for all cases.
Table 1. Parameter values. Three sets of physical properties appropriate to the
$\text{H}_{2}\text{O}$
–
$\text{KNO}_{3}$
–
$\text{NaNO}_{3}$
system solidified in the
$\text{H}_{2}\text{O}$
field (first column), the
$\text{H}_{2}\text{O}$
–
$\text{KNO}_{3}$
–
$\text{NaNO}_{3}$
system solidified in the
$\text{KNO}_{3}$
field (second column) and the
$\text{H}_{2}\text{O}$
–
$\text{NH}_{4}\text{Cl}$
–
$\text{KNO}_{3}$
solidified in the
$\text{H}_{2}\text{O}$
field (third column) of the ternary phase diagram. The values of
$m_{j}^{\ast }$
have been calculated from the available phase diagram properties. The values of
$D_{j}^{\ast }$
and
$\unicode[STIX]{x1D6FC}_{j}^{\ast }$
account for slight but essential differences in material properties of the different solute components involved. For simplicity, all the three systems are characterized by a single set of parameter values for
$\unicode[STIX]{x1D705}^{\ast }$
and
$\unicode[STIX]{x1D6FC}^{\ast }$
. The quantities distinguished by superscript
$^{\dagger }$
correspond to values estimated from the dilute binary systems (no ternary component). The principal dimensionless parameters independent of the compositional controls
$\unicode[STIX]{x0394}C_{j}^{\ast }$
are listed in the final seven rows.
a Aitta et al. (Reference Aitta, Huppert and Worster2001a ).
b Thompson et al. (Reference Thompson, Huppert, Worster and Aitta2003b ).
c Tait & Jaupart (Reference Tait and Jaupart1992).
d Yeh & Wills (Reference Yeh and Wills1970).
e Harned & Hudson (Reference Harned and Hudson1951).
f Worster & Kerr (Reference Worster and Kerr1994).
g Bennon & Incropera (Reference Bennon and Incropera1987).
The condition representing net static stability,
$\mathit{Ra}+\mathit{Ra}_{1}+\mathit{Ra}_{2}=\mathit{Ra}_{C1}+\mathit{Ra}_{C2}>0$
, can be expressed as
$$\begin{eqnarray}(1+R_{1})\unicode[STIX]{x1D6FC}_{1}^{\ast }\unicode[STIX]{x0394}C_{1}^{\ast }+(1+R_{2})\unicode[STIX]{x1D6FC}_{2}^{\ast }\unicode[STIX]{x0394}C_{2}^{\ast }>0,\end{eqnarray}$$
where
$R_{j}\equiv \unicode[STIX]{x1D6FC}^{\ast }m_{j}^{\ast }/\unicode[STIX]{x1D6FC}_{j}^{\ast }$
are the thermal-to-solutal buoyancy ratios. Since, typically,
$|R_{j}|<1$
(see table 1), convection in the primary mush is dictated by the two solute fields and, on noting that
$\mathit{Ra}_{Cj}\equiv (1+R_{j})\mathit{Ra}_{j}$
, is quantified by the two solutal Rayleigh numbers
$\mathit{Ra}_{j}$
.
Individual static stability requires
$\mathit{Ra}_{1}>0$
and
$\mathit{Ra}_{2}>0$
, i.e.
$$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{1}^{\ast }\unicode[STIX]{x0394}C_{1}^{\ast }>0\quad \text{and}\quad \unicode[STIX]{x1D6FC}_{2}^{\ast }\unicode[STIX]{x0394}C_{2}^{\ast }>0.\end{eqnarray}$$
For all the systems considered below, we associate the solutes 1 and 2 with the heavier components of a ternary alloy so that both
$\unicode[STIX]{x1D6FC}_{1}^{\ast }$
and
$\unicode[STIX]{x1D6FC}_{2}^{\ast }$
are negative, implying that the region of individual static stability corresponds to the third quadrant of the
$(\unicode[STIX]{x0394}C_{1}^{\ast },\unicode[STIX]{x0394}C_{2}^{\ast })$
-plane in all cases.
Finally, a condition for dynamic (linear) instability of the primary mush,
$\unicode[STIX]{x1D6FA}\mathit{Ra}_{C1}+\unicode[STIX]{x1D6FA}\mathit{Ra}_{C2}<-4\unicode[STIX]{x03C0}^{2}$
, can be most conveniently expressed in dimensional terms as
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{g^{\ast }\unicode[STIX]{x1D6F1}^{\ast }H^{\ast }}{\unicode[STIX]{x1D708}^{\ast }\unicode[STIX]{x1D705}^{\ast }}\unicode[STIX]{x1D6FC}_{1}^{\ast }\unicode[STIX]{x0394}C_{1}^{\ast }\left\{(1+R_{1})+\left(1+R_{2}\right)\frac{\unicode[STIX]{x1D6FC}_{2}^{\ast }}{\unicode[STIX]{x1D6FC}_{1}^{\ast }}\frac{\unicode[STIX]{x0394}C_{2}^{\ast }}{\unicode[STIX]{x0394}C_{1}^{\ast }}+\frac{1}{1-\unicode[STIX]{x1D6F7}}\left[(1+R_{1})\frac{m_{2}^{\ast }}{m_{1}^{\ast }}-(1+R_{2})\frac{\unicode[STIX]{x1D6FC}_{2}^{\ast }}{\unicode[STIX]{x1D6FC}_{1}^{\ast }}\right]\right.\nonumber\\ \displaystyle & & \displaystyle \qquad \quad \times \,\frac{\unicode[STIX]{x1D705}^{\ast }(D_{2}^{\ast }-D_{1}^{\ast })}{D_{1}^{\ast }D_{2}^{\ast }}\left(1+\frac{m_{2}^{\ast }}{m_{1}^{\ast }}\frac{\unicode[STIX]{x0394}C_{2}^{\ast }}{\unicode[STIX]{x0394}C_{1}^{\ast }}\right)^{-1}\left.\frac{\unicode[STIX]{x0394}C_{2}^{\ast }}{\unicode[STIX]{x0394}C_{1}^{\ast }}\right\}<-4\unicode[STIX]{x03C0}^{2}.\end{eqnarray}$$
Given that temperature and solute concentrations conform to the liquidus relationship, we regard the temperature drop
$\unicode[STIX]{x0394}T^{\ast }$
to be compositionally controlled through the two compositional drops
$\unicode[STIX]{x0394}C_{j}^{\ast }$
via
$\unicode[STIX]{x0394}T^{\ast }=m_{1}^{\ast }\unicode[STIX]{x0394}C_{1}^{\ast }+m_{2}^{\ast }\unicode[STIX]{x0394}C_{2}^{\ast }$
, a relationship which has been explicitly used in writing (5.3). Below, we examine only a portion of the
$(\unicode[STIX]{x0394}C_{1}^{\ast },\unicode[STIX]{x0394}C_{2}^{\ast })$
-plane which corresponds to the case of cooling from below, so that
$\unicode[STIX]{x0394}T^{\ast }$
is positive. While accessible to inspection, the case of heating from below (
$\unicode[STIX]{x0394}T^{\ast }<0$
) is not pursued since the usual thermal buoyancy for the liquids with
$\unicode[STIX]{x1D6FC}^{\ast }>0$
would then mask an examination of compositional convection under statically stable conditions.
The results for the
$\text{H}_{2}\text{O}$
–
$\text{KNO}_{3}$
–
$\text{NaNO}_{3}$
system, with
$\text{H}_{2}\text{O}$
solidified in the primary mush, are represented in figure 8. As in figure 2, the coloured regions indicate the parameter space in which the system is net statically stable (yellow), individually statically stable (cyan) and dynamically unstable (grey shaded). We find that the onset of direct modes of instability occurs inside the region of net static stability. The present system does not display the oscillatory behaviour; the values of
$\unicode[STIX]{x0394}C_{j}^{\ast }$
needed to achieve an oscillatory branch of the neutral stability curve correspond to unphysically large concentrations. The data points correspond to the experimental operating conditions investigated by Aitta et al. (Reference Aitta, Huppert and Worster2001a
) exhibiting diffusion-controlled solidification with no convection. The experimental data are consistent with our predictions from a simplified stability problem for the primary mush.
Figure 8. Regime diagram for the
$\text{H}_{2}\text{O}$
–
$\text{KNO}_{3}$
–
$\text{NaNO}_{3}$
system with primary solidification in the
$\text{H}_{2}\text{O}$
-field of the ternary phase diagram. The diagram quantifies the operating conditions in the
$(\unicode[STIX]{x0394}C_{1}^{\ast },\unicode[STIX]{x0394}C_{2}^{\ast })$
-plane at which different regimes of static and dynamic stability are predicted. Here
$C_{1}^{\ast }$
and
$C_{2}^{\ast }$
are interpreted as the compositions of
$\text{NaNO}_{3}$
and
$\text{KNO}_{3}$
, respectively. Yellow-coloured is the region of net static stability, while cyan-coloured is the region of individual static stability. The solid curve indicates the direct branch of the neutral stability curve. The system is linearly unstable in the grey-shaded area. A description is not pursued in the hatched area which corresponds to a negative temperature gradient (
$\unicode[STIX]{x0394}T^{\ast }<0$
, heating from below). The data points correspond to the experimental regimes investigated by Aitta et al. (Reference Aitta, Huppert and Worster2001a
) (consistent with their experiments, no convective instabilities are expected in this region). The inset shows in the
$(C_{1}^{\ast },C_{2}^{\ast })$
-plane the solidification paths that descend along the liquidus surface until they intersect the cotectic boundaries (dashed), corresponding to the formation of primary mushy layers. Note that for small solute diffusivities the solidification paths nearly coincide with the tie-lines connecting the origin (pure
$\text{H}_{2}\text{O}$
) to the initial compositions (data points).
Turning to a convective situation, Thompson et al. (Reference Thompson, Huppert, Worster and Aitta2003b
) have considered experiments on the same physical system as Aitta et al. (Reference Aitta, Huppert and Worster2001a
), but with the initial concentrations of the dilute alloy falling in the
$\text{KNO}_{3}$
, rather than the
$\text{H}_{2}\text{O}$
, field of the ternary phase diagram. The linear-theory predictions for this configuration are represented in figure 9. In the statically unstable regime (second quadrant), instability is predicted and observed because of a density reversal due to an unstable compositional stratification of
$\text{KNO}_{3}$
. The comparison with the experimental data, taken from Thompson et al. (Reference Thompson, Huppert, Worster and Aitta2003b
), is less satisfactory, with the data lying a little on the left of our computed theoretical neutral stability bound. A little improvement is provided by a region of subcritical behaviour identified to the left of the upper stability branch, but this region is too narrow to be discernible on the scale of the figure. The theoretical predictions would also suggest that instability might be expected in the present system even with a fluid stably stratified in both individual solute components (third quadrant), although if
$\text{KNO}_{3}$
component becomes stably stratified the corresponding liquid line of descent would markedly deviate from the tie-lines passing through the
$\text{KNO}_{3}$
corner, which appears unlikely.
Figure 9. As in figure 8, but for the
$\text{H}_{2}\text{O}$
–
$\text{KNO}_{3}$
–
$\text{NaNO}_{3}$
system with primary solidification in the
$\text{KNO}_{3}$
-field of the ternary phase diagram. Here
$C_{1}^{\ast }$
and
$C_{2}^{\ast }$
are interpreted as the compositions of
$\text{NaNO}_{3}$
and
$\text{KNO}_{3}$
, respectively. The data points correspond to the experimental regimes investigated by Thompson et al. (Reference Thompson, Huppert, Worster and Aitta2003b
).
Figure 10. As in figure 8, but for the
$\text{H}_{2}\text{O}$
–
$\text{NH}_{4}\text{Cl}$
–
$\text{KNO}_{3}$
system with primary solidification in the
$\text{H}_{2}\text{O}$
-field of the ternary phase diagram. Here
$C_{1}^{\ast }$
and
$C_{2}^{\ast }$
are interpreted as the compositions of
$\text{KNO}_{3}$
and
$\text{NH}_{4}\text{Cl}$
, respectively. The red solid curves in the main plot represent the neutral critical boundary for three different values of the mush thickness
$H^{\ast }$
(indicated in
$\text{m}^{2}$
by the numbers on the graph). Note the prediction of convective instability (grey-shaded area) in the third quadrant of the
$(\unicode[STIX]{x0394}C_{1}^{\ast },\unicode[STIX]{x0394}C_{2}^{\ast })$
-plane; i.e. inside the region of individual static stability where convective instability would not typically be anticipated (our ‘unexpected’ instability region). The cross marks a sample data point which corresponds to a proposed solidification path shown in the inset. The compositions at the ternary eutectic point were taken from a database of the calculated phase diagrams (Bale et al.
Reference Bale, Chartrand, Degterov, Eriksson, Hack, Ben Mahfoud, Melanon, Pelton and Petersen2002).
A necessary condition for the prediction of convection in the statically stable regime is that the third term in the curly brackets of (5.3), arising from the phase-change-induced double-diffusive solutal effects, is negative. Assuming that
$\unicode[STIX]{x1D6FC}_{j}^{\ast }<0$
,
$m_{j}^{\ast }<0$
and processing conditions
$\unicode[STIX]{x0394}C_{j}^{\ast }<0$
, this criterion amounts to either
$\unicode[STIX]{x1D6FC}_{2}^{\ast }/\unicode[STIX]{x1D6FC}_{1}^{\ast }<m_{2}^{\ast }/m_{1}^{\ast }$
,
$D_{2}^{\ast }/D_{1}^{\ast }<1$
or
$\unicode[STIX]{x1D6FC}_{2}^{\ast }/\unicode[STIX]{x1D6FC}_{1}^{\ast }>m_{2}^{\ast }/m_{1}^{\ast }$
,
$D_{2}^{\ast }/D_{1}^{\ast }>1$
. Thus, whether such an instability is experimentally realizable depends on the availability of ternary alloy systems with the proper ratios
$\unicode[STIX]{x1D6FC}_{2}^{\ast }/\unicode[STIX]{x1D6FC}_{1}^{\ast }$
,
$m_{2}^{\ast }/m_{1}^{\ast }$
and
$D_{2}^{\ast }/D_{1}^{\ast }$
. An example of an alloy which satisfies these requirements is provided by the
$\text{H}_{2}\text{O}$
–
$\text{NH}_{4}\text{Cl}$
–
$\text{KNO}_{3}$
system, with the linear-theory predictions shown in figure 10. The neutral critical stability boundaries (red solid curves in the main plot) are shown for three different values of the mush thickness
$H^{\ast }$
to illustrate the influence of changing
$H^{\ast }$
in (5.3). We observe a stabilization (larger values of
$|\unicode[STIX]{x0394}C_{j}^{\ast }|$
required for instability) and an ultimate suppression of the statically unexpected convective instability as indicated by shifting the stability boundary towards the left bottom corner of the diagram, suggesting that this instability is promoted by thin mushy layers. Note that the similar change in
$H^{\ast }$
has indistinguishable effects on the results displayed in figures 8 and 9. We must note that, unlike the typical experiments, our simple model omits the effects of solute rejection and background solidification. In addition, the model does not account for the effects of secondary solidification and treats the primary mushy layer in dynamic isolation from the rest of the system. On the other hand, as noted earlier, the linear stability results for the system under consideration show that the presence of these novel convective modes are fairly robust with respect to variations in the parameter values. Further, based on the experience from binary alloy studies, the inclusion of additional physical effects is expected to change the quantitative results, but not the qualitative predictions in terms of the existence and the nonlinear properties of a statically unexpected convection revealed here. We therefore expect that the present work may provide a framework within which to motivate and interpret new experiments on ternary alloy solidification in mushy layers.
6 Concluding remarks
We have examined the weakly nonlinear stability of convection in a mushy layer during primary solidification of ternary alloys. Based on a single-layer model introduced by Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010), our analysis extends to examine the nonlinear development of steady convection in the form of two-dimensional roll, and three-dimensional square and hexagonal patterns. A detailed calculation to determine the various essential coefficients in the nonlinear amplitude equations governing roll/square and roll/hexagon competitions have been performed. Consideration of their limiting forms in special cases has laid the framework for the physical interpretation of nonlinear interactions present in the system.
Guided by the linear stability results of Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010), in an idealized case with a single stratifying agency (solute 1), we have distinguished between two cases: the intuitive, statically expected onset of steady convection appropriate when
$\mathit{Le}_{1}\rightarrow \infty$
; and the counter-intuitive, statically unexpected onset when
$\mathit{Le}_{2}\rightarrow \infty$
. We have obtained a fairly complete classification of the possible bifurcations of steady-state solutions off the purely diffusive, thermostatic equilibrium state (no background solidification, no solute rejection, no latent-heat release) in both limiting cases. In the large-
$\mathit{Le}_{2}$
limit, in particular, the linear onset is characterized by a phase-change-induced double-diffusive mechanism envisaged by Anderson et al. (Reference Anderson, McFadden, Coriell and Murray2010). Away from onset, we find that this instability develops into stable roll or square pattern provided the liquidus surface is sufficiently flat with respect to a slower diffusing solutal agency (solute 2); see § 4.1. Further, hexagons rather than rolls are favoured as the stable convection pattern if, in addition, there are sufficiently strong nonlinear permeability variations in the mush; see § 4.2.
We have considered a situation that departed an
$O(\unicode[STIX]{x1D716})$
from the vertical bifurcation to hexagonal patterns. This assumption was made to facilitate the computation without resorting to the use of a more involved asymptotic analysis based on the method of multiple (time) scales. Examination of the weakly nonlinear stability theory shows that this approximation in fact allows a prediction of stable finite-amplitude solutions in the form of rolls and hexagons, particularly in the regime of large Lewis numbers, as discussed in § 4.2. Thus, the near-vertical-bifurcation limit taken is of dual interest. As a check on this calculation, the interaction coefficients in the amplitude equations were determined for the general case without any restriction on the nature of hexagonal bifurcation. By confining to the near-vertical-bifurcation case, this produced results identical with those recorded above.
A feature of the simplified model studied here is the absence of a number of phase-change effects such as background solidification (
$V=0$
), solute rejection in the local solute balances (
$k_{j}=1$
) and latent-heat release (
$S=0$
). The robustness of the novel instabilities identified in previous linear stability analyses and, in particular, their presence even in highly simplified scenarios such as the case considered here, make it reasonable to expect that the present nonlinear results carry a similar robustness. That said, the convective instability develops on a thermal diffusion time scale
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D705}^{\ast }}^{\ast }=H^{\ast 2}/\unicode[STIX]{x1D705}^{\ast }$
based on the mush thickness
$H^{\ast }$
. In their study of a distinguished limit of the primary mush equations which occurs for
$V\rightarrow 0$
, Guba & Anderson (Reference Guba and Anderson2014) have revealed that the presence of background solidification (
$V\neq 0$
) and solute rejection (
$k_{j}\neq 1$
) in conjunction similarly gives rise to a diffusive destabilization of the system but the underlying mechanism is different. The instability develops on a much longer thermal diffusion time scale
$\unicode[STIX]{x1D70F}_{V^{\ast }}^{\ast }=\unicode[STIX]{x1D705}^{\ast }/V^{\ast 2}$
dependent on the background speed
$V^{\ast }$
, noting that
$\unicode[STIX]{x1D70F}_{V^{\ast }}^{\ast }\gg \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D705}^{\ast }}^{\ast }$
. It would be desirable to extend the analysis of Guba & Anderson (Reference Guba and Anderson2014) not only to nonlinear regime but also to establish a formal link with the nonlinear results identified here. From such extensions, further insight will be gained into the multicomponent phase-change problems in which mass diffusion has a controlling, dynamical influence on the bulk fluid flow.
Acknowledgements
This work was supported by the Slovak Research and Development Agency under the contract No. APVV-14-0378, the VEGA grants 1/0319/15 and 2/0115/16 (P.G.), and by the U.S. National Science Foundation through DMS-1107848 (D.M.A.). We would like to thank M. Revallo for helpful discussions.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2017.382.
