1. Introduction
The properties of cast alloys are affected to a large extent by the local solidification conditions, which are in turn strongly coupled to the thermal and solutal transport during solidification and the associated buoyancy-driven flows in the bulk liquid. Therefore, it is desirable to understand the details of the transport phenomena that take place during solidification in order to predict the eventual properties. Owing to the complex nature of the governing equations, numerical methods, such as the mixture model detailed by Bennon & Incropera (Reference Bennon and Incropera1987a ,Reference Bennon and Incropera b ), have generally been used to predict process behaviour. However, these simulations are computationally expensive and specific to a particular material system and set of process parameters, and general functional relationships are not easily identified.
Several researchers have used analytical or semi-analytical methods to investigate temperature fields during alloy solidification under certain simplifying assumptions in order to determine relevant metallurgical quantities such as the solidification front position and velocity and the local solidification time (LST). These analyses include similarity solutions by Voller (Reference Voller1997) and Chung et al. (Reference Chung, Lee, Choi and Yoo2001) for multicomponent alloys, applications of Goodman’s (Reference Goodman1958) heat balance integral method by Tien & Geiger (Reference Tien and Geiger1967, Reference Tien and Geiger1968), Sunderland and coworkers (Cho & Sunderland Reference Cho and Sunderland1969; Muehlbauer et al. Reference Muehlbauer, Hatcher, Lyons and Sunderland1973) and Voller (Reference Voller1989), and scaling analyses by Krane & Incropera (Reference Krane and Incropera1996) and Amberg (Reference Amberg1997). These analyses make various assumptions about the boundary conditions and the initial temperature of the liquid, but none adequately model the effect of the flow of the bulk liquid during the process.
The only approximate analysis of metallic alloy solidification that considered the bulk fluid motion was Amberg’s (Reference Amberg1997) scaling analysis. Amberg developed a parameter map for determining the main solidification regimes as a function of the system parameters and material properties, although functional relationships were not found for all regimes. He assumed that a steady-state laminar natural convection boundary layer formed at the liquidus interface. A scaled form of the Nusselt number was used to predict the heat flux at the interface and therefore the time required for the removal of any superheat in the bulk liquid. Amberg’s analysis, however, does not consider the effect of the changing bulk liquid temperature on the boundary layer flow, the changing amount of liquid due to solidification or the coupling of the time-dependent liquid temperature on the solidification conditions.
The fluid flow involved in the solidification of alloys is driven by density differences caused by temperature and composition gradients. Significant effort has been put forth to understand convection under such combined driving forces, often termed double-diffusive convection, and the reader is referred to extensive reviews by Turner (Reference Turner1974) and Huppert & Turner (Reference Huppert and Turner1981). Specific studies on solidification phenomena that include these effects have been performed by Thompson & Szekely (Reference Thompson and Szekely1988) and Jarvis & Huppert (Reference Jarvis and Huppert1995). However, these analyses are applied to aqueous salt solidification or geological phenomena with large Prandtl numbers and assume an infinitely thin mushy zone, whereas small Prandtl numbers (
$\mathit{Pr}\ll 1$
) and large freezing ranges are characteristic of metallic alloys.
The purpose of the current study is to analyse the transient solidification behaviour, including the coupled effect of the bulk fluid flow, for metallic alloys with
$\mathit{Pr}\ll 1$
. Under the limiting assumption of thermally dominated convection, the flow configured is expected to be similar to that studied by Lin, Armfield & Patterson (Reference Lin, Armfield and Patterson2007). They considered the transient motion of a liquid in a rectangular cavity suddenly cooled at an isothermal sidewall and predicted the time for the fluid to cool to the wall temperature. They approximated the fluid flow with a laminar boundary layer solution at the cooling wall and tracked the volume of liquid that was cooled by passing through that boundary layer.
The liquid region in an alloy cooled from one side of a rectangular cavity has a geometry similar to that investigated by Lin et al. (Reference Lin, Armfield and Patterson2007), except that the chilled wall is at the liquidus temperature and the domain shrinks as the liquid region is consumed by the solidification front. In the present work, the analysis of Lin et al. is adapted for the fluid flow down the solidification front as the bulk liquid is cooled and allowed to solidify over the alloy’s freezing range. The system configuration is detailed in § 2. This flow solution is coupled to the solidification behaviour, which is analysed using both scaling and integral methods in §§ 3 and 4, respectively. The scaling analysis is based on the one-dimensional solution proposed by Krane & Incropera (Reference Krane and Incropera1996), but is extended to consider the effects of the liquid flow, the latent heat release in the mushy zone and the transient temperature of the bulk liquid. The goal of this analysis is to determine the functional relationships among the material properties, process conditions and the system behaviour. The results are compared to predictions of a numerical mixture model for a binary alloy in § 5. A detailed comparison of the flow field analysis to the work of Lin et al. (Reference Lin, Armfield and Patterson2007) is made in § 6. The scaling and integral analyses are derived based upon assumptions of thermally dominated flow and an impermeable mushy zone, which are evaluated and the applicability of the analyses discussed in § 7. Concluding remarks are presented in § 8.
2. Problem statement
The system under consideration is a binary alloy with nominal composition,
$C_{0}$
, solidifying in the rectangular domain shown in figure 1. The thermal boundary conditions are three insulated walls and one vertical isothermal wall suddenly changed at time
$t=0$
from the initial temperature,
$T_{0}$
, to the wall temperature,
$T_{w}$
, below the solidus temperature of the alloy,
$T_{sol}$
. The initial condition is a uniform temperature,
$T_{0}$
, in which the liquid is superheated above the alloy’s liquidus temperature,
$T_{liq}$
. The temperature averaged over the liquid region,
$T_{avg}(t)$
, decreases with time from
$T_{0}$
. The solidification front is treated as one-dimensional; its shape unaffected by the fluid flow. Advection in the mushy zone is neglected and all properties are uniform and constant.
Figure 1. (a) A simplified binary alloy phase diagram and (b) the system schematic.
Under these conditions, solid immediately forms at the chill, resulting in three distinct spatial regions for
$t>0$
: (i) solid metal in contact with the wall, (ii) a two-phase mushy region, and (iii) bulk liquid. Temporally, there are two regimes, which must be considered individually. The first includes all three spatial regions and begins at the start of the process, as shown in figure 1(b); the second occurs after the liquidus front reaches the far wall at
$x=L$
, leaving only solid and mush. The following analyses consider these two distinct regimes and determine the time at which the system transitions from the first to the second.
As solidification commences, the negative buoyancy from local cooling of the bulk fluid causes a natural convection boundary layer to form at the solidification front. This flow is assumed to be laminar. The fluid cooled by the edge of the mushy zone at the liquidus temperature flows down that interface and fills the cavity at the expense of fluid at the initial temperature, acting to extinguish the superheat of the fluid. This flow configuration is seen in temperature fields shown in figure 2, obtained from a preliminary numerical simulation (described in § 5) for an Al–4.5 wt% Cu alloy, in which cooling of the bulk fluid by the boundary layer along the solidification front is clearly in evidence. The near-vertical solidification shown in that figure and in figure 3 by the liquidus isotherm also supports a one-dimensional approximation of the solidification behaviour.
Figure 2. Numerical simulation of an Al–4.5 wt% Cu alloy in a 0.25 m square cavity initially at 1000 K and cooled from the left wall at 500 K, at times (a)
$t=10$
s, (b)
$t=30$
s, (c)
$t=50$
s, (d)
$t=70$
s, (e)
$t=90$
s and (f)
$t=110$
s. The contours are isotherms 5 K apart starting at the liquidus temperature and increasing to the right to the initial temperature.
Figure 3. Streamlines showing bulk fluid flow during solidification of an Al–4.5 wt% Cu alloy initially at 1000 K and cooled by the left isothermal wall at 500 K, at times (a)
$t=10$
s, (b)
$t=30$
s, (c)
$t=50$
s, (d)
$t=70$
s, (e)
$t=90$
s and (f)
$t=110$
s.
Figure 3 shows streamlines of the flow in the bulk fluid for the simulation in figure 2. A natural convection boundary layer clearly forms on the solidification front, moving down the entire height (
$H$
) of the cavity, finally flowing out into the bulk fluid at the bottom. This boundary layer is fed by the entrainment of the bulk liquid. The circulating fluid is cooled at the liquidus interface, decreasing the average temperature of the bulk fluid and reducing the driving force for the natural convection boundary layer.
3. Scaling analysis
3.1. Scaling the solidification behaviour
The simplified one-dimensional governing equation for energy conservation in the solid and mushy regions, including the effect of the latent heat release but neglecting advection in the mushy zone, is
$$\begin{eqnarray}\frac{\partial T}{\partial t}={\it\alpha}\frac{\partial ^{2}T}{\partial x^{2}}+\frac{L_{f}}{c_{p}}\frac{\text{d}f_{s}}{\text{d}t},\end{eqnarray}$$
where
$T$
is temperature,
$t$
is time,
${\it\alpha}$
is thermal diffusivity,
$L_{f}$
is latent heat of fusion,
$c_{p}$
is constant-pressure specific heat and
$f_{s}$
is the mass fraction of solid. The effect of the fluid flow and the heat transfer in the bulk fluid seen in figures 2 and 3 between the mushy zone and the bulk liquid may be approximated with a quasi-static natural convection boundary layer solution for an isothermal flat plate at the liquidus temperature as described in § 3.2.2. The time for the initial transient leading-edge effect to reach the bottom of the cavity is assumed small compared to the total solidification time. (This assumption is evaluated during the course of the fluid flow scaling analysis below.) The scaling solution for solidification begins for the first temporal regime, in which all three regions are present, by applying the simplified energy equation to the mushy zone and the thermal boundary layer. The following reference differences are used:
$$\begin{eqnarray}\displaystyle & {{\rm\Delta}T_{ref}\sim {\rm\Delta}T}_{sol}=T_{avg}-T_{sol}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\rm\Delta}x_{ref}\sim {\it\delta}_{m1}+\bar{{\it\delta}}_{T}. & \displaystyle\end{eqnarray}$$
${\rm\Delta}T_{ref}$
is a reference temperature difference equal to
${\rm\Delta}T_{sol}$
, which is the difference between the average bulk fluid temperature
$T_{avg}$
and the solidus temperature of the alloy
$T_{sol}$
. The characteristic reference distance
${\rm\Delta}x_{ref}$
is the width of the mushy zone during the first solidification regime
${\it\delta}_{m1}$
plus the mean thermal boundary layer thickness
$\bar{{\it\delta}}_{T}$
. For the subsequent analysis, it is temporarily assumed that the average liquid temperature is constant with time. The scaled energy equation is then 
$$\begin{eqnarray}\frac{{\rm\Delta}T_{sol}}{t}\sim \frac{{\it\alpha}{\rm\Delta}T_{sol}}{({\it\delta}_{m1}+\bar{{\it\delta}}_{T})^{2}},\frac{L_{f}}{c_{p}}\frac{{\rm\Delta}f_{s}}{t},\end{eqnarray}$$
where
${\rm\Delta}f_{s}$
is the reference difference for the solid mass fraction. The first term represents the sensible heat storage, the second the conduction over the boundary layer and mushy zone, and the third the latent heat storage. An obvious disadvantage of this scaling analysis is the use of a linear temperature profile over the mushy zone, meaning that the heat fluxes at the solid–mush and mush–boundary layer interfaces cannot both be matched unless they are forced to be equal. This shortcoming is circumvented here by assuming that the boundary layer is very thin such that the temperature at the liquidus interface is equal to the bulk fluid temperature. Accordingly, the boundary layer thickness is neglected in (3.3). The change in the fraction solid over the mushy zone is
$O(1)$
, and the third term may be rewritten in terms of the freezing range,
${\rm\Delta}T_{m}=(T_{liq}-T_{sol})$
, and the Stefan number, St. The scaled energy equation can then be simplified:
$$\begin{eqnarray}\frac{{\rm\Delta}T_{sol}}{t}\sim \frac{{\it\alpha}{\rm\Delta}T_{sol}}{{\it\delta}_{m1}^{2}},\left(\frac{1}{\mathit{St}}\right)\frac{{\rm\Delta}T_{m}}{t},\quad \text{where}~\mathit{St}=\frac{c_{p}{\rm\Delta}T_{m}}{L_{f}}.\end{eqnarray}$$
The thickness of the mushy zone,
${\it\delta}_{m1}$
, may be found by rearranging (3.4):
$$\begin{eqnarray}{\it\delta}_{m1}\sim \sqrt{Bt},\quad \text{where }B=\frac{{\it\alpha}}{\displaystyle 1+\frac{{\rm\Delta}T_{m}}{{\rm\Delta}T_{sol}}\frac{1}{\mathit{St}}}.\end{eqnarray}$$
The result in (3.5) is characteristic of conduction-dominated problems, giving a square root relationship with time. Here, however, it is seen through the parameter
$B$
that the thermal diffusivity is augmented by the superheat in
${\rm\Delta}T_{sol}$
, and by the latent heat of fusion within the Stefan number. The inclusion of these parameters in (3.5) marks an improvement over previous scaling analyses for alloy solidification (Krane & Incropera Reference Krane and Incropera1996; Amberg Reference Amberg1997). A relationship between the thicknesses of the solid and mushy regions may be found by writing a flux balance at the solid–mush interface in which linear temperature gradients are assumed for both regions:
$$\begin{eqnarray}\frac{k(T_{avg}-T_{sol})}{{\it\delta}_{m1}+\bar{{\it\delta}}_{T}}\sim \frac{k(T_{sol}-T_{w})}{{\it\delta}_{s1}}\end{eqnarray}$$
or
$$\begin{eqnarray}{\it\delta}_{s1}\sim {\it\delta}_{m1}\left(\frac{T_{sol}-T_{w}}{T_{avg}-T_{sol}}\right),\end{eqnarray}$$
where
$k$
is the thermal conductivity and
${\it\delta}_{s1}$
is the thickness of the solid region during the first solidification regime. The location of the solidification front,
${\it\delta}$
, for the first temporal regime in which all three spatial regions are present is now simply the sum of the solid and mush thicknesses:
$$\begin{eqnarray}{\it\delta}\sim {\it\delta}_{m1}+{\it\delta}_{s1}\sim \sqrt{Bt}\left(\frac{T_{sol}-T_{w}}{{\rm\Delta}T_{sol}}+1\right).\end{eqnarray}$$
The flux balance at the solid–mush interface (3.6) adds an additional relationship with the wall temperature, so that the position of the liquidus front (
${\it\delta}$
) is a function of the superheat, latent heat of fusion and the boundary condition. The time,
$t_{m}(x)$
, at which the liquidus interface passes a particular position,
$x$
, may be found by replacing
${\it\delta}$
with the position of interest,
$$\begin{eqnarray}t_{m}(x)\sim \frac{x^{2}}{B}\left[1+\frac{T_{sol}-T_{w}}{{\rm\Delta}T_{sol}}\right]^{-2}.\end{eqnarray}$$
The time at which the domain transitions from the first to the second temporal regime, which includes only solid and mush, may be approximated by evaluating (3.9) at the cavity length,
$L$
:
$$\begin{eqnarray}t_{m}(L)\sim \frac{L^{2}}{B}\left[1+\frac{T_{sol}-T_{w}}{{\rm\Delta}T_{sol}}\right]^{-2}.\end{eqnarray}$$
The time at which the solidus interface passes a particular position during the first time regime may be found by solving (3.7) for time, denoted
$t_{s1}(x)$
:
$$\begin{eqnarray}t_{s1}(x)\sim \frac{x^{2}}{B}\left[\frac{T_{sol}-T_{w}}{{\rm\Delta}T_{sol}}\right]^{-2}.\end{eqnarray}$$
After the liquid has been consumed by the advancing solidification front, the scaling analysis must be repeated, considering only the solid and mush regions. Scaling (3.1) over the mush gives
$$\begin{eqnarray}\frac{{\rm\Delta}T_{L}}{t_{s2}(L)-t}\sim \frac{{\it\alpha}{\rm\Delta}T_{L}}{{\it\delta}_{m2}^{2}},\frac{-L_{f}{\rm\Delta}f_{l}}{c_{p}(t_{s2}(L)-t)},\end{eqnarray}$$
where
${\rm\Delta}T_{L}=T_{L}(t)-T_{sol}$
is the temperature range bounding the mushy region,
$T_{L}$
is the unknown temperature of the insulated wall at
$x=L$
,
$t_{s2}(L)$
is the time at which the solidus interface reaches the length of the cavity,
${\it\delta}_{m2}$
is the thickness of the mushy zone in the second temporal regime and
${\rm\Delta}f_{l}$
is the reference difference for mass fraction liquid, which is used here, noting that the time derivative of mass fraction liquid is equal to the negative time derivative of mass fraction solid. If a linear change in fraction liquid with temperature is assumed,
$$\begin{eqnarray}{\rm\Delta}f_{l}\sim \frac{{\rm\Delta}T_{L}}{{\rm\Delta}T_{m}},\end{eqnarray}$$
then substituting (3.13) into (3.12) produces an expression for
${\it\delta}_{m2}$
:
$$\begin{eqnarray}{\it\delta}_{m2}\sim \left[\frac{{\it\alpha}(t_{s2}(L)-t)}{1+1/\mathit{St}}\right]^{1/2}.\end{eqnarray}$$
The total solidification time,
$t_{s2}(L)$
, can be found by using the end of the first temporal regime as an initial condition, that is, setting
${\it\delta}_{m2}$
equal to
${\it\delta}_{m1}$
at
$t_{m}(L)$
:
$$\begin{eqnarray}t_{s2}(L)\sim ({\it\delta}_{m1}(t_{m}(L)))^{2}\frac{(1+1/\mathit{St})}{{\it\alpha}}+t_{m}(L).\end{eqnarray}$$
The location of the solidus isotherm at a given time during the second solidification regime,
${\it\delta}_{s2}$
, is simply (3.14) subtracted from the length of the cavity:
$$\begin{eqnarray}{\it\delta}_{s2}\sim (L-{\it\delta}_{m2})\sim L-\left[\frac{{\it\alpha}(t_{s2}(L)-t)}{\displaystyle 1+\frac{1}{\mathit{St}}}\right]^{1/2}.\end{eqnarray}$$
The time,
$t_{s2}(x)$
, at which the solidus isotherm passes a particular location during the second time regime can be expressed by solving for time in (3.16) and replacing
${\it\delta}_{s2}$
with any location
$x$
:
$$\begin{eqnarray}t_{s2}(x)\sim \left(\frac{1+1/\mathit{St}}{{\it\alpha}}\right)({{\it\delta}_{m1}(t_{m}(L))}^{2}-(L-x)^{2})+t_{m}(L).\end{eqnarray}$$
The LST as a function of position,
$t_{lst}(x)$
, can be found from the difference in time between when the liquidus and solidus isotherms pass a particular point using (3.9), (3.11) and (3.17). The time at which the solidus interfaces passes a point must be split between the two solidification regimes. The location at which this occurs may be found from (3.7) by substituting in the time of the transition:
$$\begin{eqnarray}{\it\delta}_{s1}(t_{m}(L))\sim \left(\frac{T_{sol}-T_{w}}{{\rm\Delta}T_{sol}}\right)\sqrt{Bt_{m}(L)}.\end{eqnarray}$$
The LST can now be found for each regime by subtracting
$t_{m}(x)$
from
$t_{s1}(x)$
or
$t_{s2}(x)$
as appropriate:
$$\begin{eqnarray}\displaystyle \displaystyle \text{if }x<{\it\delta}_{s1}(t_{m}(L)):\quad t_{lst} & {\sim} & \displaystyle \frac{x^{2}}{B}\left[\left(\frac{{\rm\Delta}T_{sol}}{T_{sol}-T_{w}}\right)^{2}-\left(1+\frac{T_{sol}-T_{w}}{{\rm\Delta}T_{sol}}\right)^{-2}\right],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \text{if }x>{\it\delta}_{s1}(t_{m}(L)):\quad t_{lst} & {\sim} & \displaystyle \left[\left(\frac{1+1/\mathit{St}}{{\it\alpha}}\right)({{\it\delta}_{m1}(t_{m}(L))}^{2}-(L-x)^{2})+t_{m}(L)\right]\nonumber\\ \displaystyle & & \displaystyle -\left[\frac{x^{2}}{B}\left(1+\frac{T_{sol}-T_{w}}{{\rm\Delta}T_{sol}}\right)^{-2}\right].\end{eqnarray}$$
The LST monotonically increases with position during the first solidification regime (3.19a ), which is to be expected, as conduction to the cooling wall is occurring over ever longer distances. However, there are two competing terms during the second solidification regime (3.19b ), suggesting that, in this portion of the domain, there is a maximum in the LST. This effect is the result of the competition between increasing conduction resistance, slowing the solidification front as it conducts over a longer distance, and decreasing latent heat release rate as the width of the mushy zone is reduced.
3.2. Flow scaling
The scaling of the fluid flow in the cavity has proceeded under the following assumptions. First, it has been assumed that fluid does not penetrate or leave the mushy zone, which may be considered perfectly impermeable, and the ‘flat interface’ assumption (Woods & Huppert Reference Woods and Huppert1989) is used for the boundary between the mush and liquid such that the solid is assumed to be locally flat. The flow is also assumed to be thermally dominated, and composition effects are neglected. These assumptions allow flow down the solidification interface to be approximated with a quasi-static flat-plate boundary layer solution for a constant wall temperature where, for
$t>0$
,
$T(x={\it\delta})\sim T_{liq}$
. This boundary layer flows down the solidification front and fills the cavity with a fluid that is cooler than what is entrained from the bulk. The consequences of these assumptions are discussed in § 7.
3.2.1. Scaling analysis of the transient leading-edge effect propagation
The scaling procedure for the propagation of the transient leading-edge effect of the boundary layer follows that outlined by Lin et al. (Reference Lin, Armfield and Patterson2007) and begins by considering the one-dimensional transient conduction regime ahead of the leading-edge effect in the fluid immediately bordering the solidification interface (Schetz & Eichhorn Reference Schetz and Eichhorn1962). The thermal boundary layer thickness,
${\it\delta}_{T}$
, may be found for the transient conduction regime by scaling the simplified energy equation, (3.1), while neglecting the latent heat term:
$$\begin{eqnarray}{\it\delta}_{T}\sim ({\it\alpha}t)^{1/2}.\end{eqnarray}$$
For the one-dimensional conduction regime, the
$x$
direction velocity and any
$y$
derivatives are zero. The vertical momentum equation then reduces to
$$\begin{eqnarray}\frac{\partial v}{\partial t}={\it\nu}\frac{\partial ^{2}v}{\partial x^{2}}+g{\it\beta}_{T}{\rm\Delta}T,\end{eqnarray}$$
where
$v$
is the vertical velocity component,
${\it\nu}$
is the kinematic viscosity,
$g$
is the acceleration due to gravity and
${\it\beta}_{T}$
is the coefficient of thermal expansion. Equation (3.21) is scaled using
$v_{0}$
as the unknown velocity reference and the thermal boundary layer thickness as the
$x$
direction length scale:
$$\begin{eqnarray}\begin{array}{@{}lcc@{}}\displaystyle \frac{v_{0}}{t}\sim & \displaystyle \frac{{\it\nu}v_{0}}{{\it\delta}_{T}^{2}}, & g{\it\beta}_{T}{\rm\Delta}T_{liq},\\ \,I & F & B\end{array}\end{eqnarray}$$
where
${\rm\Delta}T_{liq}=T_{avg}-T_{liq}$
is the temperature difference across the boundary layer causing the thermal buoyancy. The first term in (3.22) is inertia (
$I$
), the second is viscous friction (
$F$
) and the third is buoyancy (
$B$
). The ratio of the inertia and friction terms may be used to select the appropriate balance with the buoyancy term:
$$\begin{eqnarray}\frac{I}{F}\sim \frac{{\it\delta}_{T}^{2}}{t{\it\nu}}\sim \frac{1}{\mathit{Pr}},\end{eqnarray}$$
where the Prandtl number
$\mathit{Pr}={\it\nu}/{\it\alpha}$
. For the case of liquid metals, the Prandtl number is much less than one, so (3.23) shows that the friction term may be neglected relative to inertia. Balancing the buoyancy and inertial forces in (3.22) leads to a scale for the vertical velocity:
$$\begin{eqnarray}v_{0}\sim g{\it\beta}{\rm\Delta}T_{liq}t.\end{eqnarray}$$
Behind the leading-edge effect, heat transfer is dominated by the vertical advection term. In this region, the energy equation may be simplified by neglecting the horizontal advection and conduction terms:
$$\begin{eqnarray}\frac{\partial T}{\partial t}+v\frac{\partial T}{\partial y}=0.\end{eqnarray}$$
Scaling (3.25) over the thermal boundary layer gives a relationship between the time, position and vertical velocity scales:
$$\begin{eqnarray}\frac{v_{0}{\rm\Delta}T_{liq}}{y}\sim \frac{{\rm\Delta}T_{liq}}{t}\end{eqnarray}$$
and
$$\begin{eqnarray}y\sim v_{0}t\sim (g{\it\beta}{\rm\Delta}T_{liq})t^{2}.\end{eqnarray}$$
This relationship (3.26b
) may be interpreted as the time at which the leading-edge effect reaches a particular location,
$y$
. By substituting in the velocity scale (3.24) found from the vertical momentum equation, a relationship between only position and time is found. The main quantity of interest is the time required for the boundary layer to develop over the full height of the cavity:
$$\begin{eqnarray}t_{H}\sim \left[\frac{H}{g{\it\beta}{\rm\Delta}T_{liq}}\right]^{1/2}.\end{eqnarray}$$
Comparing (3.27) to (3.10) shows that the assumption that the leading-edge effect propagates to the height of the cavity in a much shorter time than it takes for the metal to solidify is correct if
$$\begin{eqnarray}\frac{t_{H}}{t_{m}(L)}\ll 1.\end{eqnarray}$$
Substituting order-of-magnitude values into (3.28) for a metallic system with a moderate superheat shows that this ratio is approximately
$O(10^{-5})$
. Therefore, the inequality in (3.28) is satisfied and the use of a quasi-static boundary layer solution for the flow at the mush–liquid interface throughout the course of the process is reasonable.
3.2.2. Quasi-static boundary layer scaling analysis
The solution for the flow in the bulk liquid near the liquidus surface was inspired by the work of Lin et al. (Reference Lin, Armfield and Patterson2007). The quasi-static boundary layer solution follows from Bejan (Reference Bejan1995) and starts by considering the steady energy equation, where the vertical advection term is retained and conduction along the direction of flow is neglected. The equation is scaled using
$u_{0}$
and
$v_{0}$
as the velocity scales in the horizontal and vertical directions, and
${\it\delta}_{T}$
and
$H$
as the
$x$
and
$y$
direction length scales, respectively. The height of the cavity is the appropriate vertical length scale, since the boundary layer flow occurs over the full height of the cavity throughout the process, as shown by the streamlines in figure 3 (this vertical length scale is an important departure from Lin et al. (Reference Lin, Armfield and Patterson2007), as will be discussed in § 7):
$$\begin{eqnarray}\displaystyle & \displaystyle u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}={\it\alpha}\frac{\partial ^{2}T}{\partial x^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{0}\frac{{\rm\Delta}T_{liq}}{{\it\delta}_{T}},v_{0}\frac{{\rm\Delta}T_{liq}}{H}\sim {\it\alpha}\frac{{\rm\Delta}T_{liq}}{{\it\delta}_{T}^{2}}. & \displaystyle\end{eqnarray}$$
The
$x$
and
$y$
velocity scales may be related by scaling the incompressible continuity equation and using the thermal boundary layer thickness and cavity height as the
$x$
and
$y$
length scales, respectively:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{0}\sim \frac{v_{0}{\it\delta}_{T}}{H}. & \displaystyle\end{eqnarray}$$
Substituting (3.30b
) into (3.29b
) shows that the two advection terms are equal, and allows for a relationship between the vertical velocity scale,
$v_{0}$
, and the thermal boundary layer thickness,
${\it\delta}_{T}$
:
$$\begin{eqnarray}v_{0}\sim \frac{{\it\alpha}H}{{\it\delta}_{T}^{2}}.\end{eqnarray}$$
The steady vertical momentum equation is now considered, in which the transient, pressure and vertical friction terms have been neglected. The two advection terms are shown to be equal by employing (3.31):
$$\begin{eqnarray}\displaystyle & \displaystyle u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}={\it\nu}\frac{\partial ^{2}v}{\partial x^{2}}+g{\it\beta}_{T}(T-T_{avg}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \begin{array}{@{}lcc@{}}\displaystyle \frac{v_{0}^{2}}{H}\sim & \displaystyle \frac{{\it\nu}v_{0}}{{\it\delta}_{T}^{2}}, & g{\it\beta}_{T}{\rm\Delta}T_{liq}.\\ \,I & F & B\end{array} & \displaystyle\end{eqnarray}$$
The first term represents inertia, the second represents viscous friction forces and the third is the driving buoyancy force. Equation (3.31) is used to eliminate
${\it\delta}_{T}$
from the viscous term, yielding
$$\begin{eqnarray}\frac{v_{0}^{2}}{H},\mathit{Pr}\frac{v_{0}^{2}}{H}\sim g{\it\beta}_{T}{\rm\Delta}T_{liq}.\end{eqnarray}$$
In (3.33) the buoyancy term may balance with either the inertial or viscous term, the correct choice being determined by the magnitude of Prandtl number. For liquid metals, the Prandtl number is much less than one, and so the friction term may be neglected relative to inertia. Thus, a scale for the velocity in the boundary layer at the bottom of the cavity is
$$\begin{eqnarray}v_{0}\sim [g{\it\beta}_{T}{\rm\Delta}T_{liq}H]^{1/2}.\end{eqnarray}$$
Substitution of (3.34) into (3.31) yields a scale for the thermal boundary layer thickness,
$$\begin{eqnarray}{\it\delta}_{T}\sim \left[\frac{H{\it\alpha}^{2}}{g{\it\beta}_{T}{\rm\Delta}T_{liq}}\right]^{1/4}.\end{eqnarray}$$
The fluid flowing down the mush–liquid interface is estimated to exit into the bulk liquid at the bottom of the cavity at the boundary layer film temperature,
$T_{BL}(t)$
,
$$\begin{eqnarray}T_{BL}(t)\approx \frac{T_{avg}(t)+T_{liq}}{2}.\end{eqnarray}$$
So far, it has been assumed that the average temperature of the bulk liquid is constant over time, but, of course, as fluid near the liquidus temperature is advected from the boundary layer in exchange for the entrainment of warmer fluid, this temperature continuously decreases. In order to estimate the average temperature of the fluid over time, the enthalpy of the liquid is considered. It remains useful to approximate the liquid as a reservoir with a uniform temperature for the solidification and boundary layer solutions, but this assumption is not appropriate for determining the enthalpy of the fluid that is entrained into the boundary layer. It is clear from figure 2 that the bulk liquid quickly becomes thermally stratified and, from the streamlines in figure 3, that liquid moving through the boundary layer exits at the bottom of the cavity and is mostly entrained at the top. Therefore, the enthalpy of the liquid entrained into the boundary layer is found from the temperature at the top of the cavity at a given time, necessitating an approximation of the thermal stratification of the system.
The cavity begins as all liquid at the initial superheat,
$T_{avg}(t=0)=T_{0}$
. During the initial cool-down period, this volume at
$T_{0}$
is consumed at a rate equal to the flow rate out of the boundary layer. It is assumed that a discrete interface at
$y=Y_{i}$
forms between fluid that is at the initial superheat and fluid that has been advected out of the boundary layer, as shown in figure 4(a). The temperature of the lower region is equal to the volume-weighted average of the fluid that has exited the boundary layer at
$T_{BL}(t)$
and is further assumed to have a linear vertical temperature profile with a minimum at the liquidus temperature. The fluid closest to the mushy zone is very close to the liquidus temperature, leaving the vertical boundary layer at
$y=H$
and moving along the bottom of the cavity. The total average temperature of the liquid (
$T_{avg}$
) is then computed as an average of the temperature of the lower region and the initial superheat, weighted by
$Y_{i}$
.
Figure 4. Assumed temperature profiles in the bulk liquid (a) before and (b) after the initial superheat has been consumed.
The distance
$Y_{i}$
decreases as fluid from the boundary layer continues to displace the liquid at the initial temperature. At some time,
$Y_{i}$
is equal to zero and the initial superheat has been consumed. After this time, it is assumed that the fluid is stratified such that the temperature profile is linear, with a minimum at the liquidus temperature and the average temperature at
$H/2$
, as shown in figure 4(b). Fluid at
$T_{BL}(t)$
continues to enter the cavity at a rate determined by the boundary layer flow analysis, now at the expense of fluid at a temperature averaged over an equal volume at the top of the cavity, denoted
$\bar{T}_{max}$
. These approximations for the thermal stratification in the bulk liquid are in contrast to its treatment by Lin et al. (Reference Lin, Armfield and Patterson2007), who assume two isothermal regions and that all fluid exiting from the boundary layer is at the wall temperature. The effect of these changes will be discussed in § 7.
Considering the enthalpy of the fluid being advected into the cavity, and that being consumed by entrainment into the boundary layer and by the advancing solidification front yields equations for the average liquid temperature
$T_{avg}(t)$
in the two temporal regimes shown in figure 4, in which the specific heat has been assumed constant with temperature over the liquid region:
$$\begin{eqnarray}\displaystyle & \displaystyle \text{if }Y_{i}>0:\quad H(L-{\it\delta})\frac{\text{d}T_{avg}}{\text{d}t}\sim -\mathit{HT}_{avg}\frac{\text{d}{\it\delta}}{\text{d}t}+v_{0}{\it\delta}_{T}(T_{BL}-T_{0}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{if }Y_{i}=0:\quad H(L-{\it\delta})\frac{\text{d}T_{avg}}{\text{d}t}\sim -\mathit{HT}_{avg}\frac{\text{d}{\it\delta}}{\text{d}t}+v_{0}{\it\delta}_{T}(T_{BL}-\bar{T}_{max}). & \displaystyle\end{eqnarray}$$
In both cases, the first term represents the storage of energy in the bulk liquid and the second term is the reduction in the amount of liquid at the average temperature due to the advancing solidification front. The third term is the cooling of the fluid in the boundary layer from the temperature at which it is entrained (either
$T_{0}$
or
$\bar{T}_{max}$
) to that at which it leaves (
$T_{BL}$
). Equations (3.37) are complicated by the dependence of the boundary layer velocity and thickness on the average fluid temperature
$T_{avg}$
, but can be readily solved numerically. While the previous scaling analysis for the one-dimensional solidification assumed that the superheat was a constant, it is clear from the above discussion that it decreases continuously at a rate determined by the bulk fluid flow. The solidification scaling can be improved by using the transient average temperature for the superheat in (3.8) for
${\it\delta}$
. To account for the behaviour of
$T_{avg}(t)$
found in (3.37),
${\it\delta}(t+{\rm\Delta}t)$
is found by using a truncated Taylor series expansion of
${\it\delta}(t)$
,
$$\begin{eqnarray}{\it\delta}(t+{\rm\Delta}t)\sim {\it\delta}(t)+{\rm\Delta}t\left(\frac{\partial {\it\delta}(T_{avg}(t))}{\partial t}\right)_{{\it\delta}(t)}.\end{eqnarray}$$
Considering the effect of the changing superheat on the solidification adds an additional nonlinearity to (3.37), which should now be solved iteratively with (3.38).
4. Integral analysis
A more detailed solution to this solidification problem can be found by using an integral analysis of the energy equation in the mushy zone. For the first solidification regime, the temperature profile is assumed to be linear in the solid, quadratic in the mush and linear over the thermal boundary layer. In the second solidification regime, the same profiles are used for the solid and mush, in which the mush temperature profile terminates at the far wall at an unknown temperature,
$T_{L}$
, between the solidus and liquidus temperature of the alloy. These configurations are shown in figure 5.
Figure 5. A schematic of the temperature profiles over the solid, mush, boundary layer and liquid regions used for the integral solution during (a) the first solidification regime and (b) the second.
Under the same assumptions used for the solidification scaling, the integral analysis begins by constructing non-dimensional spatial and temperature variables for the mushy zone and assuming a quadratic profile for the non-dimensional temperature profile in the mushy zone:
$$\begin{eqnarray}{\it\theta}_{m}({\it\eta})=\frac{T-T_{sol}}{T_{liq}-T_{sol}}=C_{1}{\it\eta}^{2}+C_{2}{\it\eta}+C_{3},\quad \text{where}~{\it\eta}=\frac{x-{\it\delta}_{s1}}{{\it\delta}_{m1}},\end{eqnarray}$$
where
${\it\theta}_{m}$
is the non-dimensional temperature in the mushy zone,
${\it\eta}$
is the non-dimensional distance, and
$C_{1}$
,
$C_{2}$
and
$C_{3}$
are unknown constants in the assumed polynomial profile. The mushy zone is bounded by the solidus and liquidus temperatures. These boundary conditions can be expressed in terms of the non-dimensional variables in (4.1):
$$\begin{eqnarray}\displaystyle & {\it\theta}_{m}(0)=0=C_{3}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\theta}_{m}(1)=1=C_{1}+C_{2}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\left.\frac{\partial {\it\theta}_{m}}{\partial {\it\eta}}\right|_{{\it\eta}=1^{-}}=\left.\frac{\partial {\it\theta}_{m}}{\partial {\it\eta}}\right|_{{\it\eta}=1^{+}}.\end{eqnarray}$$
The left-hand side of the equation can be easily evaluated using the quadratic profile in (4.1). The right-hand side of the equation can be expanded using the chain rule:
$$\begin{eqnarray}\left.\frac{\partial {\it\theta}_{m}}{\partial {\it\eta}}\right|_{{\it\eta}=1^{-}}=\frac{\partial {\it\theta}_{m}}{\partial T}\frac{\partial T}{\partial x}\left.\frac{\partial x}{\partial {\it\eta}}\right|_{{\it\eta}=1^{+}}.\end{eqnarray}$$
It is assumed that the temperature gradient in the liquid is linear over some mean boundary layer thickness that is constant in
$y$
, denoted
$\bar{{\it\delta}}_{T}$
:
$$\begin{eqnarray}\left.\frac{\partial T}{\partial x}\right|_{{\it\eta}=1^{+}}=\frac{{\rm\Delta}T_{liq}}{\bar{{\it\delta}}_{T}}.\end{eqnarray}$$
The value of the mean boundary layer thickness is where the bulk fluid flow influences the solidification behaviour and can be obtained from the scaling solution or from a more sophisticated analysis. Evaluating (4.3) allows for the remaining coefficients in the equation for the temperature profile for the mushy zone to be found:
$$\begin{eqnarray}{\it\theta}_{m}({\it\eta})=(P{\it\delta}_{m1}-1){\it\eta}^{2}+(2-P{\it\delta}_{m1}){\it\eta},\quad \text{where}~P=\frac{{\rm\Delta}T_{liq}}{\bar{{\it\delta}}_{T}{\rm\Delta}T_{m}}.\end{eqnarray}$$
The governing one-dimensional energy equation for solidification when advection is neglected is shown in (3.1). The equation can be simplified by assuming a linear relationship between fraction solid and temperature:
$$\begin{eqnarray}\frac{\partial f_{s}}{\partial t}=\frac{\partial f_{s}}{\partial T}\frac{\partial T}{\partial t}=\frac{-1}{{\rm\Delta}T_{m}}\frac{\partial T}{\partial t}.\end{eqnarray}$$
The simplified energy equation can then be written as
$$\begin{eqnarray}\frac{\partial T}{\partial t}={\it\alpha}^{\prime }\frac{\partial ^{2}T}{\partial x^{2}},\quad \text{where }{\it\alpha}^{\prime }={\it\alpha}\left(\frac{\mathit{St}}{1+\mathit{St}}\right).\end{eqnarray}$$
Changing variables in (4.7) and integrating over
${\it\eta}$
gives
$$\begin{eqnarray}\int _{0}^{1}\frac{\partial {\it\theta}_{m}}{\partial {\it\eta}}\left[\frac{\partial {\it\eta}}{\partial {\it\delta}_{m1}}\frac{\partial {\it\delta}_{m1}}{\partial t}+\frac{\partial {\it\eta}}{\partial {\it\delta}_{s1}}\frac{\partial {\it\delta}_{s1}}{\partial t}\right]\,\text{d}{\it\eta}=\frac{{\it\alpha}^{\prime }}{{\it\delta}_{m1}^{2}}\int _{0}^{1}\frac{\partial ^{2}{\it\theta}_{m}}{\partial {\it\eta}^{2}}\,\text{d}{\it\eta}\end{eqnarray}$$
and
$$\begin{eqnarray}\dot{{\it\delta}}_{m1}\left[\frac{-2-P{\it\delta}_{m1}}{6{\it\delta}_{m1}}\right]-\frac{\dot{{\it\delta}}_{s1}}{{\it\delta}_{m1}}=\frac{{\it\alpha}^{\prime }}{{\it\delta}_{m1}^{2}}(2P{\it\delta}_{m1}-2).\end{eqnarray}$$
An additional relationship between the solid and mushy zone thicknesses is needed. This function can be constructed from a flux balance at the solid–mush interface:
$$\begin{eqnarray}\left.-k\frac{\partial T}{\partial x}\right|_{{\it\delta}_{s1}^{-}}=\left.-k\frac{\partial T}{\partial x}\right|_{{\it\delta}_{s1}^{+}}.\end{eqnarray}$$
A linear temperature profile between the solidus and wall temperatures has been assumed for the solid. The chain rule and the non-dimensional equation for the temperature in the mushy zone are used to evaluate the temperature gradient in the mush:
$$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\partial T}{\partial x}\right|_{{\it\delta}_{s1}^{-}}=\frac{T_{sol}-T_{w}}{{\it\delta}_{s}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\partial T}{\partial x}\right|_{{\it\delta}_{s1}^{+}}=\left.\frac{\partial T}{\partial {\it\theta}_{m}}\frac{\partial {\it\theta}_{m}}{\partial {\it\eta}}\frac{\partial {\it\eta}}{\partial x}\right|_{{\it\eta}=0^{+}}=\frac{{\rm\Delta}T_{m}(2-P{\it\delta}_{m})}{{\it\delta}_{m}}. & \displaystyle\end{eqnarray}$$
Substituting (4.11) into (4.10) yields the relationship between the solid and mushy zone thicknesses:
$$\begin{eqnarray}{\it\delta}_{s1}=(T_{sol}-T_{w})\left[\frac{2{\rm\Delta}T_{m}}{{\it\delta}_{m1}}-P{\rm\Delta}T_{m}\right]^{-1}.\end{eqnarray}$$
Equation (4.12) can be substituted into (4.9) to yield a first-order differential equation for the mushy zone thickness:
$$\begin{eqnarray}\displaystyle \frac{\partial {\it\delta}_{m1}}{\partial t}=\left[\frac{{\it\alpha}^{\prime }}{{\it\delta}_{m1}^{2}}(2P{\it\delta}_{m1}-2)\right]\!\left[\frac{-2-P{\it\delta}_{m1}}{6{\it\delta}_{m1}}-\frac{1}{{\it\delta}_{m1}}\left(\frac{T_{sol}-T_{w}}{{\rm\Delta}T_{m}(2-P{\it\delta}_{m1})}\right)\!\left(\frac{P{\it\delta}_{m1}}{(2-P{\it\delta}_{m1})}+1\right)\!\right]^{-1}.\hspace{-10.00002pt} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
This analysis has assumed that the average liquid temperature is constant with time. As noted in the previous section, the location of the solidification interface is coupled to the flow because the bulk fluid is cooled by the boundary layer as dictated by (3.37). To account for this nonlinearity, equation (4.13) can be combined with (4.12) to numerically predict the solidification front velocity for a given solidification front location and instantaneous average liquid temperature by noting that the location of the solidification front is simply the sum of the solid and mushy zone thicknesses. This approximation for the front velocity can be used to predict the subsequent position of the solidification front, although this procedure must be performed iteratively to account for the changing average bulk fluid temperature. In the same way as for the scaling analysis, this procedure is predicated on the observation that the front velocity is quasi-static and so is independent of the path by which the system approached the current values of the solid and mushy zone thicknesses, effective superheat and wall temperature. The instantaneous average fluid temperature is evaluated as described in § 3.2.2.
Once the liquid is consumed and the mush–liquid interface reaches the far wall (
$x=L$
), the second temporal regime begins, and the thermal problem changes to one in which the temperature of the opposing vertical wall is unknown and the heat flux there is zero. The new non-dimensional variables over the mushy zone are formulated as before and a quadratic profile is assumed in the mush:
$$\begin{eqnarray}{\it\theta}_{m}^{\ast }({\it\varepsilon})=\frac{T-T_{sol}}{T_{liq}-T_{sol}}=C_{4}{\it\varepsilon}^{2}+C_{5}{\it\varepsilon}+C_{6},\quad \text{where }{\it\varepsilon}=\frac{x-{\it\delta}_{s2}}{{\it\delta}_{m2}}.\end{eqnarray}$$
Here
${\it\theta}_{m}^{\ast }$
is the non-dimensional temperature profile for the mushy zone,
${\it\varepsilon}$
is the non-dimensional position during the second solidification regime, and
$C_{4}$
,
$C_{5}$
and
$C_{6}$
are constant polynomial coefficients for the mushy zone temperature profile. The boundary conditions at either end of the mushy zone are as follows:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\it\theta}_{m}^{\ast }({\it\varepsilon}=0)=0,\\ \displaystyle {\it\theta}_{m}^{\ast }({\it\varepsilon}=1)=\frac{T_{L}-T_{sol}}{T_{liq}-T_{sol}}={\it\theta}_{L}(t),\\ \displaystyle \left.\frac{\partial {\it\theta}_{m}^{\ast }}{\partial {\it\varepsilon}}\right|_{{\it\varepsilon}=1}=0,\end{array}\right\}\end{eqnarray}$$
where
${\it\theta}_{L}$
is the unknown non-dimensional temperature at
$x=L$
. The resulting equation for the non-dimensional temperature is then
$$\begin{eqnarray}{\it\theta}_{m}^{\ast }={\it\theta}_{L}(2{\it\varepsilon}-{\it\varepsilon}^{2}).\end{eqnarray}$$
The simplified governing equation in (4.7) can be transformed to the non-dimensional variables of (4.14) and integrated over the remaining mushy zone:
$$\begin{eqnarray}\int _{0}^{1}\frac{\partial T}{\partial {\it\theta}_{m}^{\ast }}\left[\frac{\partial {\it\theta}_{m}^{\ast }}{\partial {\it\varepsilon}}\frac{\partial {\it\varepsilon}}{\partial {\it\delta}_{m2}}\frac{\partial {\it\delta}_{m2}}{\partial t}+\frac{\partial {\it\theta}_{m}^{\ast }}{\partial {\it\theta}_{L}}\frac{\partial {\it\theta}_{L}}{\partial t}\right]\text{d}{\it\varepsilon}=\int _{0}^{1}\frac{\partial }{\partial {\it\varepsilon}}\left[\frac{\partial T}{\partial {\it\theta}_{m}^{\ast }}\frac{\partial {\it\theta}_{m}^{\ast }}{\partial {\it\varepsilon}}\left(\frac{\partial {\it\varepsilon}}{\partial x}\right)^{2}\right]\text{d}{\it\varepsilon}.\end{eqnarray}$$
Substituting for all of the derivatives and carrying out the integration yields a differential equation for the mushy zone thickness as a function of the unknown wall temperature at
$L$
and its time derivative:
$$\begin{eqnarray}\dot{{\it\delta}}_{m2}={\it\delta}_{m2}\left[\frac{-3{\it\alpha}^{\prime }}{{\it\delta}_{m2}^{2}}-\frac{\dot{{\it\theta}}_{L}}{{\it\theta}_{L}}\right].\end{eqnarray}$$
A relationship between the mushy zone thickness and the wall temperature at
$L$
can be found by equating the fluxes at the solid–mush interface, which begins as shown in (4.10). Changing variables and simplifying by noting that the sum of the solid and mush thicknesses is the length of the domain (
${\it\delta}_{m2}+{\it\delta}_{s2}=L$
) yields
$$\begin{eqnarray}{\it\theta}_{L}=\left(\frac{T_{sol}-T_{w}}{2{\rm\Delta}T_{m}}\right)\left(\displaystyle \frac{1}{\displaystyle \frac{L}{{\it\delta}_{m2}}-1}\right).\end{eqnarray}$$
Substituting (4.19) into (4.18) results in an ordinary differential equation for the thickness of the mushy zone:
$$\begin{eqnarray}\dot{{\it\delta}}_{m2}=\frac{-3{\it\alpha}^{\prime }}{{\it\delta}_{m2}\left(\displaystyle 1+\frac{L}{L-{\it\delta}_{m2}}\right)}.\end{eqnarray}$$
Here, the inverse dependence of the front velocity on the mushy zone thickness suggests that the motion of the solid–mush interface will accelerate as the mush thickness decreases towards the end of solidification. This will tend to reduce the LST, resulting in a maximum in the second solidification regime, consistent with the interpretation of the scaling results. The form of (4.20) can be readily implemented into the numerical framework described for the first temporal regime. However, this equation can also be solved by separating variables and integrating, yielding an implicit expression for the mushy zone thickness as a function of time:
$$\begin{eqnarray}\frac{{\it\delta}_{m2}^{2}}{2}-L{\it\delta}_{m2}-L^{2}\ln (L-{\it\delta}_{m2})=-3{\it\alpha}^{\prime }t+C.\end{eqnarray}$$
The integration constant in (4.21) can be found by using the condition that the thickness of the mushy zone at the time at which the liquidus surface reaches the far wall,
$t_{m}(L)$
, must be equal to the solution for the first temporal regime at that time.
5. Comparison with numerical results
5.1. Numerical model description
The numerical model used to generate predictions for comparison with the above analyses was a continuum mixture model based on the work of Bennon & Incropera (Reference Bennon and Incropera1987a ,Reference Bennon and Incropera b ). In this case, it was assumed that there was no macrosegregation within the casting, such that all of the metal solidified at the nominal composition, and that there are no solutal contributions to buoyancy (an assumption that is discussed in § 7). The material properties are assumed to be independent of temperature and phase fraction. The governing equations for continuity, horizontal and vertical momentum and energy are (Krane Reference Krane2010)
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\bar{\boldsymbol{V}}=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\frac{\partial u}{\partial t}+{\it\rho}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\bar{\boldsymbol{V}}u)={\it\mu}{\rm\nabla}^{2}u-\frac{\partial P}{\partial x}-\frac{{\it\mu}}{K}u, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\frac{\partial v}{\partial t}+{\it\rho}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\bar{\boldsymbol{V}}v)={\it\mu}{\rm\nabla}^{2}v-\frac{\partial P}{\partial y}-\frac{{\it\mu}}{K}v-g{\it\rho}{\it\beta}_{T}(T-T_{ref}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}c_{p}\frac{\partial T}{\partial t}+{\it\rho}c(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }T\bar{\boldsymbol{V}})=k{\rm\nabla}^{2}T-{\it\rho}L_{f}\frac{\partial f_{L}}{\partial t}-{\it\rho}L_{f}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\bar{\boldsymbol{V}}). & \displaystyle\end{eqnarray}$$
Here
$u$
and
$v$
are horizontal and vertical mixture velocity components,
${\it\rho}$
is density,
${\it\mu}$
is dynamic viscosity,
$P$
is pressure,
$T_{ref}$
is a reference temperature and
$\bar{\boldsymbol{V}}$
is the mixture velocity vector,
$$\begin{eqnarray}\bar{\boldsymbol{V}}=f_{s}\boldsymbol{V}_{s}+f_{l}\boldsymbol{V}_{l}.\end{eqnarray}$$
In this case, the solid is assumed stationary, so that the mixture velocity is equal to the last term in (5.5). For the full numerical model, the mushy zone was allowed to be permeable. The last term in (5.2) and the second-to-last term in (5.3) represent the Darcy drag associated with the mushy zone, where
$K$
is the isotropic permeability of the dendritic array determined using the Blake–Kozeny model,
$$\begin{eqnarray}K=\frac{{\it\lambda}^{2}(1-g_{s})^{3}}{180g_{s}^{2}},\end{eqnarray}$$
where
${\it\lambda}$
is the dendrite arm spacing and
$g_{s}$
is the volume fraction solid. The last term in (5.3) is the thermal buoyancy. The final two terms in (5.4) represent the evolution of latent heat during solidification and the advection of latent heat in the liquid, respectively. The transient latent heat term is evaluated using the linearization scheme proposed by Voller & Swaminathan (Reference Voller and Swaminathan1991) and the fraction solid is related to the temperature through a simple binary alloy lever rule. The system of equations (5.1)–(5.6) was solved using standard finite volume techniques and the SIMPLER algorithm on a structured, orthogonal, staggered grid in Cartesian coordinates (Patankar Reference Patankar1980). The general mixture model has been previously validated (Prescott, Incropera & Gaskell Reference Prescott, Incropera and Gaskell1994; Krane & Incropera Reference Krane and Incropera1997) and the temperature formulation of the energy equation was successfully compared to heat balance integral calculations by Voller & Swaminathan (Reference Voller and Swaminathan1991).
Table 1. Property values for Al–4.5 wt% Cu alloy.
Figure 6. Grid dependence study on the base case given in table 2 showing changes in the time at which the system transitions between solidification regimes,
$t_{m}(L)$
, and the total solidification time,
$t_{s2}(L)$
. The central data points (
$50\times 50$
cells on a square 0.25 m domain) was selected as the grid size for use in this study.
Table 2. Test matrix showing system variables for comparison of analyses to simulations.
Property data for Al–4.5 wt% Cu are given in table 1 (Vreeman & Incropera Reference Vreeman and Incropera2000) and the conditions of the simulations run for comparison to the preceding analyses are given in table 2. This table also includes the initial Rayleigh number for each case,
$$\begin{eqnarray}\mathit{Ra}_{H,init}=\frac{g{\it\beta}_{T}(T_{0}-T_{liq})H^{3}}{{\it\alpha}{\it\nu}}.\end{eqnarray}$$
Because the superheat in the bulk decreases during the process, the Rayleigh number relevant to the boundary layer flow is not constant during the process, but initial values are given in table 2 to represent the relative vigour of these flows. For changes in aspect ratio,
$A$
, the length of the cavity was held constant while the height was varied. This approach was used because the solidification for the current system was primarily one-dimensional in the
$x$
direction, and maintaining a constant length allows for direct comparisons of the solidification times. The domain length,
$L$
, for all simulations was 0.25 m. A grid dependence study was performed with the base case given in table 2. The results of various grid sizes for the time at which the domain transitioned between solidification regimes and the total solidification time are shown in figure 6. The central grid spacing of
$5~\text{mm}\times 5~\text{mm}$
, or
$50\times 50$
control volumes in the 0.25 m domain, beyond which little change in the solution was observed, was selected for all subsequent simulations.
5.2. Comparison of the scaling and integral analyses to simulation results
The simplest case for the scaling analysis is one in which the average liquid temperature does not change over time, which reduces the present analysis to the results of Krane & Incropera (Reference Krane and Incropera1996). The location of the solidification front,
${\it\delta}$
, as a function of time predicted by the numerical model was compared to two cases for a wall temperature of 300 K, either with an initial temperature of 1000 K, or with the initial temperature equal to the liquidus temperature (figure 7). Because the scaling results are expected to be correct only within an order of magnitude, they have been multiplied by a constant factor so that
$t_{m}(L)$
for the case with no superheat is equal to that of the simulation, allowing for a direct comparison of their behaviour. Figure 7 shows that the scaling analyses for the maximum and minimum possible liquid temperatures bound the simulation results, with the simulation results ranging from the
$T_{avg}=T_{0}$
case at early times to the
$T_{avg}=T_{liq}$
case later. This behaviour is expected because the average fluid temperature in the simulation decreases from
$T_{0}$
to
$T_{liq}$
throughout the process. At short times, the liquid temperature is approximately equal to the initial superheat, and the scaling analysis with constant non-zero superheat very closely matches the shape of the simulation results. At intermediate times, the fluid temperature is decreasing, causing an effective superheat lower than the initial value. Here, the simulation results depart from the scaling analysis for finite superheat and begins to approach the case with no superheat. Eventually, the fluid has cooled to a temperature near the liquidus of the alloy, and the simulation results begins to match the shape of the curve for scaling results with no superheat.
Figure 7. Comparison of simulation results and scaling results for the location of the solidification front with time. The alloy system is Al–4.5 wt% Cu in a
$0.25~\text{m}\times 0.25~\text{m}$
cavity and a wall temperature at 300 K. The simulation, initially at 1000 K, is compared to scaling results with constant superheat at the initial condition, and for scaling results with no superheat. The scaling results have been dilated such that
$t_{m}(L)$
for the former case is approximately equal to the simulation time.
If the average liquid temperature is allowed to vary over time in accordance with the thermal stratification of the bulk fluid and flow from the boundary layer, it is expected that the scaling analysis will more closely resemble the simulation results. The scaling analysis results for the location of the liquidus front as a function of time are compared in figure 8 with the simulation results. Again, the scaling results have been dilated so that the time at which the liquidus front reaches the opposite wall is equal to the simulation, allowing for a direct comparison. The analogous result for the integral analysis is also shown in figure 8 with no corrections. Both analyses closely predict the shape of
${\it\delta}$
as a function of time, and the integral solution does so without any adjustment relative to the simulation.
Figure 8. Comparison of liquidus position solution to simulation results over time for Al–4.5 wt% Cu alloy in 0.25 m square cavity, initially at 1000 K and cooled with a 300 K sidewall where the scaling results have been dilated so that
$t_{m}(L)$
is equal to the simulation results. The integral analysis results have not been altered.
The purpose of the scaling analysis is to predict the functional form of the relationships between the system and material properties and the solidification parameters. One way to verify that the correct trends were identified in the scaling analysis is to plot the scaling results against the simulation results. If a linear profile is found, then the trends are the same within some constant multiplier.
This approach was used to compare the trends in the scaling and integral analyses with wall temperature and initial temperature to those of the simulation results. The trends for the time required for the solidification front to reach the opposing cavity wall are shown in figure 9. Note that, because of the two-dimensional flow field, the top of the casting solidifies slower than the bottom, as seen in the isotherms in figure 2. Therefore, the solidification times from the simulations are an average over the height of the cavity. Good agreement was found for both analyses for cases where the initial temperature was held constant and the wall temperature was varied. Reasonable agreement was also found when the wall temperature was held constant and the initial temperature was varied, although it can be seen that the change in initial temperature had a smaller effect on
$t_{m}(L)$
than changing the wall temperature.
Figure 9. A comparison of
$t_{m}(L)$
trends for scaling and integral analyses relative to the simulation results for (a) various wall temperatures and (b) various initial temperatures.
Both the scaling and integral analyses are expected to match the trends in the simulations, but the scaling analysis is only expected to agree with these results within an order of magnitude. Figure 10 shows a direct comparison of the time required for the liquidus interface to propagate the length of the cavity as a function of both the wall temperature and initial temperature for the full simulation and the integral analysis. The response to changes in the initial temperature appear to be poorly predicted compared to those for changes in the wall temperature, but in fact are very similar in agreement, and are less sensitive to changes in the particular system parameter. In both cases, the predicted values of
$t_{m}$
are accurate to within approximately 5 % of the simulation results.
Figure 10. Comparison of simulation and integral analysis results for
$t_{m}(L)$
as a function of (a) wall temperature and (b) initial temperature.
Trends for the time required for the solidus interface to reach the opposite wall are shown in figure 11 as a function of both the wall temperature and initial temperature. The comparison to the integral analysis shows slightly more linear trends than to the scaling analysis in this case. The values for the integral analysis are compared directly with the simulation results in figure 12. Note the difference in scale for the vertical axes on the two plots. The integral analysis results are generally accurate within 20 % of the simulation results for all cases, although the quality of the results is comparatively poor for very long solidification times at wall temperatures near the solidus.
Figure 11. Comparison of trends for the scaling and integral analysis for
$t_{s2}(L)$
relative to the simulation results for (a) various wall temperatures and (b) various initial temperatures.
Figure 12. Comparison of
$t_{s}(L)$
between simulation and integral results as a function of (a) wall temperature and (b) initial temperature.
Next, the cool-down time for the bulk liquid is compared to the simulations. For these comparisons, the cool-down time was defined as the time at which the average liquid temperature fell 99 % of the difference between the initial temperature and the liquidus temperature. The results from varying
$T_{w}$
and
$T_{0}$
independently are shown in figure 13. Figure 14 shows the direct comparison between the integral analysis results and the simulations. The integral analysis results are within approximately 8 % of the simulation values.
Figure 13. Comparison of trends in cool-down time with respect to the simulation results for both the scaling integral analyses for (a) various wall temperatures and (b) various initial temperatures.
Figure 14. A direct comparison between the integral analysis and the simulation results for the cool-down time as a function of (a) wall temperature and (b) initial temperature.
In comparing the LST over the variety of cases in table 2, the most pertinent value is the maximum LST, which is important because slower solidification leads to a coarser microstructures, such as larger grain size and spacing of dendrite arms, which subsequently affect the mechanical properties of the alloy. The LST increases with distance from the isothermal vertical wall, but the maximum value does not occur at the insulated vertical wall. Instead, there is a local maximum value, shown in figure 15 for simulation results for the base case, and exhibited by both analyses, although the scaling analysis predicts the maximum LST at the opposing vertical wall for low
$T_{w}$
. This local maximum occurs because, after the liquid–mush interface reaches the opposite wall, further solidification decreases the total latent heat release rate. This decrease allows the solid–mush interface to accelerate, eventually reducing the LST near the insulated vertical wall. The distance over which this decrease in LST occurs is directly related to the thickness of the mushy zone, which is in turn strongly dependent on the wall temperature because this parameter controls the rate at which latent heat is pulled from the mush at the solidus interface. For a low wall temperature and a correspondingly thin mushy zone, the maximum LST occurs near the opposite wall, while for higher wall temperatures and a thicker mushy zone, the maximum LST occurs further from the insulated boundary.
Figure 16(a) shows a comparison of the trends for the scaling and integral analyses to the simulation results for various wall temperatures. The direct comparison between the present analyses and the simulation results is shown in figure 16(b). Because the maximum LST tends to occur after the superheat has been extinguished and the liquid has been consumed, its position and value are generally independent of the initial liquid temperature. The scaling analysis better predicts the trend in LST with wall temperature, as shown in figure 16(a). The actual values of the maximum LST for both analyses compare favourably to the numerical simulations, with the scaling analysis tending to overpredict and the integral analysis tending to underpredict the simulation results. Interestingly, the scaling analysis predicts the simulation values much better than the integral analysis for long solidification times when the wall temperature is near the solidus of the alloy. Other than this, both analyses are accurate to within approximately 18 % of the simulation results.
Figure 15. Local solidification time as a function of distance from the chill wall from a full numerical simulation of Al–4.5 wt% Cu in a 0.25 m square cavity, initially at 1000 K and cooled with an isothermal sidewall at 500 K. The vertical line shows the approximate thickness of the solid
$({\it\delta}_{s})$
at the time when the liquid has been completely consumed
$({\it\delta}=L)$
.
Figure 16. A comparison of the simulation results for maximum local solidification time (a) to the trends predicted by the scaling and integral analyses and (b) directly to the integral analysis.
Based on the above results, both the scaling analysis and the integral analysis are found to follow the trends predicted by the simulation with good agreement. Generally, the integral solution is found to be superior to the scaling analysis, although this is at the expense of increased complexity. While the integral analysis does not always predict the trends as closely as the scaling analysis, it is generally better at predicting the actual values, with the exception of maximum LST for long solidification times, and is within approximately 20 % of the simulation results for all of the quantities of interest.
5.3. Effect of cavity aspect ratio
In the numerical simulations,
$t_{m}(L)$
,
$t_{s}(L)$
, the cool-down time and the maximum LST were not strong functions of the aspect ratio. Additionally, the trends for each with the aspect ratio were not monotonic, as shown in figure 17. The nature of the scaling and integral analyses presented above is that the predicted trends for these values are monotonic with aspect ratio and therefore do not match the simulation results.
Figure 17. Simulation results for (a)
$t_{m}(L)$
, (b)
$t_{s2}(L)$
, (c) the cool-down time and (d) the maximum local solidification time as a function of aspect ratio.
Figure 18 shows isotherms near the liquidus temperature for various cavity heights. It is clear that the assumption of a flat solidification interface becomes less valid as the height of the cavity increases. This effect is caused by the increased thermal stratification over the height of the cavity. The collection of relatively cool fluid near the bottom of the cavity locally increases the rate of solidification. Additionally, as the height increases, the flow at the bottom of the cavity becomes more vigorous and the thermal field of the bulk liquid becomes more two-dimensional. These details of the complexity of the temperature field in the liquid as a function of aspect ratio are not adequately captured by the approximate analyses presented here, and therefore the trends predicted are not the same as those observed in the simulations.
Figure 18. Isotherms starting at the liquidus temperature and increasing to the right by 1 K for various aspect ratios, showing the change in the details of the temperature distribution in the liquid: (a)
$A=3$
, (b)
$A=2$
, (c)
$A=1$
, (d)
$A=1/2$
and (e)
$A=1/3$
. Plots are made at 100 s into the process. The length of the cavity in the
$x$
direction is 0.25 m for all cases.
While the trends with aspect ratio are not well predicted by the proposed models, the simulation results are also not strong functions of the aspect ratio, as can be seen by the range of data in figure 17. The simple models find no significant variation in
$t_{m}(L)$
or
$t_{s2}(L)~({<}4\,\%)$
over the range of aspect ratios, which is consistent with the simulation results. The agreement is not as good for the cool-down time, which varies by approximately 10 % over the various aspect ratios for the simulations, and approximately 30 % for the variable superheat scaling analysis and integral analysis, probably as a result of changes in the flow patterns. However, it is clear that this difference in the cool-down time does not significantly affect the solidification conditions. Essentially, while the aspect ratio changes the fluid flow in the cavity, these changes have little effect on the solidification behaviour, which is still primarily controlled by conduction in the
$x$
direction. Therefore, the approximate models developed above for
$A=1$
can be applied over the range of aspect ratios presented.
6. Comparison to Lin et al.
The flow analysis in this study was inspired by Lin et al. (Reference Lin, Armfield and Patterson2007), who investigated the same configuration without solidification. They assumed that the fluid exiting the boundary layer did not mix with the fluid at the initial temperature. The flow entered the boundary layer at the initial temperature and exited at the wall temperature. In considering two isothermal reservoirs divided at
$Y_{i}$
, this configuration results in the cool-down time being equal to the time at which
$Y_{i}$
equals zero, the point when the liquid is uniformly at the wall temperature. With this temperature profile in the liquid (figure 19
a), there is no driving force for the boundary layer below
$Y_{i}$
. Thus, Lin et al. assumed that the boundary layer flow does not travel down the full height of the cavity, but instead exits into the bulk liquid at the interface between the fluid reservoirs. A comparison of the assumed temperature profiles and flow configuration of that work and the present study is shown in figure 19.
Figure 19. Comparison of assumed temperature profiles in the liquid for (a) the study by Lin et al. (Reference Lin, Armfield and Patterson2007) and (b) this work.
In order to investigate the validity of the flow pattern assumed by Lin et al. (Reference Lin, Armfield and Patterson2007), a simulation was done for a test case with no solidification. The system parameters for this case are given in table 3 and are similar to case 1 performed by Lin et al. The Rayleigh number
$\mathit{Ra}_{H,w}$
is redefined using the wall temperature by
$$\begin{eqnarray}\mathit{Ra}_{H,w}=\frac{g{\it\beta}(T_{0}-T_{w})H^{3}}{{\it\alpha}{\it\nu}}.\end{eqnarray}$$
Figure 20. Streamlines at various times for the test case in table 3: (a)
$t=5$
s, (b)
$t=15$
s, (c)
$t=25$
s, (d)
$t=35$
s, (e)
$t=45$
s and (f)
$t=55$
s.
Table 3. System parameters used for comparison to work by Lin et al.
Lin et al. (Reference Lin, Armfield and Patterson2007) based their assumptions on the simulation of a test case similar to that described in table 3 by observing isotherms similar to those shown in figure 2, which show the stratification of the liquid developing over time. However, as shown by the streamlines (figure 20) for the test case in table 3, the flow does not exit the boundary layer at some intermediate height, but instead flows to the bottom of the cavity throughout the process.
The expected consequence of Lin et al.’s (Reference Lin, Armfield and Patterson2007) assumed flow configuration is that it exaggerates the rate at which the bulk fluid cools by always subtracting fluid from the bulk at the highest temperature in the system in exchange for the addition of fluid at the lowest temperature, rather than accounting for the range of fluid temperatures in the boundary layer and the developing stratified layers. The difference between Lin et al.’s model and the scaling analysis presented here can be more closely examined by rearranging (3.34) and (3.35) to develop a scale for the volumetric flow rate out of the boundary layer, equal to the product of the boundary layer velocity and thickness:
$$\begin{eqnarray}v_{0}{\it\delta}_{T}\sim \sqrt{{\it\alpha}}y^{3/4}(g{\it\beta}{\rm\Delta}T)^{1/4}.\end{eqnarray}$$
Using the assumptions from Lin et al.’s (Reference Lin, Armfield and Patterson2007) work, the temperature difference that drives the flow is held constant throughout the process, while
$y$
is equal to
$Y_{i}$
, which is initially set to
$H$
and decreases with time. However, if the liquid in the boundary layer is taken as the average of the bulk and wall temperatures, then flow occurs over the full length of the wall,
$y$
is equal to
$H$
, while the temperature difference that drives the flow decreases with time. It is clear that, in both cases, the result is a decrease in the flow rate from the boundary layer as the process continues, although it is not immediately clear if the trends are similar considering that the changes in
$Y_{i}$
and
${\rm\Delta}T$
with time are not known explicitly. However, both models do predict very similar boundary layer behaviours at short times, where
$Y_{i}$
is similar in value to
$H$
and the temperature difference is close to the initial condition. Lin et al.’s model only becomes less accurate as
$Y_{i}$
becomes small and the assumption of two distinct fluid reservoirs breaks down.
The two models were investigated more closely by considering the average temperature of the bulk liquid over time in comparison to simulation results. This comparison was performed for the test case in table 3 as shown in figure 21. As expected, Lin et al.’s (Reference Lin, Armfield and Patterson2007) model reasonably resembles the cooling behaviour of the bulk liquid at short times, which, because this period accounts for the most significant cooling rate, explains the reported agreement in the trends between this approximate model and simulation results (Lin et al. Reference Lin, Armfield and Patterson2007). However, this model poorly predicts the cooling behaviour at long times, while the analysis in the present work is accurate throughout the process.
Figure 21. Comparison of average bulk fluid temperature over time for the case in table 3 for simulation results, the present scaling analysis and the analyses performed by Lin et al. (Reference Lin, Armfield and Patterson2007) and Amberg (Reference Amberg1997).
Amberg (Reference Amberg1997) also developed a simple model for the cooling behaviour of the cavity (without solidification) using a lumped capacitance approach for the bulk liquid and a scaling analysis for the convective heat flux from the boundary layer. The boundary layer is assumed to extend the full height of the cavity throughout the process, and the resulting heat flux is a function of the instantaneous fluid temperature. The equation for the average liquid temperature is
$$\begin{eqnarray}{\it\theta}=\frac{1}{(1+{\it\tau}/4)^{4}},\end{eqnarray}$$
where the following non-dimensional variables are used:
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\theta}=\frac{T_{avg}-T_{w}}{T_{0}-T_{w}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}=\left(\frac{t{\it\alpha}}{HL}\right)Bo_{H}^{1/4}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \mathit{Bo}_{H}=\mathit{Pr}\mathit{Ra}_{H,w}. & \displaystyle\end{eqnarray}$$
Amberg’s (Reference Amberg1997) solution for the test case in table 3 is also shown on the plot in figure 21, and agrees quite well with the simulation results. However, Amberg did not expand his approach to the case of solidification in which the volume of liquid is changing in time.
7. Discussion of assumptions and applicability of analysis
Several assumptions were made during the scaling of the flow field that limit the applicability of the present analysis. The first is the decision to neglect the effects of solutal buoyancy, which necessarily develop during the solidification of binary or multicomponent alloys. The relative contributions of solutal and thermally driven buoyancy may be expressed by the buoyancy ratio,
$N$
:
$$\begin{eqnarray}N=\frac{{\it\beta}_{C}(C-C_{0})}{{\it\beta}_{T}(T-T_{avg})}.\end{eqnarray}$$
For
$|N|\gg 1$
, composition effects are dominant, whereas thermal effects predominate for
$|N|\ll 1$
, and a mixed condition exists for
$|N|\sim 1$
. The sign of
$N$
describes the relative direction of the two forces, where positive values indicate aiding effects, and negative values opposing. For negative values of
$N$
, counterflow may develop with two boundary layers flowing in opposite directions (Nilson Reference Nilson1985). Counterflow is a common concern in magmatic solidification (Jarvis & Huppert Reference Jarvis and Huppert1995), and is present during processing of some metallic systems, such as lead-rich Pb–Sn alloys (Krane & Incropera Reference Krane and Incropera1997), where the liquid near the solidification front is enriched in tin, causing an inner solutally driven upward-flowing boundary layer bordering an outer thermally driven downward flow. Evaluating
$N$
for metallic systems, such as the one considered here, can be difficult because the instantaneous value of the superheat may change dramatically throughout the process, commonly starting at values of the order of 100 K, and necessarily vanishing as solidification is completed. The magnitude of the composition differences may be as large as a few weight per cent, but these are almost entirely confined within the mushy zone, in which the flow is greatly restricted by the limited permeability, as will be discussed shortly.
In the case of the Al–Cu system discussed here, where the solutal expansion coefficient is
$-0.73$
per unit mass fraction Cu, the sign of
$N$
is positive and if the compositional difference is estimated as 1 wt% combined with a conservative estimate of the superheat at 10 K,
$N\sim O(10)$
, suggesting that solutal effects are dominant near the edge of the mushy zone. However, the influence of solutal buoyancy on the bulk fluid flow is significantly limited by the ratio of thermal to mass diffusivities, which is described with the Lewis number, Le:
$$\begin{eqnarray}\mathit{Le}=\frac{\mathit{Sc}}{\mathit{Pr}}=\frac{{\it\alpha}}{D},\quad \text{where }\mathit{Sc}=\frac{{\it\nu}}{D}.\end{eqnarray}$$
Generally, for metallic alloys,
$\mathit{Le}\sim O(10^{3})$
or greater, and for this Al–Cu system,
$\mathit{Le}=10\,500$
. This high Lewis number suggests that the flow will only be strongly affected by solutal effects within the mushy zone, as there will be very little solutal diffusion into the bulk liquid, and so the bulk fluid flow will primarily be a function of the thermal effects, and the analysis presented here is valid. For cases with opposing buoyancy forces, the present analysis may still apply as long as the flow in the mushy zone is much smaller than
$v_{0}$
.
The analysis further assumes that the mush is impermeable, which may be considered valid if the fluid flow in this region, driven primarily by solutal effects, is slow in comparison to the adjacent thermally driven boundary layer flow in the bulk. This assumption may be evaluated by considering the vertical momentum equation within the mush:
$$\begin{eqnarray}u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}={\it\nu}\left(\frac{\partial ^{2}v}{\partial x^{2}}+\frac{\partial ^{2}v}{\partial y^{2}}\right)-g({\it\beta}_{T}{\rm\Delta}T+{\it\beta}_{C}{\rm\Delta}C)-\frac{{\it\nu}v}{K}.\end{eqnarray}$$
The appropriate length scale,
${\it\lambda}$
, for the permeability (
$K$
) in metallic alloys is the order of the dendrite arm spacing, and is approximately
$O(10^{-5}~\text{m})$
. If shrinkage-induced flow is neglected, then the horizontal velocity component may be taken as zero, since no buoyancy forces are acting in this direction. Subsequently, by continuity, the vertical velocity must be approximately constant. Under these conditions and assuming solutally driven flow (
$N>1$
), the vertical moment equation may be reduced and differentials replaced by reference differences:
$$\begin{eqnarray}\begin{array}{@{}lcc@{}}\displaystyle \frac{{\it\nu}v}{K}, & \displaystyle {\it\nu}\frac{v_{ref}}{{\rm\Delta}x_{ref}^{2}}\sim & g{\it\beta}_{C}{\rm\Delta}C_{ref}.\\ \,P & F\,\, & B\end{array}\end{eqnarray}$$
The driving buoyancy force (
$B$
) is always retained and is balanced by a restraining force, either permeability (
$P$
) or friction (
$F$
). The ratio of these two terms may be used to determine their relative importance:
$$\begin{eqnarray}\frac{P}{F}\sim \frac{{\rm\Delta}x_{ref}^{2}}{K}.\end{eqnarray}$$
For liquid metals, the appropriate length scale for the numerator is the width of the mush, which is approximately
$O(10^{-2}~\text{m})$
for the cases described here. Evaluating the permeability near the edge of the mushy zone by using a relatively small value for the volume fraction solid (e.g. 0.1) and an approximate dendrite arm spacing
${\it\lambda}\sim O(10^{-5}~\text{m})$
give that
$K\sim O(10^{-10}~\text{m}^{2})$
. The ratio in (7.5), then, is
$O(10^{8})$
, which shows that the viscous forces are small relative to the permeability forces in the mushy zone, even near the liquidus interface. The viscous term may then be dropped from (7.4) and a velocity scale found for the edge of the mushy zone:
$$\begin{eqnarray}v_{ref}\sim \frac{K}{{\it\nu}}(g{\it\beta}_{C}{\rm\Delta}C_{ref}).\end{eqnarray}$$
If the compositional difference is approximately 1 wt%, and considering that
${\it\beta}_{C}\sim O(1)$
, the resulting order of magnitude for the velocity is
$v_{ref}\sim O(10^{-5}~\text{m}~\text{s}^{-1})$
. If the same scales are used to calculate the boundary layer velocity outside of the mushy zone in (3.34), it is found that
$v_{0}\sim O(10^{-1.5}~\text{m}~\text{s}^{-1})$
, which is several orders of magnitude larger. Based on these results, it can be safely concluded that the flow in the outer region of the mushy zone and its direction has only a very small influence on the flow in the liquid region boundary layer.
The validity of these assumptions was tested by running the base case (table 2) with and without compositional effects. The positions of the liquidus and solidus interface in time are shown for both cases in figure 22. The position of the liquidus interface deviates between the two cases for long times. This result may be attributed to two causes. The first is that the superheat necessarily approaches zero over time, slowly reducing the magnitude of the thermal buoyancy. As the superheat becomes small, the relative importance of the solutal effects increases, associated with an increase in
$N$
. The second effect is that enrichment of the bulk liquid over time also reduces the liquidus temperature, as shown in the phase diagram in figure 1, which will tend to delay solidification. The differences between the solidus positions are small at all times.
Figure 22. Comparison of progression of solidus and liquidus fronts with and without solutal effects for the base case in table 2.
8. Conclusions
The unidirectional solidification of a binary metallic alloy in a rectangular cavity with an isothermal sidewall boundary condition has been investigated using scaling and integral analyses. This study differed from previous solutions by coupling the solidification behaviour to the bulk fluid flow, including the extinguishing of the superheat during the process. The flow was approximated using a boundary layer solution for small Prandtl numbers characteristic of liquid metals. These models were compared to numerical simulations for an Al–4.5 wt% Cu alloy with various boundary and initial conditions and cavity aspect ratios. The scaling analysis generally achieved good agreement in the trends observed in the numerical model. The integral analysis showed similar agreement with the simulation trends, and also correctly predicted several important results, including the maximum LST of the casting. The inclusion of the effect of fluid flow on solidification also resulted in an extension of previous work on cooling of a cavity in the absence of solidification, particularly the development of a more accurate depiction of the transient thermal stratification in the bulk fluid.
The results presented here were based on important assumptions that neglected compositional buoyancy effects and the permeability of the mushy zone. Neglecting solutal effects outside of the mushy zone was justified by the very high Lewis number. The strength of the solutally driven flow that might occur in the mushy zone was shown to be negligible compared to the thermally driven boundary layer in the bulk liquid.
The isothermal boundary condition considered in this study is convenient from an analytical perspective, but is often not a good approximation of actual casting situations. An improvement upon this work would be to consider the general case of a uniform heat transfer coefficient boundary condition. Complications in this analysis include the necessary consideration of time regimes containing only liquid or a combination or liquid and mush in cases in which any initial superheat is present. Additionally, inclusion of the rejected solute at the solidification front and its effect on the fluid flow and solidification behaviour would improve the predictive power of these models, particularly in cases where counterflow develops a separate opposing flow cell driven by solutal buoyancy effects. The limited mass diffusion in liquid metals suggests that consideration of the solutal effect may also require accounting for the permeability of the mush near the liquidus interface. This type of analysis would allow for predictions of the macro-segregation in the final ingot.
