1. Introduction
The one-dimensional problem for the solidification of a cooling binary alloy melt which simultaneously undergoes solidification shrinkage is the simplest configuration in which to examine the role of convection in a binary system; as such, it can be used as a benchmark for testing higher-dimensional numerical codes. Furthermore, the fact that there is relative motion between solid and liquid phases in the so-called mushy zone, where the two phases coexist, leads to the possibility of severe macrosegregation [Reference Chen and Tsai10, Reference Diao and Tsai13, Reference Diao and Tsai14, Reference Minakawa, Samarasekera and Weinberg25, Reference Mo26], i.e. variations in composition that occur in alloy castings or ingots and range in scale from several millimetres to centimetres or even metres; this is a central problem in industry, since it strongly influences the further workability of the cast products and their mechanical properties.
The mathematical problem to be solved was recently formulated in [Reference Assunção, Vynnycky and Mitchell6]; however, only similarity solutions available for short times were computed, although it was suggested that these solutions could be used as initial conditions for solving the problem for all time, when the problem no longer has a similarity solution. Thus, the purpose of this paper is to tackle the remainder of the problem, and thereby solve it up to complete solidification; doing so would then yield the macrosegregation profile over the entire domain. The novelty of the work is that, in so doing, we are able to shed light on earlier experimental and theoretical results that have been obtained for close variants of this problem [Reference Chen and Tsai10, Reference Kato and Cahoon21, Reference Kirkaldy and Youdelis22, Reference Minakawa, Samarasekera and Weinberg25, Reference Rousset, Rappaz and Hannart29, Reference Voller31, Reference Voller and Sundarraj34], in particular as regards the possibility of severe macrosegregation as complete solidification is approached. Also, the solution to this problem, which is time-dependent in nature, turns out to be of some relevance to the phenomenon of centreline segregation in steady-state continuous casting processes. We will carry out our study using a moving mesh formulation, rather than the more commonly used enthalpy method, which employs a fixed mesh. As a consequence, this avoids the need to use arbitrary artificial parameters to ensure numerical convergence [Reference Amberg2, Reference Assunção and Vynnycky5]. Moreover, by tracking the solidus and liquid interfaces explicitly, it is hoped to be able to avoid the usual issues of numerical dispersion and diffusion that are present in the numerically simulated macrosegregation profiles reported in the literature [Reference Du, Eskin and Katgerman15]. In fact, to the best of our knowledge, a moving mesh formulation for this type of problem has been used only once [Reference Rousset, Rappaz and Hannart29]; however, the mathematical details given were comparatively sparse and there is considerable doubt about the accuracy of the results exactly in the region of complete solidification.
The layout of the paper is as follows. In Section 2, we re-cap the nondimensionalized governing equations from [Reference Assunção, Vynnycky and Mitchell6]. In Section 3, we discuss the numerical implementation of these equations; the results are presented in Section 4. In Section 5, the results are discussed in the context of earlier attempts to solve similar problems and their relevance to the modelling of macrosegregation in continuous casting processes; conclusions are then drawn in Section 6.
2. Mathematical formulation
We consider a transient one-dimensional problem, as shown in Figure 1, in which a region of binary alloy melt of initial extent
$l,$
temperature
$T_{\text{cast}}$
and composition
$C_{0}$
, is cooled by a boundary at
$y=0$
that is held at temperature
$T_{\mathrm{w}}.$
To fix ideas more precisely, we suppose that the alloy in question has a phase diagram with a eutectic point, with eutectic temperature and concentration
$T_{\mathrm{eut}}$
and
$C_{\mathrm{eut}}$
respectively, as indeed is the case for many alloys of technological interest. Figure 2 shows the phase diagram for the Al-Cu system, which was the one implicitly considered by [Reference Voller31] and the one that we will consider here, and we will also assume that
$T_{\mathrm{w}}\leq T_{\mathrm{eut}},$
which ensures that complete solidification occurs for any choice of
$C_{0}$
, such that
$0\leq C_{0}\leq C_{\mathrm{eut}}.$
As indicated in Figure 1, at some time
$t\gt 0,$
the region that initially occupied
$0\leq y\leq l$
will shrink to occupy
$0\leq y\leq y_{\infty } ( t ),\,$
which consists of a solid region occupying
$0\leq y\leq y_{\mathrm{s}} ( t ),$
a mushy zone occupying
$y_{\mathrm{s}} ( t ) \leq y\leq y_{\mathrm{l}} ( t )$
and the melt region occupying
$y_{\mathrm{l}} ( t ) \leq y\leq y_{\mathrm{\infty }} ( t ) ;$
the shrinkage occurs as a result of phase change, as molten melt begins to solidify, with the density of the solid being greater than that of the molten phase. In addition, it is assumed that the solid is stationary, corresponding to the case of columnar or consolidated equiaxed dendritic solidification.
Solidification with shrinkage: (a) initial configuration at.
$t=0$
; (b) after time
$t\gt 0$
.
Al-rich side of the linearised phase diagram for the Al–Cu system.
$L$
denotes liquid phase, and
$\alpha$
and
$\theta$
denote two solid phases. The eutectic point in this figure,
$(C_{\mathrm{eut}},T_{\mathrm{eut}})$
, is at (33.2 Wt% Cu, 821 K).
We re-cap the model equations briefly here; more complete details, including a discussion of how the less obvious ones are motivated, are given in [Reference Assunção, Vynnycky and Mitchell6]. We have
-
in the solid, the conservation of energy is given by
(2.1)where
\begin{align} \rho _{\mathrm{s}}c_{\mathrm{ps}}\frac{\partial T}{\partial t}=K_{\mathrm{s}}\frac{\partial ^{2}T}{\partial y^{2}}\text{,} \end{align}
$T$
denotes the temperature, and
$\rho _{\mathrm{s}},c_{\mathrm{ps}}$
and
$K_{\mathrm{s}}$
the density, specific heat capacity at constant pressure and thermal conductivity, respectively, of the solid phase, all of which are assumed to be constant;
-
in the mush, the total mass balance over solid and liquid phases in the mush, derived from summing the phase conservation equations as in [Reference Bennon and Incropera8], is given by
(2.2)with
\begin{align} \left ( 1-R\right ) \frac{\partial \chi }{\partial t}+\frac{\partial U}{\partial y}=0, \end{align}
$R=\rho _{\mathrm{s}}/\rho _{\mathrm{l}}$
and
$U=\chi v,$
where
$\rho _{\mathrm{l}}$
is the density of the liquid phase,
$\chi$
is the liquid fraction and
$v$
is the liquid phase velocity. The conservation of energy is given by(2.3)where
\begin{align}& c_{\mathrm {ps}}\frac {\partial }{\partial t}\left ( \rho _{\mathrm {s}}\!\left ( 1-\chi \right ) \left ( T-T_{\mathrm {ref}}\right ) \right ) +c_{\mathrm {pl}}\!\left \{ {\frac {\partial }{\partial t}\left ( \rho _{\mathrm {l}}\chi \!\left ( T-T_{\mathrm {ref}}\right ) \right ) +\frac {\partial }{\partial y}\left ( \rho _{\mathrm {l}}U\!\left ( T-T_{\mathrm {ref}}\right ) \right ) }\right \}\nonumber\\& \quad =\frac{\partial }{\partial y}\left ( K\frac{\partial T}{\partial y}\right ) -\Delta H_{\mathrm{f}}\frac{\partial }{\partial t}\left ( \rho _{\mathrm{l}}\chi \right ), \end{align}
$c_{\mathrm{pl}}$
and
$K_{\mathrm{l}}$
denote the specific heat capacity at constant pressure and thermal conductivity, respectively, of the liquid phase,
$\Delta H_{\mathrm{f}}$
is the latent heat of fusion and
$K$
is the mixture thermal conductivity, given by(2.4)Also,
\begin{align} K{=\chi K_{\mathrm{l}}+\left ( 1-\chi \right ) K_{\mathrm{s}}.} \end{align}
$T_{\mathrm{ref}}$
is a reference temperature, to which thermal energies in all phases must refer consistently; a convenient choice for
$T_{\mathrm{ref}}$
is
$T_{\mathrm{m}},$
the melting point of the solvent element. The equation for the conservation of solute, taken over solid and liquid phases, is given by(2.5)where the form for the mixture concentration,
\begin{align}{\frac{\partial }{\partial t}\left ( \rho C\right ) +\frac{\partial }{\partial y}\left ( \rho _{\mathrm{l}}C_{\mathrm{l}}U\right ) =0,}\quad \end{align}
$C$
, depends on the assumption made regarding solute transport at the microscale. To take into account all possibilities, we write(2.6)with
\begin{align}{\rho C=\chi \rho _{\mathrm{l}}C_{\mathrm{l}}+\rho _{\mathrm{s}}\int _{0}^{1-\chi }\left ( C_{\mathrm{s}}+\beta \chi ^{\prime }\frac{\mathrm{d}C_{\mathrm{s}}}{\mathrm{d}\chi ^{\prime }}\right ) \mathrm{d}\chi ^{\prime },} \end{align}
$C_{\mathrm{l}}$
and
$C_{\mathrm{s}}$
as the concentrations of the solute in the liquid and solid phases, respectively, which are related by(2.7)where
\begin{align}{C_{\mathrm{s}}=k_{0}C_{\mathrm{l}},} \end{align}
$k_{0}$
is the partition coefficient. The parameter
$\beta,$
where
$0\leq \beta \leq 1,$
determines the extent of solute back diffusion at the microscale; thus,
$\beta =1$
corresponds to the lever rule, i.e. complete back diffusion, whereas
$\beta =0$
corresponds to the Scheil equation, i.e. no back diffusion. Consequently, equation (2.5) becomes(2.8)this is demonstrated in Appendix A. We also have local thermodynamic equilibrium in the mushy region, so that, from the liquidus line in the phase diagram,
\begin{align} \rho _{\mathrm{l}}\frac{\partial }{\partial t}\left ( \chi C_{\mathrm{l}}\right ) -k_{0}\rho _{\mathrm{s}}C_{\mathrm{l}}\frac{\partial \chi }{\partial t}+\rho _{\mathrm{s}}\beta \!\left ( 1-\chi \right ) k_{0}\frac{\partial C_{\mathrm{l}}}{\partial t}+\frac{\partial }{\partial y}\left ( \rho _{\mathrm{l}}C_{\mathrm{l}}U\right ) =0; \end{align}
(2.9)where
\begin{align} T=T_{\mathrm{m}}-mC_{\mathrm{l}}\quad \text{for }0\leq \chi \leq 1, \end{align}
$m\gt 0;$
note here that, in view of the minus sign in equation (2.9). Note also that we have not mentioned any equation for the conservation of momentum; this was discussed in appendix A in [Reference Assunção, Vynnycky and Mitchell6] and reflects the fact that, in this one-dimensional flow, the velocity can be determined ahead of the pressure.
-
in the melt, here, the mass balance can be written as just
(2.10)which can be thought of as the result of setting
\begin{align} \frac{\partial U}{\partial y}=0, \end{align}
$\chi =1$
in equation (2.2). The conservation of energy is given by(2.11)also derived considering
\begin{align}{\rho _{\mathrm{l}}c_{\mathrm{pl}}\left ( \frac{\partial }{\partial t}\left ( T-T_{\mathrm{ref}}\right ) +\frac{\partial }{\partial y}\left ( U\!\left ( T-T_{\mathrm{ref}}\right ) \right ) \right ) =K_{\mathrm{l}}\frac{\partial ^{2}T}{\partial y^{2}}\text{,}}\quad \end{align}
$\chi =1$
in equation (2.3), while neglecting changes in material properties inside the melt region, whereas the conservation of solute, as in [Reference Flemings17], is given by(2.12)
\begin{align} \frac{\partial }{\partial t}\left ( \rho _{\mathrm{l}}C_{\mathrm{l}}\right ) +\frac{\partial }{\partial y}\left ( \rho _{\mathrm{l}}C_{\mathrm{l}}U\right ) =0. \end{align}
For simplicity, we henceforth set
$c_{\mathrm{ps}}=c_{\mathrm{pl}}$
and
$K_{\mathrm{s}}=K_{\mathrm{l}}.$
This is a convenient reduction, since all four quantities vary with temperature, as well as alloy composition: for example, for Al-3.4 wt% Cu, thermal conductivity decreases by around 30% and specific heat capacity increases by around 60% as the temperature increases from 300 K to 1073 K [Reference Rousset, Rappaz and Hannart29]. Note also that this will imply that the mixture thermal conductivity is independent of
$\chi ;$
however, this is not expected to impact the qualitative features of the resulting macrosegregation profile at leading order and is indeed a common enough assumption in models of this type [Reference Swaminathan and Voller30–Reference Voller, Mouchmov and Cross33]. More importantly from the point of view of macrosegregation, however, we keep
$\rho _{\mathrm{s}}\neq \rho _{\mathrm{l}}$
throughout.
2.1. Boundary and interface conditions
2.1.1. Mass
The boundary conditions for
$U$
are, at
$y=y_{\mathrm{s}} ( t ),$
\begin{align} U=\left ( 1-R\right ) \chi _{\mathrm{s}}\dot{y}_{\mathrm{s}}\left ( t\right ) \quad \end{align}
and, at
$y=y_{\mathrm{l}} ( t )$
\begin{align} U=\dot{y}_{\infty }\left ( t\right ) .\quad \end{align}
2.1.2. Heat and solute
At
$y=0,$
we have
\begin{align} T=T_{\text{w}}. \end{align}
On the other hand, at the upper surface of the melt, which would typically be in contact with air [Reference Rousset, Rappaz and Hannart29], it is reasonable to expect the continuity of heat flux, which will consist of both diffusive and convective components. Since the temperature and the product of normal component of velocity and material density will both be continuous across this interface, and since the specific heat capacity of air is similar to that of the melt, we will have that the convective flux will be continuous to a good approximation; hence, the diffusive flux is continuous. However, since the thermal conductivity of air is many orders of magnitude of smaller than that of the melt, we will obtain just
\begin{align} \frac{\partial T}{\partial y}=0\text{ at }y=y_{\infty }\!\left ( t\right ) . \end{align}
Note that there is actually no need for a boundary condition on
$C_{\mathrm{l}}$
at
$y=y_{\infty } ( t ),$
as it is governed by a first-order hyperbolic partial differential equation (PDEs). However, we show in Section 4 that it turns out that
\begin{align}{\frac{\partial C_{\mathrm{l}}}{\partial y}=0} \end{align}
when there is melt at
$y=y_{\infty } ( t ) ;$
on the other hand, when there is mush, then we also arrive at (2.17), via (2.9) and (2.16).
At
$y=y_{\mathrm{s}} ( t ),$
we have
\begin{align} \left [ T\right ] _{-}^{+} & =0\text{ at }y=y_{\mathrm{s}}\!\left ( t\right ), \end{align}
\begin{align}{\left [ K\frac{\partial T}{\partial y}\right ] _{-}^{+}} & ={\displaystyle -\rho _{\mathrm{l}}\Delta H_{\mathrm{f}}\chi \frac{\mathrm{d}y_{\mathrm{s}}}{\mathrm{d}t}\text{ at }y=y_{\mathrm{s}}\!\left ( t\right ),}\quad \end{align}
with
$ [ \text{ } ] _{-}^{+}$
denoting the difference in the value of that function above and below
$y=y_{\mathrm{s}} ( t )$
, noting that the value of
$K$
on the minus side in (2.19) is implicitly taken to be
$K_{\mathrm{s}}$
, and
\begin{align}{\left \{ \begin{array} [c]{l@{\quad}l}\chi =0,\quad T=T_{\mathrm{m}}-mC_{\mathrm{l}}, & \text{if }C_{\mathrm{l}}\lt C_{\mathrm{eut}}\\[5pt] T=T_{\mathrm{eut}}, & \text{if }C_{\mathrm{l}}=C_{\mathrm{eut}}\end{array} \right . \quad \text{at }y=y_{\mathrm{s}}\!\left ( t\right ) .} \end{align}
Furthermore, we shall henceforth write
$\chi _{\text{s}}= \chi \vert _{y=y_{\mathrm{s}} ( t ) }.$
At
$y=y_{\mathrm{l}} ( t ),$
we have
\begin{align} \left [ T\right ] _{-}^{+} & =0, \end{align}
\begin{align}{\left [ K\frac{\partial T}{\partial y}\right ] _{-}^{+}} =0, \end{align}
where the value of
$K$
on the plus side in (2.22) is implicitly taken to be
$K_{\mathrm{l}}$
, and, analogous to the first alternative in (2.20),
\begin{align}{\chi =1,\quad T=T_{\mathrm{m}}-mC_{\mathrm{l}}.}\quad \end{align}
2.2. Initial conditions
For
$y\gt 0,$
\begin{align} T & =T_{\text{cast}}, \end{align}
\begin{align} C_{\mathrm{l}} & =C_{0}, \end{align}
as well as
\begin{align} y_{\mathrm{l}}\!\left ( 0\right ) & =0, \end{align}
\begin{align} y_{\mathrm{s}}\!\left ( 0\right ) & =0, \end{align}
\begin{align} y_{\infty }\!\left ( 0\right ) & =l. \end{align}
Furthermore, these yield a further condition that describes the conservation of solute over the whole system:
\begin{align} \int _{0}^{y_{\infty }\left ( t\right ) }\rho C\mathrm{d}y=\rho _{\mathrm{l}}C_{0}l. \end{align}
2.3. Nondimensionalized governing equations
In order to nondimensionalize the equations, we set
\begin{equation*} Y=\frac {y}{l},\quad \tau =\frac {t}{\left [ t\right ] },\quad \theta =\frac {T-T_{\text {w}}}{T_{\text {cast}}-T_{\text {w}}},\quad \hat {U}=\frac {U}{l/\left [ t\right ] }, \end{equation*}
\begin{align} \varrho =\frac{\rho }{\rho _{\mathrm{l}}},\quad \mathcal{C}_{\mathrm{l}}=\frac{C_{\text{l}}}{C_{0}},\quad \mathcal{C}_{\mathrm{s}}=\frac{C_{\text{s}}}{C_{0}},\quad \quad Y_{\mathrm{s}}=\frac{y_{\mathrm{s}}}{l},\quad Y_{\mathrm{l}}=\frac{y_{\mathrm{l}}}{l},\quad Y_{\infty }=\frac{y_{\infty }}{l}, \end{align}
where
$ [ t ]$
is a time scale, which we shall take as being given by
$ [ t ] =\rho _{\mathrm{l}}c_{\mathrm{pl}}l^{2}/K_{\mathrm{l}}.$
In the solid, where
$0\leq Y\leq Y_{\mathrm{s}} ( \tau )$
,
\begin{align} R\frac{\partial \theta }{\partial \tau }=\frac{\partial ^{2}\theta }{\partial Y^{2}}. \end{align}
In the mush, where
$Y_{\mathrm{s}} ( \tau ) \leq Y\leq Y_{\mathrm{l}} ( \tau ),$
\begin{align} \frac{\partial }{\partial \tau }\left ( \chi +\left ( 1-\chi \right ) R\right ) +\frac{\partial \hat{U}}{\partial Y}=0, \end{align}
\begin{align} \left [ \chi +\left ( 1-\chi \right ) R\right ] \frac{\partial \theta }{\partial \tau }+\hat{U}\frac{\partial \theta }{\partial Y}=\frac{\partial ^{2}\theta }{\partial Y^{2}}-\frac{R}{St}\frac{\partial \chi }{\partial \tau }, \end{align}
\begin{align} \frac{\partial }{\partial \tau }\left ( \chi \mathcal{C}_{\mathrm{l}}\right ) -k_{0}R\mathcal{C}_{\mathrm{l}}\frac{\partial \chi }{\partial \tau }+\beta R\!\left ( 1-\chi \right ) k_{0}\frac{\partial \mathcal{C}_{\mathrm{l}}}{\partial \tau }+\frac{\partial }{\partial Y}\left ( \mathcal{C}_{\mathrm{l}}\hat{U}\right ) =0, \end{align}
\begin{align} \theta =\theta _{\mathrm{m}}-\tilde{m}\mathcal{C}_{\mathrm{l}}, \end{align}
where
$St=c_{\mathrm{pl}} ( T_{\text{cast}}-T_{\text{w}} )/\Delta H_{\mathrm{f}}$
and
\begin{equation*} \theta _{\mathrm {m}}=\frac {T_{\mathrm {m}}-T_{\text {w}}}{T_{\text {cast}}-T_{\text {w}}},\text { }\tilde {m}=\frac {mC_{0}}{T_{\text {cast}}-T_{\text {w}}}. \end{equation*}
In the melt, where
$Y_{\mathrm{l}} ( \tau ) \leq Y\leq Y_{\infty } ( \tau ),$
\begin{align} \frac{\partial \hat{U}}{\partial Y} & =0, \end{align}
\begin{align} \frac{\partial \theta }{\partial \tau }+\hat{U}\frac{\partial \theta }{\partial Y} & =\frac{\partial ^{2}\theta }{\partial Y^{2}}, \end{align}
\begin{align} \frac{\partial \mathcal{C}_{\mathrm{l}}}{\partial \tau }+\frac{\partial }{\partial Y}\left ( \mathcal{C}_{\mathrm{l}}\hat{U}\right ) & =0. \end{align}
Note, incidentally, that there is no trace of any terms involving the reference temperature in equations (2.33) and (2.37), as these vanish as a consequence of (2.32) and (2.36) and the fact that we have taken
$c_{\mathrm{ps}}=c_{\mathrm{pl}}$
.
The boundary and interfacial conditions are then at
$Y=0$
,
\begin{align} \theta =0; \end{align}
at
$Y=Y_{\mathrm{s}}(\tau ),$
\begin{align} \left \{ \begin{array}{l@{\quad}l}\chi =0,\quad \theta =\theta _{\mathrm{m}}-\tilde{m}\mathcal{C}_{\mathrm{l}} & \text{if }\mathcal{C}_{\mathrm{l}}\lt \mathcal{C}_{\mathrm{eut}}\\[5pt] \theta =\theta _{\mathrm{eut}} & \text{if }\mathcal{C}_{\mathrm{l}}=\mathcal{C}_{\mathrm{eut}}\end{array} \right ., \end{align}
\begin{align} \left [ \frac{\partial \theta }{\partial Y}\right ] _{-}^{+} & =-\frac{\chi _{\text{s}}}{St}\dot{Y}_{\mathrm{s}}(\tau ), \end{align}
\begin{align} \hat{U} & =\left ( 1-R\right ) \chi _{\text{s}}\dot{Y}_{\mathrm{s}}(\tau ), \end{align}
where
\begin{equation*} \theta _{\mathrm {eut}}=\frac {T_{\mathrm {eut}}-T_{\text {w}}}{T_{\text {cast}}-T_{\text {w}}},\quad \mathcal {C}_{\mathrm {eut}}=\frac {C_{\mathrm {eut}}}{C_{0}}, \end{equation*}
with
$\theta _{\mathrm{eut}}=\theta _{\mathrm{m}}-\tilde{m}\mathcal{C}_{\mathrm{eut}};$
at
$Y=Y_{\mathrm{l}}(\tau ),$
\begin{align} \theta & =\theta _{\mathrm{m}}-\tilde{m}\mathcal{C}_{\mathrm{l}}, \end{align}
\begin{align} \left [ \frac{\partial \theta }{\partial Y}\right ] _{-}^{+} & =0, \end{align}
\begin{align} \chi & =1; \end{align}
at
$Y=Y_{\infty } ( \tau ),$
\begin{align} \frac{\partial \theta }{\partial Y}=0. \end{align}
As for the initial conditions, we have, for
$0\lt Y\leq 1,$
\begin{align} \theta & =1, \end{align}
\begin{align} \mathcal{C}_{\mathrm{l}} & =1, \end{align}
as well as
\begin{align} Y_{\mathrm{l}}\left ( 0\right ) & =0, \end{align}
\begin{align} Y_{\mathrm{s}}\!\left ( 0\right ) & =0, \end{align}
\begin{align} Y_{\infty }\!\left ( 0\right ) & =1. \end{align}
Also, we have
\begin{align} \int _{0}^{Y_{\infty }\left ( \tau \right ) }\varrho \mathcal{C}\mathrm{d}Y=1.\quad \end{align}
2.4. Dimensionless parameters
Before proceeding further, it is worth noting that there are now seven dimensionless parameters left:
$k_{0},\tilde{m},R,St,\beta,\theta _{\mathrm{m}}$
and
$\theta _{\mathrm{eut}}.$
Note also that none of the remaining dimensionless parameters depend on the initial domain size,
$l$
, meaning that the results will not depend on it either. Moreover, using the data in Table 1 for the same binary alloy as in [Reference Assunção, Vynnycky and Mitchell6], i.e. Al-5 wt% Cu, we have
\begin{equation*} \tilde {m}\approx 0.05,\quad R\approx 1.06,\quad St\approx 0.91,\quad \theta _{\mathrm {m}}\approx 0.97,\quad \theta _{\mathrm {eut}}\approx 0.63.\quad \end{equation*}
Since most of these parameters, as well as
$k_{0}$
and
$\beta,$
are all of
$O ( 1 ),$
there appear to be no meaningful simplifications to be made; for example,
$\tilde{m}=0$
is merely the case of pure solvent. We can note, however, that
$R$
is typically no greater than around 1.1 even for other alloys, meaning that one might think to treat
$(R-1)$
as a small parameter and consider a regular perturbation expansion in this parameter. This was indeed done in [Reference Vynnycky and Saleem40] for a problem similar to that described in Section 5.2 for the particular case of
$\beta =1$
, albeit with rather minimal advantage, as the algebra soon became rather lengthy; in fact, although it was not pursued, it would have become even lengthier still for
$\beta \lt 1$
. Observe also that the equations do become much simpler when
$R=1;$
however, there will be no macrosegregation in this case. Moreover, once a particular binary alloy with a particular initial composition has been chosen, then
$T_{\mathrm{m}},T_{\mathrm{eut}},m,C_{0}$
and
$k_{0}$
are prescribed; in addition, if
$T_{\text{w}}$
and
$T_{\text{cast}}$
are also given, then
$\theta _{\mathrm{m}},\theta _{\mathrm{eut}},\tilde{m},k_{0}$
and
$St$
will all be known, leaving just
$\beta$
and
$R$
as parameters that can be varied. Thus, in what follows, we will focus primarily on varying these two parameters.
Parameters for computations for Al-5 wt% Cu
3. Numerical considerations and implementation
A principal difficulty with solving these equations is that there is initially only molten phase; however, in view of the isothermal boundary condition (2.39), the fully solid and mush regions begin to grow instantaneously. To account for this, Assunção et al. [Reference Assunção, Vynnycky and Mitchell6] introduced a triple Landau transformation, and considered the resulting PDEs in the limit as
$\tau \rightarrow 0.$
This gave ordinary differential equations (ODEs), the solution of which yielded initial conditions for the Landau-transformed system of PDEs. A feature of this system was the presence of a boundary layer of width
$\tau ^{1/2}$
in temperature on the liquid side of the mush-liquid region for
$R\neq 1$
in the Landau coordinates
$,$
which is precisely the case of interest for macrosegregation; in addition, note that since it was also found that
$Y_{\mathrm{s}},Y_{\mathrm{l}}\sim \tau ^{1/2}$
in the limit as
$\tau \rightarrow 0,$
this implies that the aforementioned layer is of width
$\tau$
in physical coordinates.
To continue the computation for
$\tau \gt 0,$
there are two possibilities, each with its own advantages and drawbacks. One can continue the integration in the Landau space; however, this means that it would not be possible to extend the method to two or three dimensions, which is the ultimate goal. The alternative, and the one we adopt here, is to use the small-time solution obtained from the Landau-based approach to set an initial condition at
$\tau =\tau _{0}\gt 0$
, where
$\tau _{0}$
is taken to be arbitrarily small, and to integrate the equations in the physical space, using a deformed-mesh approach. What we have in mind is the Arbitrary Lagrangian-Eulerian (ALE) formulation within the finite-element software Comsol Multiphysics, which was used in [Reference Vynnycky37, Reference Vynnycky and Kimura38] to solve the problem of solidification of a pure material in the presence of natural convection in a two-dimensional cavity; although that problem involved only one moving boundary, it was still demanding in its own right, since the full Navier Stokes equations had to be solved in the liquid region.
Even so, the situation is, in some ways, more complicated than that in [Reference Vynnycky37, Reference Vynnycky and Kimura38]. Here, there are initially three domains: solid, mush and melt. After some time, only solid and mush is left; after that, only solid is left. Thus, although we start by solving equations (2.31)–(2.38), subject to (2.39)–(2.52), we transition to solving (2.31)–(2.35), subject to (2.39)–(2.42) and (2.46), which is when only solid and mush are left. Once the mush disappears, the macrosegregation profile is determined and there is little point in solving any further.
The curves in (
$\beta,C_{0})$
-space which determine whether, initially,
$\ \chi _{\text{s}}=0$
or
$\chi _{\text{s}}\gt 0$
when
$R=1,3$
for: (a)
$\Delta T=1$
K; (b)
$\Delta T=101$
K. Results obtained using the numerical method in [Reference Assunção, Vynnycky and Mitchell6].
One might question why, if one is going to integrate into the physical space, it was necessary to perform the small-time analysis in the Landau space at all; after all,
$\tau _{0}$
has to be chosen arbitrarily anyway. The subtle reason for this is that, for a given set of parameters (
$\beta,C_{0},R$
), it is not in general obvious
$a$
$priori$
which of the two choices in (2.40) is the appropriate one at
$\tau =0.$
When
$R=1,$
however, it is possible to delineate analytically in
$ ( \beta,C_{0} )$
-space which is the appropriate condition; this was given in Figure 4 in [Reference Assunção, Vynnycky and Mitchell6]. By way of example, Figure 3 shows the curves in
$ ( \beta,C_{0} )$
-space which determine whether
$\chi _{\text{s}}=0$
or
$\chi _{\text{s}}\gt 0$
for
$R=1$
and 3, for two different values of
$\Delta T:=T_{\text{cast}}-T_{\text{liq}},$
which is often termed as the superheat; here, as in earlier work [Reference Voller31], one of the values of
$R$
is chosen to be unrealistically high in order to demonstrate the trend for
$R\gt 1$
. In addition, note that the location of this curve can also be affected by the other dimensionless parameters in the problem, i.e.
$\tilde{m},St,\theta _{\text{m}};$
however, once a particular binary alloy has been chosen, this can only mean a dependency on the dimensional parameters
$T_{\text{cast}}$
and
$T_{\text{w}}.$
We also remark that, if
$R=1,$
none of the other parameters can affect the curves in Figure 3, which are given for this case by
\begin{align} C_{0}=C_{\mathrm{eut}}\!\left ( \beta k_{0}\right ) ^{\frac{1-k_{0}}{1-\beta k_{0}}}. \end{align}
From another point of view, Figure 4 shows
$\chi _{\mathrm{s}}$
as a function of
$\beta$
for three values of
$C_{0}$
,
$R=1,3$
and
$\Delta T=101$
K; when
$R=1,$
we have, from [Reference Assunção, Vynnycky and Mitchell6],
\begin{align}{\chi _{\mathrm{s}}=\max \!\left ( 0,\frac{1}{1-\beta k_{0}}\left \{ \left ( \frac{C_{0}}{C_{\mathrm{eut}}}\right ) ^{\frac{1-\beta k_{0}}{1-k_{0}}}-\beta k_{0}\right \} \right ) .} \end{align}
Moreover, it is perfectly possible that, as the integration in
$\tau$
proceeds, the appropriate condition at
$Y=Y_{\mathrm{s}} ( \tau )$
changes from the first option in (2.40) to the second, i.e. from non-eutectic solidification to eutectic solidification, although never the other way around.
$\chi _{\mathrm{s}}$
at
$\tau =0$
as a function of
$\beta$
for three values of
$C_{0}$
: (a)
$R=1$
; (b)
$R=3$
. Results obtained using the numerical method in [Reference Assunção, Vynnycky and Mitchell6].
Details of the numerical implementation, as well as code verification tests, are documented in Appendix B. Here, instead, we proceed to the results concerning solidification and macrosegregation profiles.
4. Results
Here, we focus on results obtained using the data in Table 1. Figure 5 shows
$Y_{\mathrm{s}},Y_{\mathrm{l}}$
and
$Y_{\infty }$
as functions of
$\tau$
for
$R=1$
and
$\beta =0,1;$
this gives upper and lower bounds, with respect to
$\beta,$
on
$Y_{\mathrm{s}}$
and
$Y_{\mathrm{l}}$
for the case when there is no solidification shrinkage. Figure 6 shows the corresponding solutions when
$R=$
1.3; this was the unrealistically high value chosen in earlier work [Reference Assunção, Vynnycky and Mitchell6, Reference Voller31], with a view to highlighting the effect of solidification shrinkage. In both cases, the computation is stopped when complete solidification is reached, at which point
$Y_{\mathrm{l}}=Y_{\infty }=1/R.$
From these, we can note that our numerical method is adequately able to capture the fact that, first, the melt-only region disappears, i.e when
$Y_{\mathrm{l}}=Y_{\infty }$
, and also when the mush region disappears, i.e. when
$Y_{\mathrm{s}}=Y_{\infty }$
.
$Y_{\mathrm{s}},Y_{\mathrm{l}}$
and
$Y_{\infty }$
as functions of
$ \tau$
for
$R=1$
and: (a)
$ \beta =0$
; (b)
$ \beta =1$
.
$Y_{\mathrm{s}},Y_{\mathrm{l}}$
and
$Y_{\infty }$
as functions of
$ \tau$
for
$R=1.3$
and: (a)
$ \beta =0$
; (b)
$ \beta =1$
.
Of primary metallurgical interest is the final solute profile after complete solidification, which we will denote by
$C_{\text{solid}};$
using (2.6), this is extracted from the computations by calculating
$\mathcal{C}_{\text{solid}},$
where
\begin{align}{R\mathcal{C}_{\text{solid}}=\chi _{\text{s}}\!\left . \mathcal{C}_{\mathrm{l}}\right \vert _{Y=Y_{\text{s}}}+R\int _{0}^{1-\chi _{\text{s}}}\!\left ( \mathcal{C}_{\mathrm{s}}+\beta \chi ^{\prime }\frac{\mathrm{d}\mathcal{C}_{\mathrm{s}}}{\mathrm{d}\chi ^{\prime }}\right ) \mathrm{d}\chi ^{\prime },} \end{align}
with
$C_{\text{solid}}=C_{0}{\mathcal{C}_{\text{solid}}.}$
Figure 7 shows plots of
$C_{\text{solid}}$
as a function of
$Y$
for
$\beta =0$
and 1, with
$R=1,1.1,1.2,1.3$
. These display a number of salient features. First of all, as expected, the profiles become less uniform as
$R$
is increased. Moreover, the plots show the well-known phenomenon of inverse macrosegregation, whereby the solute concentration is higher at the cooling surface at
$Y=0$
. However, it is perhaps surprising to see that the macrosegregation level remains fairly constant for a sizeable interval in
$Y;$
moreover, this interval corresponds to the time during which there is still a melt-only region, i.e.
$Y_{\mathrm{l}}\lt Y_{\infty }$
, as verified by cross-checking this figure against Figure 6. This is done in Figure 8, although it is evident that the trend is much more distinct for the case when
$\beta =0$
. A likely factor in this is that
$\mathcal{C}{_{\text{l}}}$
in the melt-only region retains its initial value until this point in time, as is readily demonstrated using the method of characteristics, as follows. From (2.36) and (2.38), we obtain, for
$Y_{\mathrm{l}}(\tau )\leq Y\leq Y_{\infty }(\tau ),$
\begin{align} \frac{\partial \mathcal{C}_{\mathrm{l}}}{\partial \tau }+\hat{U}\frac{\partial \mathcal{C}_{\mathrm{l}}}{\partial Y}=0, \end{align}
with
$\hat{U}=\hat{U} ( \tau ) .$
Thence,
$\mathcal{C}_{\mathrm{l}}=\mathcal{C}_{\mathrm{l}} ( Y-\int _{0}^{\tau }\hat{U} ({\tau ^{\prime }} ){\mathrm{d}\tau ^{\prime }} ) ;$
applying equation (2.48) then gives
$\mathcal{C}_{\mathrm{l}}=1,$
and hence
$C_{\mathrm{l}}=C_{0},$
as claimed, In turn, this gives equation (2.17). This will mean that
$C_{\mathrm{l}}=C_{0}$
on all characteristics that enter the mush, which makes it more likely that the values of
$C_{\mathrm{l}}$
at
$Y=Y_{\mathrm{s}}(\tau )$
will coincide with each other; this will mean that the values of
$C_{\text{solid}}$
also coincide with each other. Of course, this is not completely guaranteed, since the value of
$C_{\mathrm{l}}$
at
$Y=Y_{\mathrm{s}}(\tau )$
is determined by the solution to equation (2.34), with
$\chi$
and
$\hat{U}$
both varying on different characteristics. This line of argument is summarised in Figure 9.
$C_{\mathrm{solid}}$
as a function of
$Y$
for
$R=1,1.1,1.2,1.3$
and: (a)
$\beta =0$
; (b)
$\beta =1$
.
Superposition of
$Y_{\mathrm{l}}(\tau )$
and
$C_{\mathrm{solid}}$
at
$Y=Y_{\mathrm{s}}$
as functions of
$\tau$
for
$R=1.3$
and: (a)
$ \beta =0$
; (b)
$ \beta =1$
.
Moreover, the more or less constant value for
$C_{\text{solid}}$
suggests that the similarity solution given in [Reference Voller31] would give an accurate value for the macrosegregation level until this point. However, we should note that for
$\beta =0$
and for higher values of
$R,$
if anything, the value of
$C_{\text{solid}}$
actually increases slightly from its value at
$Y=0.$
This is also seen in Figures B2 and B3; thereafter,
$C_{\text{solid}}$
decreases. However, what happens after that depends on the value of
$\beta$
: for
$\beta =1,$
we have negative segregation as complete solidification is approached, i.e.
$C_{\mathrm{solid}}$
is locally less than the initial composition,
$C_{0}$
, but the profile for
$C_{\text{solid}}$
starts to increase again for
$\beta =0,$
leading to positive segregation, i.e.
$C_{\mathrm{solid}}\gt C_{0}$
. The differing trends may be considered reasonable on physical grounds.
$\beta =1$
corresponds to the assumption of solute diffusion in the solid at the microscale, and this has resulted in solute depletion as full solidification is reached; hence,
$C_{\text{solid}}$
decreases to ensure global solute conservation. On other hand,
$\beta =0$
corresponds to the assumption of zero solute diffusion in the solid at the microscale; this results in the build-up of solute in the remaining mush, which manifests itself in the increase in
$C_{\text{solid}}$
as full solidification is approached.
Finally, we comment that, although the Al–Cu system is often modelled with the Scheil equation at the microscale, i.e.
$\beta =0,$
it is nevertheless instructive to produce results with the lever rule, i.e.
$\beta =1,$
in order to see what effect this assumption has on the solution. Moreover, the parameter
$\beta$
is merely a way to parametrise microsegregation, with 0 and 1 being the extreme values. Since the appropriate value of
$\beta$
for any given binary alloy system is not actually known, results for the extreme values of
$\beta$
can serve to provide bounds on what may happen in reality; however, additional results for intermediate values of
$\beta$
are given in Appendix C.
5. Discussion
We now discuss the above results in two contexts: earlier attempts to consider this problem; relevance for the modelling of macrosegregation in continuous casting.
5.1. Earlier results on inverse segregation
Inverse segregation induced by solidification shrinkage has previously been considered both theoretically and experimentally on a number of occasions [Reference Chen and Tsai10, Reference Kato and Cahoon21, Reference Kirkaldy and Youdelis22, Reference Minakawa, Samarasekera and Weinberg25, Reference Rousset, Rappaz and Hannart29, Reference Voller31, Reference Voller and Sundarraj34], although not always resulting in a self-consistent macrosegregation profile over the entire solidified sample, in the sense of solute conservation. One of the exceptions is [Reference Rousset, Rappaz and Hannart29], wherein the solidification of an Al-3.4 wt% Cu ingot, for which
$R=$
1.07 and
$\beta =0,$
was considered theoretically and experimentally, although comparatively few mathematical details were given; the results were recently re-capped by Dantzig and Rappaz [Reference Dantzig and Rappaz12, p. 582] and are, for convenience, shown in Figure 10. However, we have not attempted to reproduce these results with our method, for a number of reasons. First of all, the effect of microporosity was considered simultaneously with solidification shrinkage, whereas we have only focused on the effect of the latter. More significantly, however, the situation considered was not that for isothermal cooling, which is in practice difficult to achieve experimentally anyway, but where the cooling boundary condition that would replace (2.15) would be rather on the heat flux,
$-k\partial T/\partial y.$
This change would lead to qualitatively different behaviour for
$y_{\mathrm{l}},y_{\mathrm{s}}$
and
$y_{\infty },$
as compared to that shown in Figure 6; the new behaviour is depicted in Figure 11. This time, the mushy zone will not start to form instantaneously, but only after a delay time, during which the temperature at
$y=0$
decreases to
$T_{\mathrm{m}}-mC_{0};$
prior to this time, there cannot be any solidification shrinkage, meaning that
$Y_{\infty }$
does not decrease from its original value, if the cooling shrinkage is negligible. In fact, the numerical method presented here would not be able to capture this short-time evolution, as it has been developed for the case when solid and mush form instantaneously, and an alternative small-time asymptotic solution would be required. Next, after the temperature at
$y=0$
has decreased to
$T_{\mathrm{m}}-mC_{0},$
there will be a mushy zone and a melt-only region, until the temperature at
$y=0$
is low enough for only solid to form there; at this stage, the temperature there may either be
$T_{\text{eut}}$
or a higher temperature that lies on the solidus line in Figure 2. In view of these differences, note that a plot corresponding to Figure 7 would be generated from results obtained only after solid starts to form at
$y=0,$
and we now proceed to a detailed discussion of this.
Experimental points and simulated profile of copper concentration in an Al-3.4 wt % Cu ingot, as measured and computed in [Reference Rousset, Rappaz and Hannart29], from which the figure has been adapted.
Qualitative schematic for
$Y_{\mathrm{l}},Y_{\mathrm{s}}$
and
$Y_{\infty }$
as functions of
$ \tau$
when a heat flux condition is applied at the cooling boundary.
Earlier work by Flemings and Nereo [Reference Flemings and Nereo18, Figure 6, p. 1455] suggested that the profile for
$C_{\text{solid}}$
would consist of a central plateau at which
$C_{\text{solid}}=C_{0},$
surrounded by a region of inverse segregation near
$Y=0,$
where
${C_{\text{solid}}}\gt C_{0},$
and a region of negative segregation where
$C_{\text{solid}}\lt C_{0},$
corresponding to the last part of the sample to solidify. Dantzig and Rappaz [Reference Dantzig and Rappaz12, p. 582] comment that the plateau obtained, both experimentally and theoretically, by Rousset et al. [Reference Rousset, Rappaz and Hannart29] is slightly higher than that of the nominal composition, as seen in Figure 10; note, however, that the curve for initial composition was plotted in [Reference Dantzig and Rappaz12, Reference Rousset, Rappaz and Hannart29] as being much higher than the actual value of 3.4 wt% Cu, and we have therefore corrected this in Figure 10. Of relevance to this discussion are also the results of Voller [Reference Voller31], mentioned earlier, which were for a semi-infinite sample, but which qualitatively agree better with our own findings, at least as far as the region near
$Y=0$
is concerned, in that the plateau is considerably higher than the nominal composition,
${C}_{0}$
. Moreover, as regards the last part to solidify, Rousset et al. [Reference Rousset, Rappaz and Hannart29] find that the theoretically obtained profile for
$C_{\text{solid}}$
decreases rapidly, whereas the experimental results show a decrease, followed by an increase; this resembles the results that we have obtained for a variety of values of
$R$
for
$\beta =1,$
but not for
$\beta =0,$
although we should once again emphasise that our computations were for a different cooling boundary condition at
$Y=0.$
Moreover, we point out that we have been very thorough in verifying our numerical method, as seen in Appendix B.
To summarise the discussion, we show in Figure 12, in as far as it is possible, a qualitative schematic of how we expect the profiles for
$C_{\text{solid}}$
to differ when fixed-temperature and fixed-flux conditions are used; plot (a) is for when
$\beta$
is close to zero, where plot (b) is for when
$\beta$
is close to 1. Needless to say, these plots would vary depending on the numerical values of the model parameters, but it is nevertheless possible to make some quite definitive qualitative predictions and observations:
Qualitative schematic for
$C_{\text{solid}}$
as a function of
$Y$
when fixed-temperature and fixed-flux boundary conditions are used, with: (a)
$ \beta \approx 0$
; (b)
$ \beta \approx 1.$
Note that
$C_{\text{solid}}$
at
$Y=0$
is intentionally greater in (b)
$.$
Moreover,
$C_{\text{solid}}$
at
$Y=0$
need not in practice be the same for the two types of cooling conditions.
-
1. In both cases,
$C_{\text{solid}}\gt C_{0}$
at
$Y=0,$
since inverse segregation occurs regardless of the thermal boundary condition. -
2. In the region immediately adjacent to
$Y=0,$
$C_{\text{solid}}$
will be more or less constant for the fixed-temperature condition, but will decrease from its value at
$Y=0.$
-
3. As
$Y$
increases further,
$C_{\text{solid}}$
will decrease from its more or less constant value in the case of the fixed-temperature condition. For the fixed-flux condition,
$C_{\text{solid}}$
may either plateau out or still continue to decrease. -
4. Thereafter, it becomes harder to draw conclusions. In both cases, solute must be conserved. For the fixed-temperature condition,
$C_{\text{solid}}$
has hitherto had a value greater than
$C_{0},$
which means that it has to decrease below this value; this is seen in both plots in Figure 7, although the values of
$C_{\text{solid}}$
become so low for
$\beta =0$
that the profile has to compensate so that
$C_{\text{solid}}\gt C_{0}$
near
$Y=1/R$
in order to achieve solute conservation, all the more so for higher values of
$R.$
For the fixed-flux condition, one would not expect as dramatic a difference in the behaviour of
$C_{\text{solid}}$
near
$Y=1/R$
for different values of
$\beta,$
since
$C_{\text{solid}}$
has not strayed from
$C_{0}$
for as large an interval in
$Y$
as was the case for the fixed-temperature condition. Note, for example, that there is negative segregation near
$Y=1/R$
even for
$\beta =0$
in Figure 10, which is a fixed-flux case, whereas we obtained this kind of negative segregation for a fixed-temperature case for
$\beta =1.$
5.2. Relevance to continuous casting
An important question is whether this model can be used to predict what happens as regards centreline segregation in steady-state continuous casting processes, in particular those for steel [Reference Fredriksson and Åkerlind20, Reference Rogberg28]; this is a severe form of positive segregation near the axis of symmetry of the solidified casting. A simplified schematic for the continuous casting process is given in Figure 13. Typically, molten steel enters a cooling mould region of width
$2W$
at
$z=0$
, solidifies and is withdrawn with a constant casting speed
$V_{\mathrm{cast}}$
. Whereas
$W$
is of order of tens of centimetres, the distance from the top of the cooling mould to the end of solidification at the centreline can be as great as 20 m; consequently, the geometry for continuous casting can be considered as slender. One reason to suppose that the two situations might be related is that, because the geometry for continuous casting is slender, the governing equations for a two-dimensional (2D) steady-state model for the process have a great number of similarities with those analysed in this paper [Reference Vynnycky and Kimura39, Reference Vynnycky and Saleem40]; in particular, the coordinate in the casting direction,
$z,$
acts as a time-like variable via the relation
$z=V_{\mathrm{cast}}t$
. Furthermore, the fact that
$y=W$
is an axis of symmetry means that only the left-half of the geometry need be considered and that the boundary conditions there for
$T$
and
$C_{\mathrm{l}}$
will be the same as those given in (2.16) and (2.17), respectively. Moreover, the region near the centreline in continuous casting is always the last to solidify and, as seen from Figure 7, the region that solidifies last tends to undergo the greatest changes in solute concentration. Indeed, the analogy between the two situations is appropriate, but only as regards the equation for the conservation of heat, (2.3), and when
$R=1$
. On the other hand, there are significant differences:
A schematic of a continuous casting process for an alloy.
-
the extent of the spatial domain in
$y$
is fixed in the 2D model for continuous casting; -
in the latter, there are two components of velocity, i.e. in the
$y$
- and
$z$
-directions, whereas there is only one component in the one-dimensional (1D) model; -
Darcy’s law for flow in the mushy zone is a key part of the 2D continuous casting problem, and the pressure cannot be solved for a posteriori. In particular, this means that the expression used for the mush permeability will affect the solution, and hence the macrosegregation, which was not the case for the time-dependent 1D model.
6. Conclusions
In this paper, we have extended our earlier work for the one-dimensional solidification of a binary alloy undergoing shrinkage [Reference Assunção, Vynnycky and Mitchell6]. Whereas the earlier work considered the small-time similarity solution to this problem, here we have used that solution as the initial condition for the full problem, which we compute numerically until complete solidification. This starts as a triple moving-boundary problem, which turns into a dual moving-boundary problem when the melt-only phase depletes. Of particular interest is the final macrosegregation profile that is obtained. Near the cooled surface, there is an elevated solute concentration, corresponding to the expected phenomenon of inverse macrosegregation. Thereafter, this elevated level is maintained until the melt-only phase depletes; thereafter, the solute concentration decreases appreciably, resulting in a region of negative segregation, i.e. where the local concentration is lower than the starting concentration,
$C_{0}.$
What happens after that depends on the value of the microsegregation parameter,
$\beta :$
greater values of
$\beta,$
with
$\beta =1$
corresponding to the lever rule, result in continued negative segregation, whereas lower values of
$\beta,$
with
$\beta =0$
corresponding to the Scheil equation, result in positive segregation, i.e. the local concentration is higher than
$C_{0}$
, as complete solidification is approached, particularly for higher values of the shrinkage parameter,
$R.$
The formulation presented in this paper can serve as a starting point for extension in several directions. Perhaps the most obvious would be the implementation of a convective (Robin) boundary condition for the temperature, instead of the Dirichlet condition used in equation (2.15), as this is more likely the case in practice. Furthermore, our formulation can be useful in informing on the development of a model for macrosegregation in continuous casting processes. Also of significance is that our approach gives a framework for studying macrosegregation in ingot casting using a moving mesh formulation, rather than the more commonly used enthalpy method. As a consequence, this avoids the need to use arbitrary artificial parameters to ensure numerical convergence [Reference Amberg2, Reference Assunção and Vynnycky5]; indeed, the only arbitrary parameter that we have had to introduce was
$\tau _{0},$
which is related to the time after the start of solidification.
Another direction for future work would be to use the current formulation to shed further light on the formation of channel segregates in casting processes, i.e. freckles, A-segregates, V-segregates [Reference Vynnycky, Saleem and Fredriksson41]; whereas the first two are believed to form as a consequence of the enrichment of the interdendritic melt with light solute elements, leading to a decrease in the local melt density and the onset of thermosolutal convection [Reference Auburtin, Wang, Cockcroft and Mitchell7], a mechanism for the third has still to be determined. Although the formation of channel segregates is a comparatively old problem [Reference Anderson and Guba3, Reference Emms and Fowler16, Reference Fowler19, Reference Mehrabian, Keane and Flemings24, Reference Worster43], there is still considerable doubt as to whether existing numerical simulations are able to compute them correctly [Reference Vušanović and Voller35, Reference Vušanović and Voller36], as regards mesh independence; this refers primarily to their width, their length and the spacing between them. In particular, the approach adopted here could elucidate when and exactly where in the mushy zone they are initiated; a precursor to this is believed to be when remelting first occurs, i.e. where
$\partial \chi/\partial t\gt 0$
locally. Although solidification shrinkage is not responsible for this, it may nevertheless act to moderate it, by delaying the onset of thermosolutal convection, if not preventing it completely [Reference Monde, Shrivastava, Jakhar and Chakraborty27].
Acknowledgements
The authors would also like to thank the anonymous reviewers for their important comments and suggestions.
Conflicts of interest
None.
Appendix A: derivation of equation (2.8)
First, note that the two extremes for solute transport at the microscale are the lever rule and the Scheil equation [Reference Swaminathan and Voller30, Reference Voller and Sundarraj34], where
\begin{align} \rho C=\chi \rho _{\mathrm{l}}C_{\mathrm{l}}+\quad \left \{ \begin{array}{l@{\quad}l}\rho _{\mathrm{s}}\!\left ( 1-\chi \right ) C_{\mathrm{s}},\quad & \text{lever rule}\\[5pt] \rho _{\mathrm{s}}\int _{0}^{1-\chi }C_{\mathrm{s}}\mathrm{d}\chi ^{\prime },\quad & \text{Scheil equation}\quad \end{array} \right .\!\!,\quad \quad \end{align}
with
$C_{\mathrm{l}}$
and
$C_{\mathrm{s}}$
as the concentrations of the solute in the liquid and solid phases, respectively, which are related by
\begin{align} C_{\mathrm{s}}=k_{0}C_{\mathrm{l}},\quad \end{align}
where
$k_{0}$
is the partition coefficient; moreover, the lever rule assumes thermodynamic equilibrium between the phases, whereas the Scheil equation assumes no solute diffusion in the solid and perfect mixing in the liquid [Reference Alexiades and Solomon1]. For the purposes of generalisation, it is possible to introduce a parameter
$\beta,\,$
where
$0\leq \beta \leq 1, $
that allows for a back-diffusion treatment, i.e. partial solute diffusion into the solid, that lies between the limits of zero back diffusion (
$\beta =0$
, the Scheil assumption) and complete back diffusion (
$\beta =1$
, the lever rule); as indicated in [Reference Swaminathan and Voller30], this treatment of back diffusion is equivalent to using the Clyne and Kurz correction [Reference Clyne and Kurz11] of the well-known back diffusion model of Brody and Flemings [Reference Brody and Flemings9]. Following [Reference Swaminathan and Voller30] in first setting
\begin{align} \int _{0}^{1-\chi }\frac{\partial C_{\mathrm{s}}}{\partial t}\mathrm{d}\chi ^{\prime }=\beta \!\left ( 1-\chi \right ) k_{0}\frac{\partial C_{\mathrm{l}}}{\partial t},\quad \end{align}
wherein equation (A2) has been used to eliminate
$C_{\mathrm{s}}$
, equation (2.5) becomes (2.8); furthermore, we note that (A1) can now be generalised for
$0\leq \beta \leq 1$
to equation (2.6).
Appendix B: code verification
Equations (2.31)–(2.38), subject to (2.39)–(2.52), were solved using the transient solver in the finite element-based software, Comsol Multiphysics, which has the inbuilt infrastructure to cope with this system of parabolic and hyperbolic PDEs. Quadratic Lagrangian elements were used to discretize the domain, in combination with the software’s deformed mesh mode, whereby an arbitrary Lagrangian–Eulerian (ALE) formulation is employed in order to solve free or moving boundary problems. As explained in Section 3, different equations are solved at different times of the computation, and the calculation is stopped when the mush disappears, as the final macrosegregation profile is known by that stage. Numerically, this evolution means that mesh-folding occurs, as first the melt region disappears, and then later so does the mush. Denoting the syntax of the software with typewriter font below, the specification of these models requires the use of customised equations, available through three General Form PDE interfaces, to solve the second-order PDEs (2.31), (2.33) and (2.37), and 5 Weak Form PDE interfaces, to solve first-order PDEs and algebraic equations (2.32), (2.34)–(2.36) and (2.38). In addition, two Global ODEs and DAEs interfaces, where DAE denotes differential-algebraic equation, are used to determine the time-dependent position of interfaces
$Y=Y_{\mathrm{s}}(\tau )$
and
$Y=Y_{\mathrm{l}}(\tau )$
, by solving equations (2.41) and (2.44). Note also that
$Y_\infty (\tau )$
is determined indirectly from (2.14), which gives
$\dot{Y}_\infty (\tau )$
. The remaining initial, boundary and interfacial conditions are passed directly through the three General Form PDE interfaces. As regards the geometry deformation caused by the moving boundaries at
$Y=Y_{\mathrm{s}}(\tau ),Y_{\mathrm{l}}(\tau )$
and
$Y_\infty (\tau )$
, this is handled using a Moving Mesh interface. For this, Comsol Multiphysics offers the options of Laplace smoothing or Winslow smoothing [Reference Winslow42]; however, it is well-established that the latter is more effective [Reference Knupp23], and that was our experience also, in the sense that solutions could be obtained until the melt and mush regions were practically extinct, as seen from the results in Figures 5-7.
For all cases, the same convergence criterion at each time step was used, namely,
\begin{align} \left ( \frac{1}{N_{\mathrm{dof}}}\sum _{i=1}^{N_{\mathrm{dof}}}\left ( \frac{\left \vert E_{i}\right \vert }{A+\mathcal{R}\!\left \vert U_{i}\right \vert }\right ) ^{2}\right ) ^{\frac{1}{2}}\lt 1, \end{align}
where
$E_{i}$
is the solver’s estimate of the absolute error in the latest approximation to the
$i^{\text{th}}$
component of the scaled solution vector,
$U_{i}$
, at each time step,
$A$
is the absolute tolerance,
$\mathcal{R}$
is the relative tolerance and
$N_{\mathrm{dof}}$
is the number of degrees of freedom (DOF); for simplicity, we take
$\mathcal{R}=A,$
and denote this common value by
$\varepsilon .$
Note also that
$N_{\mathrm{dof}}$
is related to the number of elements,
$N,$
that are used to discretize the computational domain for values of time when solid, mush and melt are all present; thus, in the course of a computation, the values of
$N$
and
$N_{\mathrm{dof}}$
decrease, as first the melt and then the mush disappear.
Before performing batch runs, there were a number of issues to resolve:
-
1. what values of
$\tau _{0}$
are permissible? -
2. how many elements are required for mesh-independent results?
-
3. what values of
$\varepsilon$
should be used? -
4. how well is solute conserved in the final macrosegregation profile
$?$
To some extent, the issues are interlinked, but as a starting point for numerical experiments, we began with
\begin{equation*} \tau _{0}=10^{-2},\quad N=100,\quad \varepsilon =10^{-2}, \end{equation*}
which is a combination that gives least accuracy, but also requires least computing time; this value of
$N$
corresponds to
$N_{\mathrm{dof}}=1511.$
We start by considering the effect of
$\tau _{0}.$
Since
$Y_{\mathrm{l}},Y_{\mathrm{s}}\sim \tau ^{1/2}$
for small
$\tau,$
it is clear that even if we take
$\tau _{0}=10^{-2}\ll 1,$
then of the order of 10%, or much more if
$R\gt 1,$
of the macrosegregation profile will be poorly computed. Figure B1 shows
$Y_{\mathrm{l}}$
and
$Y_{\mathrm{s}}$
as functions of
$\tau$
for different values of
$\tau _{0};$
these results were obtained for
$\beta =0$
and
$R=1.3,$
with
$\varepsilon =10^{-2}.$
Moreover, although Figure B1 may suggest that the choice of
$\tau _{0}$
does not particularly affect the profiles of
$Y_{\mathrm{l}}$
and
$Y_{\mathrm{s}}$
, it is nevertheless advisable to use the lowest value possible of
$\tau _{0},$
in view of the need to compute
$C_{\text{solid}}$
correctly
$;$
thus, all subsequent computations were carried out with
$\tau _{0}=10^{-4}.$
$Y_{\mathrm{s}}$
and
$Y_{\mathrm{l}}$
as functions of
$\tau$
for
$R=1.3,\beta =0,\varepsilon =10^{-2},N=100$
for
$\tau _{0}=10^{-4},10^{-3},10^{-2}$
: (a)
$Y_{\mathrm{s}}$
; (b)
$Y_{\mathrm{l}}$
.
$C_{\mathrm{solid}}$
as a function of
$Y$
for five different meshes
$,$
with
$R=1.3$
,
$\beta =0,\tau _{0}=10^{-4},\varepsilon =10^{-3}$
.
$C_{\mathrm{solid}}$
as a function of
$Y$
for
$\varepsilon =10^{-4},10^{-3},10^{-2},$
with
$R=1.3,\beta =0,\tau _{0}=10^{-4},N=1600$
.
Moving onto the issue of mesh independence, Figure B2 shows
$C_{\text{solid}}$
as a function of
$Y$
for
$\beta =0$
and
$R=1.3,$
and five different meshes; for these runs,
$\tau _{0}=10^{-4},\varepsilon =10^{-3}$
. At this lower value of
$\varepsilon,$
the profiles appear to be more or less on top of each other, indicating that
$N=100$
is large enough.
Figure B3 shows
$C_{\text{solid}}$
as a function of
$Y$
for
$\beta =0$
and
$R=1.3,$
and three different values of
$\varepsilon ;$
for these runs,
$\tau _{0}=10^{-4}$
and
$N=1600,$
which corresponds to
$N_{\mathrm{dof}}=24011.$
Here also, the profiles are more or less on top of each other, indicating that
$\varepsilon =10^{-2}$
is small enough, albeit at such elevated values for
$N;$
however, lower values
$N$
did not lead to such good agreement.
Lastly, we consider how well solute is conserved numerically by the end of the computation. Defining
\begin{align} C_{\mathrm{total}}\!\left ( Y\right ) =R^{-1}\int _{0}^{Y}{C_{\text{solid}}\mathrm{d}Y}^{\prime }, \end{align}
we note that we should expect from (2.52) that
$C_{\mathrm{total}} ( R^{-1} ) =C_{0}.$
Figure B4 shows
$C_{\mathrm{total}}$
for different values of
$R,$
using
$\tau _{0}=10^{-4},N=100,\varepsilon =10^{-2};$
$C_{\mathrm{total}}$
has been calculated using (4.1)
$.$
From Figure B4, we thus find that solute is conserved adequately enough, even though the values of
$N$
and
$\varepsilon$
are relatively small and large, respectively; we surmise that this is because
$C_{\mathrm{total}}$
is a global quantity, and would therefore not be as subject to local discrepancies as
$C_{\text{solid}}.$
$C_{\mathrm{total}}$
as a function of
$Y$
for
$\tau _{0}=10^{-4},N=100,\varepsilon =10^{-2},$
with
$R=1.3$
,
$\beta =0$
.
Lastly, we note that computations for which
$\tau _{0}=10^{-4},N=100,\varepsilon =10^{-2}$
were found to require less than 300 s on a laptop with an Intel Core i7-9750HQ CPU 2.60 GHz processor and 16 GB RAM. More details are given in [Reference Assunção4].
Appendix C: results for intermediate values of
$\beta$
To complement the results of Figure 7, which were for
$\beta =0$
and 1, Figure C1 shows
$C_{\mathrm{solid}}$
as a function of
$Y$
for
$R=1,1.1,1.2,1.3$
for the intermediate value of
$\beta =1/2$
. Over the two figures, it is evident how, as the value of
$\beta$
is increases, the value of
$C_{\mathrm{solid}}$
in the plateau region extending from
$Y=0$
decreases; in addition, the profiles transition from displaying negative segregation in the interior and positive segregation at far right towards just negative segregation at far right. These trends persist for all three values of
$R\gt 1$
. Note, however, the degree of positive segregation at far right is higher for
$\beta =1/2$
than for
$\beta =0$
; one might have expected it to have been lower in the transition towards negative segregation for
$\beta =1$
. The explanation for this seems to come from the fact that the negative segregation in the central portion, in relation to the positive (inverse) segregation near
$Y=0$
, is much greater for
$\beta =1/2$
than for
$\beta =0$
; thus, a greater positive segregation is necessary at far right in order to ensure global solute conservation.
$C_{\mathrm{solid}}$
as a function of
$Y$
for
$R=1,1.1,1.2,1.3$
and
$\beta =1/2$
.
Open access
