Solidification Simulation of Al-Si Alloys with Dendrite Tip Undercooling

A novel solution approach is proposed for the numerical simulation of the solidification process of binary Al-Si hypoeutectic alloys during upward and downward solidification modes. Undercooling is always observed during solidification, but the phenomenon could not be considered in the present-day numerical solution. In this approach, the temperature distribution in the mushy zone was used to define the fraction of solid, which enabled the evaluation of the effect of dendrite tip undercooling on the characteristics of the binary alloy solidification. The present numerical algorithm was found to significantly reduce the computation time. Transient temperature distribution and solidification time from the numerical analysis, with consideration of natural convection due to temperature and concentration gradients, have been successfully simulated and validated with experiment results. Numerical results with consideration of dendrite tip undercooling have better agreement with experimental results. The effect of dendrite tip undercooling on the fluid flow (velocity profile), G, R and λ1 for both upward and downward solidification modes of Al-Si alloys have been investigated and discussed. Consideration of undercooling was found to increase G and reduce R in both solidification modes. During downward solidification, considering undercooling significantly increased flow velocity and decreased λ1. The primary dendrite arm spacing could be validated with results from uni-directional solidification experiments only when dendrite tip undercooling was considered.


Abstract:
A novel solution approach is proposed for the numerical simulation of the solidification process of binary Al-Si hypoeutectic alloys during upward and downward solidification modes. Undercooling is always observed during solidification, but the phenomenon could not be considered in the present-day numerical solution. In this approach, the temperature distribution in the mushy zone was used to define the fraction of solid, which enabled the evaluation of the effect of dendrite tip undercooling on the characteristics of the binary alloy solidification. The present numerical algorithm was found to significantly reduce the computation time. Transient temperature distribution and solidification time from the numerical analysis, with consideration of natural convection due to temperature and concentration gradients, have been successfully simulated and validated with experiment results. Numerical results with consideration of dendrite tip undercooling have better agreement with experimental results. The effect of dendrite tip undercooling on the fluid flow (velocity profile), G, R and λ 1 for both upward and downward solidification modes of Al-Si alloys have been investigated and discussed. Consideration of undercooling was found to increase G and reduce R in both solidification modes. During downward solidification, considering undercooling significantly increased flow velocity and decreased λ 1 . The primary dendrite arm spacing could be validated with results from uni-directional solidification experiments only when dendrite tip undercooling was considered.

Introduction
Solidification is the phenomenon of phase transformation from liquid phase to solid phase. Solidification mechanisms for pure metals and alloys are different from each other. Solidification for pure metals is mainly dominated by heat-transfer conditions [5,6]. However, both heat and mass transfer conditions affect the solidification of metallic alloys. During solidification of binary alloys, there are three distinct zones, namely, the liquid, solid and the intermediate two-phase (solid and liquid) region termed the mushy zone [7,8]. During the past few decades, there have been several attempts to simulate solidification of binary alloys [5,6,[9][10][11][12][13][14]. These simulation attempts can be broadly classified into three categories: • Solution to energy equation (thermal profiles) alone [14,15]. • Solution to both energy equation and concentration equation (solute concentration profiles) with the assumption of equilibrium solidification conditions; the liquid fraction φ in the two-phase mushy zone [8,9] would be evaluated using the Lever rule and/or Scheil-Gulliver equation [7]. • Solution to both energy and concentration equations (solute concentration profiles) with the assumption of equilibrium solidification condition and evaluate φ from the governing equations without lever rule or Scheil-Gulliver equation assumptions [10,12,16].
All the above-mentioned three categories of approaches could not circumvent the instability caused by the strong non-linearity arising from the phase transformation term in the energy equation. There has been a lack of validation of these models with experiments; validation of transient temperature distribution and time during solidification, because Metals 2022, 12, 608 3 of 32 non-linearity from the phase transformation term impeded the successful solution of the coupled energy and concentration equations (to evaluate the transient liquid fraction in the two-phase mushy zone). Some of the common assumptions to circumvent the non-linearity include linearity in variation of temperature with φ [14,17] and absence of dendrite tip undercooling [7] of the dendrite tip temperature.
Undercooling at the dendrite tip can be caused by three factors: kinetics of solidification, rate of solute diffusion in the liquid and radius of the dendrite tip [7,18]. If the effect of dendrite tip undercooling was considered, the temperature and solute concentration gradients ahead of the dendrite tip would be higher because the temperature at the liquid/two-phase mushy zone interface would be lower than the T liq [19]. Hence, assuming no undercooling will lead to erroneous evaluation of the dendrite arm spacing in the microstructure [20][21][22][23].
In solidification simulation of binary alloy, the dendrite tip undercooling could be included only if φ is determined directly from solving the energy equation. In the present work, a new approach to numerically simulate binary alloy solidification considering the effect of dendrite tip undercooling has been proposed wherein the energy equation alone is used to evaluate φ and a new coupling scheme for T and C L inside the mushy zone has been utilized to avoid instability problems in the energy equation caused by the phase transformation term.
Presently, solidification simulations in the macro-scale (Control Volume, CV~1 mm) do not predict the microstructural features of the solidified alloy such as primary and secondary dendrite arm spacings, λ 1 and λ 2 , respectively [9][10][11][12][13][14]. One of the main reasons for this is the lack of accuracy in predicting the transient temperature distribution and hence, the progressive locations of the mushy zone/liquid interface (solidification time). The algorithm presented below enables predictions of λ 1 because the transient temperature distribution and the solidification times have been accurately predicted and validated with experiments: the reasons being that the new approach included the effect of dendrite tip undercooling prior to the solidification event, ∆T, and further optimized the materials properties as a function of temperature. Traditionally, the undercooling cannot be considered in solidification simulations of binary alloys because the order of magnitude of the diffusion coefficient of the solute in the liquid phase is four orders less than that of the thermal diffusivity and hence, the simulation with ∆T could not be practically completed. In this new approach, the liquid fraction φ in the mushy zone is evaluated by the energy equation with the heat transfer conditions that would enable successful completion of the simulation with the inclusion of the ∆T term. Evaluation of transient G and R would further enable prediction of λ 1 , which would be a valuable tool to predict the solidification microstructure from a macro scale perspective. Further, the new simulation approach has been shown to be faster than the traditional counterpart [10][11][12].

Background
Typically, during dendritic solidification in the mushy zone of a binary alloy, the interface temperature of the dendrite tip and the liquid will require an undercooling (∆T) from the T liq to sustain growth [7]. ∆T depends on G, R, Γ, k, D, m and C o .
Burden et al. [24,25] presented an expression for ∆T as given in Equation (1). The magnitude of ∆T in Equation (1) was validated by experiments for cooling rates ranged up to 3 • C/s [24,25]. In this work, the magnitude of the dendrite tip undercooling is transient for the entirety of the solidification; it is updated after every time step for each control volume (CV).
If numerical simulations do not consider ∆T, the estimation of G will be lower than the actual values in an experiment and this may lead to a prediction of shorter solidification times. The early efforts to simulate binary alloy solidification by Voller et al. [14] considered only the energy equation (enthalpy balance) wherein the instability problem caused by strong non-linearity introduced by the phase transformation term was alleviated by assuming that the temperature in the mushy zone varied linearly with φ. The algorithm was not validated with any experiments. Subsequently, Bennon et al. [9] proposed a numerical algorithm to simulate alloy solidification using one form of the concentration and energy equations in all three zones (solid, liquid and mushy). The volume fraction of liquid was evaluated by applying equilibrium conditions (lever rule) and complete diffusion of solute in the solid was assumed, contrary to reality in a solidification experiment. Moreover, the ∆T value and the volumetric shrinkage (β) during phase transformation were also ignored. No experimental validation of the algorithm was presented.
A numerical algorithm for solidification simulation of binary alloys was proposed by Felicelli et al. [10]; wherein, the liquid fraction φ in the mushy zone was directly evaluated from the energy and concentration equations instead of the lever rule or Scheil-Gulliver paradigm [9]. One form of concentration equation was used to define the solute concentration in the solid, mushy zone and liquid regions alike. One form of energy equation was used in the solid and liquid regions. However, in the mushy zone, to alleviate the instability caused by the phase transformation term in the energy equation, the concentration equation substituted this term to obtain a modified form of the energy equation [9]; φ was evaluated from the concentration equation until the eutectic temperature was attained. Subsequently, at the eutectic temperature, φ has to be evaluated by the energy equation. ∆T was not used in this algorithm because concentration equation was used to determine φ in the bulk of the mushy zone and if ∆T was to be considered at the dendrite tip, then the completion of the solidification simulation would be practically impossible due to the very small order of D (~10 −9 m.s −1 ). An outline of the numerical algorithm used by Felicelli et al. [10] is shown in Figure 1. This algorithm has been one of the most popular methods this date to simulate binary alloy solidification [10,12]. This algorithm did not consider ∆T or β. However, McBride et al. [12] further developed this algorithm by considering β and the fluid flow occurring as a result of volumetric shrinkage. Heinrich et al. [11] corrected the flow fields during solidification at eutectic temperatures. No experimental validation in transient temperature distribution or solidification time of these numerical simulations has been carried out.   The dimensions of the geometry are identical in Figure 2a,b. Equation (4) assumes the local equilibrium condition in each CV of the mushy zone to determine the temperature of the same CV [7].
If re-melting of solid does not occur, the average solute concentration of the solid in each CV was evaluated by Equation (5). Re-melting of the solidified phase occurs in a CV when the convective current with a higher-temperature liquid flows into this CV; if re-melting of solid occurs, C S must be taken from the saved history of C S values during the simulation [10]. In the model, as solid phase inside a CV grows without re-melting, the magnitude of C S has to be recorded with the decreasing of φ from the first piece of solid formed at the beginning of solidification according to Equation (5). When re-melting occurs, C S would have to be obtained from the recorded history according to the relevant value of φ instead of using Equation (5).
The solute mass balance is presented in Equation (6).
The momentum equations in presented in Equations (7) and (8) are in both the r and y directions, respectively.
) using the Boussinesq assumption, which means that all properties in all terms in Equation (8) are assumed constant except for the source term, in which case the density is a function of temperature and concentration.
The transient dendrite tip undercooling was estimated by Equation (1), following the model presented in Burden et al. [24,25]. The permeability of the two-phase mushy zone was defined as a function of the φ by a simplified version of the Kozeny-Carman model [26,27] as proposed by Asai et al. [28] and shown in Equation (9). In Equation (9), the dendrites were assumed to have a conical morphology with a high aspect ratio (length to diameter at the base) and the magnitude of permeability in both the K n (the direction normal to dendrite growth) and K p (direction parallel to dendritic growth) was assumed to be identical.
Equation (10) presents the expression to determine the primary dendrite arm spacing, λ 1 in the solidified structure; this was proposed by Bouchard et al. [20] from a semiempirical analysis of solidification under unsteady conditions. Further, Equation (10) was coupled with the expression to determine λ 1 while considering the effect of fluid flow within the two-phase mushy zone, as proposed by Lehmann et al. [29] and shown in Equation (11).

Domain Definition
Thermophysical properties of the materials considered in this study are provided in Table 1. The initial and the boundary conditions of the problem are presented for two different solidification conditions as follows: Upward Solidification → heat extraction and gravity vector are in the same direction. Figure 2a presents the schematic and boundary conditions for C L and T. There is no natural convection of the liquid because the density gradient is positive in the direction of the gravity vector [31]. The heat transfer coefficient (h) at heat extracting side for the Al-5wt%Si [18]; while the same for the Al-3wt%Si alloy is presented in Equation (12), which is an empirical correlation of the transient temperature ( • C) as a function of time (s), extracted from the experiments by Peres et al. [18]: Downward Solidification → heat extraction is in the direction opposite to the gravity vector. Figure 2b presents the boundary conditions for C L and T. There is natural convection of the liquid phase because the density gradient is positive in the direction opposite to the gravity vector [31]. The heat transfer coefficient at heat extracting side for the Al-5wt%Si was h = 2400 (W m −2 K −1 ) [30]; while the same for the Al-3wt%Si alloy is presented in Equation (13), which is an empirical correlation of the transient temperature ( • C) as a function of time (s), extracted from the experiments by Spinelli et al. [30]: Equation (14) defines the concentration at the heat extraction boundary (y = 0) for both the upward and downward solidification.
In Figure 2, the initial condition of the liquid alloy in the domain has a marginal superheat temperature over the T liq . All dimensions and boundary conditions in the simulation were the same as used in the experimental study [18,30]. Heat flux → q in heat extracting boundary would be progressively smaller during solidification according to experimental conditions [18,30] for both upward and downward cases. The simulation geometry is a cylinder with a diameter of 55 mm and height of 145 mm. As shown in Metals 2022, 12, 608 8 of 32 Figure 2, the computational domain is a two-dimensional axi-symmetric cross-section of the geometry with cylinder of 55 mm diameter and 145 mm height.
The velocities at the boundaries of the computing domain are presented in Figure 3, wherein, the velocity, u y (r) at y = 145 mm is a parabolic function of r [32] and evaluated by Equation (15): In Equation (15), the maximum velocity at r = 0, u max is evaluated from the solidification shrinkage. Volume flow rate (VFR) due to solidification shrinkage at y = 145 mm was calculated at each time step by Equation (16): where, R i is the radius and H is the height of the domain. The boundary conditions at each time step is determined by evaluating u max , which is obtained from the updated values of VFR at each time step. The initial condition for velocity at t = 0 for both the domains in Figure 3 is u r = u y = 0. The initial liquid temperature for upward solidification domain in Figure 3a is 2 • C superheat above the respective liquidus temperatures for the alloys (644.0 • C for Al-3wt%Si, and 617.3 • C for Al-7wt%Si) and that for the downward solidification domain in Figure 3b is 5 • C superheat above the respective liquidus temperatures for the alloys (633.7 • C for Al-5wt%Si and 620.3 • C for Al-7wt%Si). Initial conditions for C S , λ 1 and φ are: C S = kC o , λ 1 = 300 µm and φ = 1, respectively.
The initial conditions for the solute concentrations for the upward solidification were 3 wt% Si and 7 wt% Si and those for downward solidification simulations were 5 wt% Si and 7 wt% Si. The initial conditions for the temperature and solute concentration were identical to those reported in the experimental work [18,30]. The material properties were obtained from two sources [2,4] and presented in the earlier nomenclature section of this publication.

Numerical Simulation Procedure
One of the novelties of the numerical algorithm proposed in this study is the incorporation of ∆T during solidification in our simulations. Before the details of the numerical Metals 2022, 12, 608 9 of 32 simulation procedure are presented, an overview of the methods developed to incorporate ∆T in the simulation is presented below. Figure 4 presents a schematic of the steps during uni-directional solidification of the liquid from an initial melt superheat temperature of T ini to the final eutectic temperature of T eut .  Table 2 shows the steps involved during solidification of a binary hypoeutectic alloy. In this study, all these steps were considered in each CV during solidification simulation. Table 2. The various steps in the solidification of a CV of a binary alloy.

Steps
Time T C S (New Layer) C L Figure Step 1 0 to t 1 T ini to T liq 0 C o Figure 4b,c Step 2 Step Step Step 5 Step 6 t n T eut C eut C eut Figure 4h Step 7 t final T eut C eut C eut Figure 4i A detailed description of each step is also given below: Step 1. T ini decrease to T liq (Figure 4b,c) mainly due to the heat extraction during solidification. Concentration of the alloy melt is maintained at C o . The time taken for this step is t 1 .
Step 2. No solidification event occurs at T liq . The first solid appears after an undercooling of ∆T as shown in Figure 4d. The extent of undercooling is evaluated by Equation (1). The concentration of the first solid phase is represented by kC L1 wherein the liquid concentration C L1 is marginally greater than the average alloy composition C o . The time taken for this step is t 2 .
Step 3. During the next time step (t 2 + ∆t), the temperature remains constant as that in Step 2 (T liq − ∆T). Evolution of the primary phase occurs during this step and the temperature will continue to remain constant (latent heat of fusion balances the heat extraction) until the solute concentration in the liquid at the mushy zone/liquid interface reaches C * L as shown in Figure 4a,e. Step 4. The solute concentration of the liquid at the mushy zone/liquid interface reaches C * L in time t* and the solute concentration in the solid is given by kC * L as shown in Figure 4a,f.
Step 5. The temperature begins to decrease again and layers of solid of the primary phase are formed at composition and the solute re-distribution in the liquid phases is defined by the concentration equation given by Equation (3), as shown in Figure 4g. The solute concentration in the all layers of solid of primary phase is evaluated by Equation (5) for each time step if there is no re-melting of the solid phase due to convection of the liquid. In case re-melting occurs, the history of the solute concentration in the solid is used to evaluate the average concentration of the solid inside CV [10,11].
Step 6. Temperature of the liquid reaches the eutectic temperature and a layer of solid with eutectic solute concentration is formed marking the end of growth of the primary phase, as shown in Figure 4h.
Step 7. The remaining liquid at the eutectic temperature is transformed into a twophase mixture to complete the solidification of the CV, as shown in Figure 4i.
The unique feature in this model is the Steps 2 to 4 where the temperature of the CV is maintained at T liq − ∆T until the solute concentration in the liquid inside CV reaches the equilibrium concentration of the solute, C * L , as suggested in Figure 4a. This is observed during the solidification process of all binary alloys; however, the event has never been accounted for in any solidification simulation of binary alloys until this study. Introducing Steps 2 to 4 will identify and account for the time period of thermal recalescence during solidification. Table 3 presents the identification and details of the simulations in this study. Table 3. Notations and details of the solidification simulations.

Simulation Identification
Type Equation (1) Typical deliverables from the simulations listed in Table 3 include transient temperature distribution and solidification time for the mushy zone/liquid interface, as well as other quantified parameters such as λ 1 , G and R. In the governing equations (Equations (2), (3) and (6)-(8)), the central differencing scheme was employed for all the diffusion terms; the convection terms were discretized by hybrid difference scheme. The implicit method was used for the transient terms. Convergent criteria for all parameters were determined by the two orders of mesh size, since the numerical simulation in the present study is in two orders of accuracy. Figure 5 presents the process flow chart for the numerical simulations: there are four procedural loops, L1, L2, L3 and L4. At every time step, Equations (2)-(5) are processed in L1 to evaluate the transient temperature distribution, solute concentration distributions in both the liquid and solid phases, and φ in the mushy zone. Equations (6)- (8) are processed in L2 to obtain the velocity distribution of liquid (fluid flow) using the algorithm defined by the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) [33]. Momentum equations (i.e., Equations (7) and (8)) are solved using an initial pressure field. The pressure field is corrected using the mass balance equation (Equation (6)). By iterating these two major steps, a convergent solution for velocity and pressure fields can be obtained. Each time step is then advanced with L3 and ∆T and λ 1 are evaluated using the loop L4. The loop L4 enables the evaluation of λ 1 while incorporating the undercooling, ∆T at the macro-level (CV1 mm) solidification simulation. Previous researchers could not introduce loop L4 to predict microstructure features such as λ 1 because of the lack of valid correlations of the transient temperature distribution and solidification time with experiments. Prior to incorporating the loop L4, the optimization of critical parameters such as thermo-physical properties as a function of temperature, shrinkage, fluid flow from natural convection and undercooling were carried out with several preliminary simulations to enable experimental validation of the transient temperature distribution and solidification times during solidification. The main feature of the algorithm presented in Figure 1 is the procedural logic used in loop L1 shown in Figure 5: • At time step t = t m , Equation (3) was used to evaluate the distribution of solute concentration in all of liquid, solid and the two-phase mushy zone. • Solving Equation (2) and comparing the evaluated T of each CV to (T liq − ∆T), the CV was determined to be in either the liquid phase or mushy zone. • Equation 4 was used to evaluate temperature distribution inside the mushy zone or assumed as (T liq − ∆T) if C L < C * L (Figure 4d-f). Equation (2) was used to evaluate the temperature distribution inside the liquid and solid phases.

•
The temperature distribution, evaluated from Equation (2) was used to evaluate φ in the mushy zone.

•
If there is no re-melting of the solid in any given CV, Equation (5) was used to calculate the mean solute concentration of the solid phase in the CV, and in case of solid phase re-melting in any given CV, the history of prior mean solute concentration in the solid phase was used for the current time step [10]. • All the above steps were repeated until convergence.
In any given CV in the domain, the temperature of the CV to determine if it is presented in either the entirely liquid or solid phase regimes of the domain was evaluated by the energy equation, Equation (2); meanwhile, the temperature in the two-phase mushy zone was evaluated by the Equation (4) with local equilibrium assumption. To determine, if the given CV falls within the mushy zone, the temperature of the liquid in that CV was evaluated using Equation (2) and if this temperature value is less than T liq , then the temperature within that CV is re-evaluated using Equation (4). Further, this temperature of the liquid in a CV in the mushy zone was input back in the Equation (2) to calculate φ in the mushy zone. This procedure in the first iterative process does not include the undercooling (∆T = 0) and hence, the first iteration of loop L4 would have been carried out with ∆T = 0; meanwhile, subsequent iterations of L4 will yield a positive value for ∆T and updated values of the temperature distribution in the CV.
The novelty of the algorithm presented in this study, wherein the energy and solute concentration equations were coupled, is in the loop L1. The advantages of loop L1 in this new algorithm ( Figure 5) over that by Felicilli et al. [10] (Figure 1) that enabled the incorporation of the undercooling term, ∆T are in two aspects: firstly, the former solves the solute concentration equation first for all the three regimes of liquid, solid and mushy zone to circumvent the instabilities caused by the phase transformation term in the energy equation, Equation (2), and secondly, the new algorithm uses the heat transfer condition, solely, to evaluate the φ in the mushy zone. As a result, a significant reduction in the computational time was realized from using this newly proposed algorithm as well.
In summary, the advantages of the new algorithm for solidification simulation of binary alloys over the existing ones are the ability to incorporate the undercooling at the liquid/mushy zone interface, ∆T, overcome the instabilities caused by the phase transformation term in the energy equation due to strong non-linearity and significant reduction in computation cost. Further, a successful and complete correlation of the simulations by experiment results was also carried out to ensure its validity.

Results and Discussion
All the results and discussion are in reference to Figures 2 and 5. The results and discussion are presented for the following topics:

Comparison of New Algorithm with the Popular Version
The results at the end of the first iteration of L4 in Figure 5 are presented as simulation SA and were used to compare with the results from the popular algorithm (simulation SB) [10][11][12] to determine the reliability of this newly proposed algorithm (shown in Figures 1 and 5).
Comparison of the transient temperature distribution for SA obtained in the present study at the end of the first iteration of L4 in Figure 5 and that SB simulated by the algorithm presented in Heinrich et al. [11] (Figure 1) is reported in Figure 6.  Figure 6 shows the temperature distribution at r = 0 during solidification without the effect of ∆T for seven locations between y = 0 and y = 84 mm, respectively; also shown is the comparison between the present work and that of Heinrich et al. [11] showing that using the new algorithm proposed in this study is as accurate as the popular method which indicate that the alternative logic used to determine φ from solely the heat transfer conditions rather than the concentration equation is valid. Further, the computing time for the new algorithm to obtain the results shown in Figure 6 was reduced by a maximum of 50% depending on the initial solute concentration of the alloy and boundary conditions. The average cooling rates investigated in Figure 6 varied between 0.72 • C/s to 2.05 • C/s.

Transient Temperature Distribution
Transient temperature distribution is critical to evaluate solidification time, G, R and λ 1 . The transient temperature distribution in the computing domain was validated with experiment results presented by two groups of researchers [18,30]. Figure 7a,b show the validation of the present algorithm using the experiments by Peres et al. [18] and Spinelli et al. [30]. Figure 7a,b are for simulations S2 and S8, respectively. Simulations S2 and S8 show a good agreement with experiments and validate the numerical method. In Figure 7a, the maximum deviation in temperature is less than 10 • C located at y = 64 mm and t ≈ 100 s for upward solidification. When compared with the initial melt temperature of 644.05 • C, the deviation is less than 2%. In Figure 7b, the maximum deviation is less than 13 • C (2%) located at y = 64 mm and t ≈ 100 s for downward solidification.

Effect of Undercooling on Solidification Time
Peres et al. [18] have carried out experiments on upward solidification (Figure 2a) with Al-3wt%Si and Al-7wt%Si alloys, respectively. They have reported the solidification time (location of the mushy zone/liquid interface) by estimating the time taken for the interface to reach the location of each thermocouple in the experiments. These experiments were used to validate simulations S2 and S4. Figure 8 shows the results of these validations.  Figure 8a shows the validation of the solidification times from simulation S2 (with ∆T) with the experiments [18] along with results from S1 (no ∆T). Figure 8b shows the validation of the solidification time from simulation S4 (with ∆T) with the experiments [18] along with results from S3 (no ∆T). Figure 8 shows that the simulations that include the effect of ∆T on solidification has an excellent agreement with the experiments with a maximum time deviation of <1.5% for the farthest thermocouple (y = 84 mm) in both Figure 8a,b (refer to Figure 2a). The deviation may be due to the lack of uni-directional heat transfer in the experiments towards the later stages of solidification. However, the negligible deviation shows that experiments were nearly uni-directional for most of the solidification time. It can be seen in Figure 8a,b that simulations S1 and S3 with no consideration of ∆T during solidification have a large deviation from the respective experiments. The maximum time deviation in S1 and S3 from experiments is 11.4% and 19.8%, respectively, at y = 84 mm. In Al-3wt%Si alloy, the total time taken for the mushy zone/liquid interface to travel between y = 0 to y = 84 mm is 88 s from experiments, 86.7 s for S2; and 78 s for S1. In Al-7wt%Si alloy, the total time taken for the mushy zone/liquid interface to travel between y = 0 to y = 84 mm is 96 s from experiments, 94.7 s for S4; and 77 s for S3. The consideration of ∆T in the simulation reduces the velocity of the liquid/mushy zone interface, R because it requires longer time for each CV to reach an undercooled temperature state (T liq − ∆T) before any solidification event.
Spinelli et al. [30] have carried out experiments on downward solidification with Al-5wt%Si and Al-7wt%Si alloys. The results from these experiments were used to validate simulations S6 and S8. Figure 9a,b shows the validation of simulation S6 and S8. Additionally shown in Figure 9a,b are the solidification times for simulations S5 and S7. Figure 9 shows that simulations S6 and S8 (with ∆T) have better agreement with the experiment results. The maximum deviation of S6 and S8 from experiments at y = 85 mm is 3.6% and 16%, respectively. The larger deviations from the experiments observed for the downward solidification in Figure 9 as compared with the upward solidification presented in Figure 8 may be due to the fact that in experiments [18,30], the temperature gradient between the side wall of the container (right side of computing domain) and the liquid alloy inside increases with time in the experimental condition. In downward solidification, the fluid flow is more pronounced due to the natural convective currents induced by this temperature gradient in experiments, which results in a more significantly increased temperature gradient in the liquid side of mushy zone/liquid interface causing a smaller value of R and longer solidification time in downward solidification in experiments. Another reason for the marginal under prediction of solidification times in Figures 8 and 9 may be due to the added latent heat of fusion of the eutectic phases which will result in an over prediction of R and an under-prediction of the solidification times. In Al-5wt%Si alloy, the total solidification times for experiment, S6 and S5 are 113 s, 109 s and 97 s, respectively. In Al-7wt%Si alloy, the total solidification times for experiment, S8 and S7 are 130 s, 108 s and 96 s, respectively.    Figure 10d-f present that for simulation S1. In these two simulations (upward solidification) the positive density gradient is in the direction of the gravity vector (refer to Figure 3a) and hence, there will be no natural convection of the liquid which is evident from fluid flow from top to bottom and the absence of circular flow. Since only shrinkage affects the fluid flow velocities and shrinkage is almost identical in both simulations, S2 and S1, there is a marginal increase in the velocity of S1 (without ∆T) as compared to that in S2 (with ∆T) due to the larger value of R in S1. A similar observation was made when simulation S4 was compared with S3 as well.

Effect of Undercooling on Fluid Flow-Downward Solidification
The development of fluid flow with the effect of undercooling for simulation S6 is shown in Figure 11.  Figure 11f shows typical domains found in the simulation results: solid, mushy zone and liquid. The solid lines inside the mushy zone in Figure 11 represent the constant φ lines from φ = 0.9 to φ = 0.1.  Table 4 presents a summary of the results shown in Figure 11. Table 4 suggests that effect of shrinkage is dominant at the beginning of solidification due to the lack of natural convection. The effect of natural convection gradually begins to dominate and the effect of shrinkage on fluid flow velocity becomes negligible as solidification time increases. At the beginning of solidification, R is high and induces an instantaneous flow field in the liquid which is primarily caused by the solidification shrinkage. However, the flow due to the density gradient in the liquid develops from an initial rest to subsequently dominate the fluid flow regime over the solidification shrinkage. Higher velocities of fluid flow caused by natural convection will result in greater circular flow regions, as seen in Figure 11. There are two mechanisms that affect strength of fluid flow due to natural convection. One of them is the density gradient caused by the temperature and solute concentration gradients in the liquid which is directly proportional to the fluid velocity and the other is the height of the liquid (space available for fluid flow) which is inversely proportional to fluid velocity. Between times t = 5 s and 75 s, the effect of the density gradient greatly increases the fluid flow velocity due to natural convection. However, between t = 75 s and 110 s the effect of the reducing height of the liquid marginally decreases the fluid flow velocity due to natural convection. The maximum fluid flow velocity is always in the region close to r = 0 and is in the direction of heat extraction. Figure 11 and Table 4 also show that the length of the mushy zone is continuously increasing indicating an unsteady solidification condition at all times because of the continuously decreasing heat flux → q at y = 0 mm. The maximum flow in the computing domain can always be observed at the left top of the liquid region from Figure 11b,f (center of the solidifying cylinder). Fluid flow with a higher temperature from liquid region is driven into mushy zone in this region. Therefore, solidification time will increase in this region and hence the constant liquid fraction lines are compressed along the y direction at the liquid/mushy zone interface (Figure 11b-f). The direction of the fluid flow velocity at the liquid/mushy zone interface is predominantly in the positive r direction due to the development of clockwise fluid flow: resulting in the accumulation of enthalpy on the top ride side of the liquid region (along the wall of the solidifying domain). As a result, solidification in this region is slowed down and the constant liquid fraction lines along the wall are also compressed as shown in Figure 11d-f.
The normalized iso-stream function lines with the effect of undercooling in simulation S6 are shown in Figure 12. Figure 12a-f shows snapshots of the normalized iso-streamlines in the computing domain at increasing solidification times of 5 s, 25 s, 50 s, 75 s, 100 s and 110 s, respectively. Maximum magnitude of stream function is shown below each snapshot for the respective times. Similar to Figure 11, the top solid lines inside the mushy zone in Figure 12 represent the constant liquid fraction lines from φ = 0.9 to φ = 0.1. The strength of fluid flow can also be implied by the maximum magnitude of stream function. Maximum magnitude of the stream function in Figure 12 changes from a lower value (Figure 12a) to a maximum value (Figure 12d) and back again to a lower value (Figure 12f), which indicates the strength of flow develops from weaker to stronger and back to weaker as the solidification proceeds. Arrows in the streamlines in Figure 12 illustrate the flow direction, which indicates that flow develops from bottom to top at t = 5 s. Circular flow developed after t = 25 s implies that the flow is dominated by natural convection. Clockwise cell can be observed for times between t = 25 and 75 s. After t = 100 s, a second flow cell is created at the right bottom side of the computing domain due to the higher velocity at t = 75 s during the development of the fluid flow. A notable feature is the development and strengthening of a secondary circular flow starting from t = 75 s as shown in Figures 11 and 12.   Figure 13a,b for velocities at t = 5 s, the maximum magnitude of u y is more than 10 times that of u r , which indicates that the flow is dominated by shrinkage.
Positive direction of u y located at the region close to the right side of computing domain (wall of solidifying cylinder) implies that natural convection has already affected fluid flow even if it has not dominated the flow pattern at t = 5 s. The flow in the middle of the solidifying cylinder in the y direction is nearly fully developed as indicated by the lack of variation in the magnitude of u y and u r at locations (H-y) = 0.04 m, 0.06 m and 0.08 m in Figure 13a,b, respectively. The higher variation in the magnitudes of velocities in Figure 13a,b at locations (H-y) = 0.01 and 0.113 m is due to the effect of the cavity end walls. In Figure 13c,d, the lack of significant variation between the maximum magnitudes of u y and u r at time t = 50 s show that the flow is circular and caused predominantly by natural convection. At t = 50 s, the maximum velocity is located near the liquid/mushy zone interface at a location close to the left side of the computing domain. Figure 14 shows the development of fluid flow of simulation S5 (no effect of ∆T on solidification). Comparing Figure 11 with Figure 14 shows that the nature of development of fluid flow dominated by shrinkage initially and subsequently dominated by natural convection are similar. However, the only difference is in the magnitude of the fluid flow velocities at the various time steps. Table 5 shows a comparison of the magnitude of fluid flow velocity between simulations S6 and S5.   Due to the absence of ∆T in S5, the temperature gradient in the liquid ahead of the mushy zone/liquid interface is always less than that in S6. Hence, the value of R in S5 is always greater than that in S6. In the beginning of solidification, the higher value of R will cause a larger effect of shrinkage on the fluid flow velocity as reflected by the higher magnitude of velocity at t = 5 s as shown in Figure 14 and Table 5. However, due to the lower temperature gradient in S5, the solute density gradient in the liquid is lower than that in S6 at t ≥ 25 s, which will result in a lower fluid flow velocity due to natural convection in the liquid in S5 as shown in Figure 14 and Table 5. This underestimation of fluid velocity due to natural convection in S5 as compared to S6 will greatly affect the evaluation of λ 1 as will be discussed in subsequent sections of this manuscript.
The normalized iso-stream function lines without the effect of ∆T for simulation S5 are shown in Figure 15. Figure 15a-e shows snapshots of the normalized iso-streamlines in the computing domain at increasing solidification times of 5 s, 25 s, 50 s, 75 s, 100 s, respectively. Maximum magnitude of stream function is shown below each snapshot. The difference between the flow field shown in Figure 15 (without effect of ∆T) and that shown in Figure 12 (with effect of ∆T) is best indicated by the difference in the maximum magnitude of the respective stream functions. As shown in Figures 12a and 15a, a higher magnitude of stream function for simulation S5 implies stronger fluid flow at t = 5 s. After the initial stage of solidification, results from simulation S6 always have higher magnitudes of maximum stream function due to the stronger buoyancy flow caused by the higher-density gradient in the liquid region with ∆T. The weaker fluid flow due to the negligence of the effect of ∆T in solidification (Figure 15) results in a delay in the development of the secondary circular flow cell in the liquid when compared to the simulation S6 with the effect of ∆T (Figure 12).   Figure 16 is similar to the results shown in Figure 13. It is found that higher magnitudes of velocity could be obtained in simulation S5 at t = 5 s. However, stronger fluid flow can be observed at t = 50 s in simulation S6.   Figure 17e-h show that for simulation S7. The comparison of velocity distributions between S8 and S7 is similar to that observed in S6 and S5. With the same line of reasoning as presented for S6 and S5, it can be observed that when ∆T is not considered during solidification (S7) the maximum velocity at the beginning of solidification is almost the same as that with ∆T (S8) (compare Figure 17a,e) because fluid flow induced by shrinkage dominates over that induced by natural convection. However, at larger time intervals, the velocities of S7 are lower than that of S8 because the effect of natural convection is more on the fluid flow velocities due to ∆T (compare Figure 17b-d with Figure 17f-h, respectively). ure 17e-h show that for simulation S7. The comparison of velocity distributions between S8 and S7 is similar to that observed in S6 and S5. With the same line of reasoning as presented for S6 and S5, it can be observed that when ΔT is not considered during solidification (S7) the maximum velocity at the beginning of solidification is almost the same as that with ΔT (S8) (compare Figure 17a,e) because fluid flow induced by shrinkage dominates over that induced by natural convection. However, at larger time intervals, the velocities of S7 are lower than that of S8 because the effect of natural convection is more on the fluid flow velocities due to ΔT (compare Figure 17b-d with Figure 17f-h, respectively).  Figure 18 shows the distribution of the temperature gradient, G in the liquid ahead of the mushy zone/liquid interface. There are two main reasons for the decrease in the value of G with increasing solidification time: one is the decrease in the heat flux at the heat extraction boundary with time and the other is the decrease in the bulk temperature of liquid with time. Figure 18a,b shows the results for simulations Al-3wt%Si (S1 and S2) and Al-7wt%Si (S3 and S4) alloys, respectively; and Figure 18c,d are for Al-5wt%Si (S5 and S6) and Al-7wt%Si (S7 and S8) alloys, respectively. Since the temperature of the mushy zone/liquid interface will be lower when we consider the undercooling, ∆T, the value of G is higher for the solidification simulations with the effect of ∆T included (S2, S4, S6 and S8) than that in those without the effect of ∆T (S1, S3, S5 and S7) for the respective alloys. In Figure 18a,b (upward solidification), G decreases during the solidification process. In Figure 18c,d (downward solidification), there seem to be three distinct regions from the beginning to end of solidification. These regions are marked by line segments AB, BC and CD in all the curves in Figure 18c,d. Table 6 presents the four locations A, B, C and D of the mushy zone/liquid interface, solidification times, maximum fluid flow velocity, variation of fluid flow velocity from previous location and the temperature gradient, G in the liquid at the mushy zone/liquid interface at each location.  As discussed in the previous section on the effect of ∆T on fluid flow during downward solidification (refer to   Table 6), the fluid flow in the initial stages of solidification as defined by line segment AB is dominated by shrinkage. Between points A and B, the decrease in G is only due to the decrease in the heat extraction rate between these locations. Between points B and C, the strong positive density gradient develops against the direction of the gravity vector and causes a significant increase in the fluid flow velocity due to natural convection.

Effect of Undercooling on G
In simulation S5 (no ∆T) the change in the velocity between B and C is 557% and that in S6 (with ∆T) is much higher at 1109% because the consideration of ∆T will enhance the fluid flow velocities due to larger density gradients. This increase in velocity between B and C in both S5 and S6 will cause a tendency to increase the magnitude of the gradient of G (Figure 18c and Table 6). In S6, the velocity change is so drastic due to the added effect of ∆T and there is an increase in the value of G between locations B and C. The same argument can be made for the change in G between locations B and C for simulations S7 and S8 as well (Figure 18d and Table 6). Figure 19 shows the distribution of the solidification velocity, R at r = 0 and various locations of y = 0 to 84 mm in the computing domain ( Figure 3). The value of R was estimated by the velocity of the mushy zone/liquid interface during solidification. Figure 19 shows the value of the velocity of the mushy zone/liquid interface (φ~0.999) velocity, R. Figure 19a,d shows the results for solidification simulation for Al-3wt% (upward), Al-7wt%Si (upward), Al-5wt%Si (downward) and Al-7wt%Si (downward), respectively. Figure 19 shows that for all the simulations transient value of R decreases overall with increasing distance from the heat extraction interface because the magnitude of heat flux at the heat extraction boundary and the computing domain decreases with time as the volume of the solid phase increases. In all results shown in Figure 19, the value of R evaluated with the effect of ∆T is always lower than that without the effect of ∆T because the introduction of ∆T will require a lower temperature for the solidification than the T liq , which will cause a delay in the solidification event. The value of R is high at the beginning of solidification for all the simulations in Figure 19. The effect of undercooling on R is more pronounced at the earlier stages of solidification. In Figure 19a,d, the maximum deviation in the estimated value of R at 8 mm from the heat extraction interface is 18%, 89%, 33% and 26%, respectively. In Figure 19c,d, the value of R for simulations S5 and S7 (no ∆T) increase towards the end of solidification of the computing domain because, without the effect of ∆T for downward solidification, the value of G will decrease more significantly towards the end of solidification and result in an increase in R (between y = 58 mm to y = 84 mm). For upward solidification in simulations S1 and S3 (Figure 19a,b), the value of G does not decrease enough to cause an increase in R. If one considers the effect of ∆T during downward solidification as in simulations S6 and S8 in Figure 19c,d, the value of G will be greater than those without the effect of ∆T, as shown in Figure 18c,d, and will not decrease so low to as to increase the value of R. Figure 20a,b show the distribution of the primary dendrite arm spacing, λ 1 as a function of cooling rate (G.R) of the liquid at the mushy zone/liquid interface. It can be observed that there is no significant change in the value of (G.R) between the simulations with the effect of ∆T (simulations S2 and S4) and that without (S1 and S3), respectively, because the effect of ∆T will increase the value of G but simultaneously decrease the value of R and hence the product (G.R) will remain fairly unchanged. Additionally shown in Figure 20a,b are the values of λ 1 from the experiments in upward solidification from Peres et al. [18]. It can be seen that simulations S1 to S4 predict the respective values of λ 1 with a reasonable accuracy. There is a slight under-prediction of the value of λ 1 because of a combination of two reasons: the model (Equations (1) and (2)) itself is semi-empirical for an unsteady state solidification process and hence may tend to over-predict or the cooling rates of the mushy zone/liquid interfaces were measured after the liquidus temperature was recorded by the thermocouples in the experiment [18] and the magnitude of this value will be higher than the actual cooling rate at the mushy zone/liquid interface and will result in an over-prediction of (G.R) in the experiments, as shown in Figure 20a,b. If the value of (G.R) is being over-predicted by the experiment results, then the value of λ 1 should be under-predicted (Equation (1)) as it is shown as function of location (Figure 21a) for simulations S1 and S2. In Figure 21a, the value of (G.R) at any specific location (Figure 3) will be higher when the effect of ∆T is considered than that without the effect. Hence, the value of λ 1 obtained at any specific location will be lower in simulations with the effect of ∆T than in those without the effect. Hence in Figure 21a, the simulation S2 (with ∆T) predicts the value of λ 1 more accurately than that by S1 (without ∆T).  [18]; (b) 7% upward compared to experiments [18]; (c) 5% downward compared to experiments [30]; (d) 7% downward compared to experiments [30].

Figure 21
. λ 1 distribution at r = 0 mm and y ranged from 8 mm to 84 mm for simulations (a) S1 and S2 compared to experiments [18]; (b) S5 and S6 compared to experiments [30]. Figure 20c,d shows the distribution of λ 1 as a function of (G.R) for the cases of downward solidification in simulations S5 to S8. As discussed previously about the effect of ∆T on G in the downward solidification, there are three distinct regions wherein the magnitude of the value of gradient of G decreases at the first region and increases at the second region (fluid flow dominated by the emergence of density gradient causing significant increase in fluid flow velocities than that caused by shrinkage in the first region) and again decreases at the third regions. The trend of λ 1 will also follow the trend of G in these three regions as in Figure 20c,d. The value of λ 1 increases (G decreases) in the first region between points A and B, then λ 1 decreases when the magnitude of gradient of G increases in the second region between points B and C and subsequently when the fluid flow velocities due to density gradient is stabilized in the third region, the value of λ 1 again increases between points C and D in both simulations with and without the effect of ∆T. However, the effect of undercooling, ∆T, on the magnitude of the fluid flow velocities in the three regions will have a greater effect on the changes in the value of G in the three regions and hence a greater effect on the value of λ 1 in these regions as well (refer to Equation (11)). Figure 20c,d also shows that the λ 1 estimated in the experiments by Spinelli et al. [30] seems to be closer to the values predicted by the simulations with the effect of ∆T included.
The values of λ 1 in the first region in alloy Al-5wt%Si (simulations S5 and S6) between points A and B seem to be under-estimated by the experiments and this may be due to various errors in the experiment evaluation of the cooling rate at the mushy zone/liquid interface and λ 1 caused by the high value of R in this region. Further, the decrease in λ 1 between points B and C can be observed from the experiment results as predicted by the simulations in Figure 20c [30]. Experiments will have to be carried out to obtain more data points in this region to validate the simulation results. Figure 21a,b show the distribution of the value of λ 1 as a function of location (refer to Figure 3) for Al-3wt%Si upward and Al-5wt%Si downward solidification, respectively. These show a better prediction of the λ 1 by the simulation with the effect of ∆T than the one without. The anomalous behavior of the trend in λ 1 can also be observed in the second region between points B and C in Figure 21b. It is due to the development of the higher fluid flow velocities caused by the emergence of the density gradients causing natural convection in the computing domain.

Summary and Conclusions
A new and more efficient algorithm to simulate the solidification of binary alloys that includes the effect of undercooling at the dendrite tip during solidification has been proposed and validated. The present algorithm circumvented the problem of instabilities of the energy equation caused by the phase transformation term due to the strong non-linearity by solving the solute concentration equation prior to the energy equation. Further, the liquid fraction in the two-phase mushy zone during solidification was evaluated entirely from the energy equation without the use of the solute concentration equation. The transient temperature distribution in liquid, mushy zone and solid regions during solidification was validated by experiments. The solidification time, which indicated the instantaneous location of the mushy zone/liquid interface, was evaluated and validated with experiments as well. This algorithm was further used to evaluate solidification parameters, such as fluid flow, G, R and λ 1 of Al-Si hypoeutectic alloys in both upward and downward solidification modes.
The main conclusions of this study are presented as follows: • The undercooling has been successfully quantified in numerical simulations and has a significant impact on the selection of the primary arm spacing during solidification, especially in modes that involve development of significant fluid flow during the process arising from notable density gradients in liquid from solute redistribution. The effect of undercooling, ∆T has a significant effect on the velocity of the fluid flow in the computing domain only in the downward solidification and there is no significant effect in the upward solidification mode. There are three distinct regions in the downward solidification mode because of significant changes in the flow velocity arising from gravity-assisted density gradients in the liquid in these regions.

•
Simulations with the effect of undercooling were found to increase solidification time and be far more agreeable with the experiment results for solidification times than those without the inclusion of the undercooling term.

•
The effect of ∆T is pronounced on the distribution of the temperature gradient of the liquid at mushy zone/liquid interface, G, during both upward and downward solidification modes.

•
The effect of ∆T is pronounced on the distribution of the velocity of the mushy zone/liquid interface, R, during both upward and downward solidification modes. The effect of ∆T is more pronounced at higher values of R and reduces with decreasing R because ∆T is high at higher values of R as a result of increased heat extraction during initial stages of solidification.
This manuscript describes the impact of undercooling on flow behavior and temperature distribution and its effect on the primary dendrite arm spacing. Results of this macro-scale simulations model provides valid input such as boundary and initial conditions for thermal, solute and fluid flow regimes for the subsequent micro-scale models that further evaluate the morphology and details of microstructure developments during binary alloy solidification [34]. It is recommended in the future to investigate the impact of including the dendrite tip undercooling on the secondary dendrite arm spacings and the local solid fraction development in the mushy zone. The secondary dendrite arm spacing selection would also be affected by the presence of the dendrite tip undercooling [35]; however, the magnitude of fluid flow regimes and density gradients arising from solute redistribution within the mushy zone regions of the primary arms for the development of secondary dendrite arms would be notably different (less) than those developed ahead of the primary dendrite tip.
This numerical model is anticipated to generate a valid simulation for solidification of binary metallic alloys with considerations of dendrite tip undercooling; the caveat being that the constitutive expressions from the critical solidification parameters such as undercooling and thermophysical properties of the system are known.