Prediction of Limiting Casting Speed in a Horizontal Direct-Chill Casting through Numerical Modeling and Simulation

: Limiting casting expression speed was obtained and the ﬂ ow redistribution and thermal history in a horizontal direct-chill (HDC) casting was predicted using the numerical modeling approach. The governing solidi ﬁ cation equations were non-dimensionalized to understand the relevant contribution of each term in the solidi ﬁ cation processes in the HDC system. The e ﬀ ect of an increase in the casting speed on the ﬂ ow characteristics and sump length was represented by the Péclet number Pe. Details of the simulation reveal that at a low Pe, the natural convective ﬂ ow creates minor counter-clockwise recirculating cells in the lower half of the HDC domain. However, at a Pe above 82.75, the minor recirculating cells disappear due to the strong forced convective ﬂ ow from the upstream. Additionally, an increase in the Pe increases the sump length, strength, and spread of the turbulence ﬁ eld within and beyond the mold region. The limiting casting conditions are computed by predicting the sump length over which the alloy temperature is above the solidus temperature. This gives a simple relation for the casting speed as a function of the geometrical data and the alloy properties. The current work is useful to casting engineers who always rely on trial and error in choosing a new casting speed whenever a new alloy is to be produced. Hence, with the new information and the casting speed relations, it is possible and easy to predict the operating window over which melt break-out can occur during HDC.


Introduction
Among the several casting processes used in processing aluminum alloys, direct-chill (DC) casting is a major one. The DC casting permits the hot liquid metal to moderately solidify in a water-cooled hollow mold [1,2]. Successive freezing (cooling) is by direct water quenching at the semi-solidified billet surface. The freezing in the mold region is called primary cooling, while the region where water directly contacts the semi-solidified shell is called secondary cooling [1,2]. The maximum heat removal rate is in the secondary region.
DC technology can operate in a vertical mode referred to as vertical direct-chill (VDC) or in a horizontal mode refer to as horizontal direct-chill (HDC) [3]. In the VDC casting arrangement, the semi-solidified ingot is withdrawn vertically downward by a hydraulic ram (Figure 1). The casting speed for the VDC is lower than the HDC casting speed. This is because in the VDC system gravity acts downward in the casting direction and aids the withdrawal of the semi-solidified shell from the water-cooled mold. In the HDC casting, the surface billet cooling is asymmetric while that of the VDC is symmetric. The HDC casting system involves the withdrawal of the ingot in a horizontal direction ( Figure 1) and gravity acts perpendicular to the direction of casting, with a consequential call for higher casting speed. The flow of the liquid metals could become turbulent as the casting speed increases. The transport of heat within the solidification interfaces can be informed by the fluctuating (or turbulence) flow [1,[4][5][6][7].
Recently, researchers [2,[8][9][10][11][12][13][14][15] have numerically simulated the VDC casting (HDC casting has received less attention) via the low-Reynolds number kinetic energy (k)-energy dissipation rate (ε) (low-Re k-ε) and velocity variance-elliptic relaxation ( − or V2F) turbulence model. The V2F turbulence model consists of four equations, the standard k and ε equations, and the additional and equations. The addition of and is to handle the anisotropic wall-blocking effect, thereby acting as a natural damper closed to the walls [14]. It has been reported that the agitation of the liquid melt, as a result of turbulence eddies, can change the solidification interfaces [14]. Few researchers have numerically investigated flow and temperature distribution in HDC casting. Larouche et al. [16], Krane and Vušanović [17], and Ansari Dezfoli [18] disregarded the turbulence effect in the bulk liquid region of the HDC casting system. Additionally, Han et al. [19] numerically solved the thermal and structural field for HDC casting of copper pipes used in refrigeration and air conditioning systems, but, ignores the fluid flow effect on the solidification interfaces. Recently, the performance of some turbulence models has been investigated as reported by Nzebuka et al. [1] and Waheed et al. Their [1,2] conclusion is that the V2F model is the best turbulence model for casting simulation using Al-Cu alloy.
Al alloys, such as the Al-Cu and Al-Si binary alloys, have been extensively investigated in the DC casting using numerical approaches. One reason appears to be that the fraction of liquid and solid phases for the binary alloy can be easily determined from the phase diagram using a simple linear liquid (solid) fraction-temperature relation [20,21]. Numerical modeling approaches employed in DC casting are two-fold; models that are based on the continuum mixture model of Bennon and Incropera [22] and the volumeaveraged multiphase phase model of Ni and Beckermann [23]. The latter is used when it is necessary to explicitly track the interface between phases and account for the grain motions [23]. The difficulties associated with the volume-averaged multiphase modelling centered on tracking the interface which moves with time and space and the exact locations of the interface during solidification are not known beforehand. Therefore, assumptions about the geometric regularity of the interface phases and moving numerical grids are made. However, experimental observations have indicated that the interface separating the liquid and the mushy region is highly irregular [24,25]. Such irregularity can be attributed to growth rate localization variations and penetration of the bulk liquid to the mushy region by natural convection, thus, causing remelting at the interfaces. Nevertheless, some researchers have successfully implemented the volume-averaged model in DC casting [26][27][28]. The continuum formulation eliminates the need for consideration of conservation equations in each phase separately, in addition, it also eliminates the need for explicitly specifying boundary conditions at the interface of the solidifying domain. The continuum mixture model is valid in the liquid, slurry, mushy, and solid regions of alloy solidification systems [22], and some researchers have recently used this model in DC casting and other metal solidification simulations [1,18,29].
In the present work, the Al-Cu system is used as a working fluid, and the V2F turbulence model of previous work [1] is employed. It is important to restate that the previous work [1] only investigated and qualify different low-Re turbulence models for HDC casting under fixed operating conditions, i.e., no parametric study of casting operating conditions was investigated. In the current study, different HDC operating parameters were investigated, such as casting speed and cooling conditions, and how they influence the temperature distribution, sump length, and flow characteristics during the production of rectangular ingots. This is because cooling conditions and casting speed in the HDC casting majorly influenced the structural formation of the final cast product and the efficiency of the overall production process. The structural formation affects the physical and mechanical properties of the as-cast product. Insufficient cooling and high casting speed lead to an increase in the liquid sump length (depth) and remelting of already solidified structures as a consequence of melt break-out. Therefore, there is a strong demand to determine the relationship that can aid to predict optimum casting parameters, such as the casting speed, that can lead to melt break-out.
The motivation for this research is due to the recent call for industrial optimization of HDC casting using the CFD approach [18]. The current study will address the following research questions: (1) how does an increase in casting speed affect temperature history, liquid sump length ( ), and flow characteristics? (2) how do the cooling conditions at the secondary cooling zones affect the temperature history and liquid sump length ( ), and what is the limiting casting speed over which melt break-out can occur? The second research question has not been addressed for HDC casting in the literature to the best knowledge of the current authors. The model protocol adopted in the current work is useful to casting engineers who always rely on trial and error in choosing a new casting speed whenever a new alloy is to be produced. Hence, with the new information and the casting speed relations, it is possible and easy to predict the operating window over which melt break-out can occur during HDC and further optimized HDC technological process parameters. Additionally, it is imperative to emphasize that, to the best knowledge of the current researchers, no work in the DC casting numerical modeling and simulation has reported a protocol to obtain limiting casting speed relation for melt breakout and the effect of HDC casting operating parameters on the flow field with consideration of the V2F turbulence model coupled with turbulence energy production due to buoyancy.

Numerical Modeling Protocol
The equations governing the fluid flow and heat transfer in the solidification process were based on a continuum mixture model of Bennon and Incropera [22]. The liquid and solid phases were represented with only one fluid flow governing equation [1,2]. The respective phases of the physical properties are represented by the mixture equations, and the liquid fraction relation is a linear function of temperature as represented in [1]. The liquid melt is incompressible, and the turbulence effects are approximated using the RANS model. To simplify the computational domain and minimized the computational cost, the domain simulated is in 2D x-y coordinates. Though, the cooling of the billet is asymmetric and may require 3D simulation, however, the effect of the domain approximation to 2D was negligible as can be seen later in the model validation section (Section 4.2). The solutal buoyancy was not considered because the species equation was not solved in the present work and such an approach has been recently adopted without sacrificing an accurately predicted sump depth and shape as presented in [18]. However, thermal buoyancy was considered and is incorporated in the governing equation using the Boussinesq approximation. The dimensional governing solidification equations were presented in Table 1.

Equations
Continuity: Momentum: Energy: Kinetic turbulent energy: Turbulent dissipation rate of energy: Velocity variance: Elliptic relaxation: where: The details of the meaning of each term (except the terms in the k and ) can be found in [1,2].
The ,which is the turbulence energy production due to buoyancy, was represented as follows [30]: The gravity term g in Equation (8) is taken to be positive (9.81 m/s ). The buoyancy contribution in Equation (6) is governed by and Henkes [30] suggested that should be modeled as follows: = ℎ | ⁄ | (9) where and are the horizontal and vertical components of velocities. In a situation where the flow is normal to the gravity direction = 0, and if the flow is parallel to the gravity direction = 1.

Protocol for Non-Dimensionalization of the Governing Equations
The protocol adopted is as follows; first, formulate dimensionless parameters and numbers; second, use the parameters to transform the governing equations into non-dimensional form; third, dimensionless numbers that are relevant to the HDC casting were used as an input to perform parametric simulation and results analysis; and fourth, the results obtained from the inputted dimensionless numbers were used to derive the limiting casting speed in HDC casting.
The non-dimensional forms of Equations (1)-(7) are obtained by using the following non-dimensional parameters: Energy: * + ∇ * ⃗ * = ∇ (1 + * )∇ * − ∇ * ⃗ − * − ∇ * ⃗ − + * * Kinetic turbulent energy: * + * ⃗ ∇ * = ∇ 1 + µ * ∇ * + µ * * − µ * * − * + * Turbulent dissipation rate of energy: * + * ⃗ ∇ * = ∇ 1 + µ * ∇ * + µ * * * * − µ * * * + * (15) Velocity variance: * + * ⃗ ∇ * = ∇ 1 + µ * ∇ * + * * − * * * − * + * Elliptic relaxation: * − * ∇ * = ( − 1) * * * + * * + * * * Equation (12) indicates that the strength of the viscous forces which is the first term on the right-hand side (RHS) is reduced by increasing the Reynolds number (Re). Increasing Re is equivalent to increasing the inlet . Additionally, the pressure force (second term in the RHS) and the drag force (third term in the RHS) are reduced by increasing the Re. The increase in the Darcy number (Da) is equivalent to increasing the permeability of the mushy zone, therefore its increase in combination with the increase in Re can reduce the third term to zero. Therefore, in the vicinity of the bulk liquid zone, the drag forces are negligible because of the high permeability and velocity encountered. The buoyance force (fourth term in the RHS) and the turbulent flux (fifth term in the RHS) are weakly affected by the Re increase. This is because the Re is balanced with the Grashof number (Gr) and the turbulence viscosity ratio (µ * ) in their respective terms. Therefore, the strength of the buoyance forces is reduced by the Re. The µ * is the ratio of the turbulence viscosity to the molecular viscosity. The molecular viscosity has an inverse relationship with the Re since high Re reduces the viscous shear forces, consequently increasing the µ * .
The first term in the RHS of Equation (13) is reduced with an increase in the Péclet number (Pe). The Pe effect on the energy equation is similar to the effect of Re on the momentum equation. The increase in the Stefan number (St) decreases from the second to the fourth term on the RHS of Equation (13) and vice versa. Before the beginning of solidification, the solid fraction ( = 0) and the second term on RHS are zero and can be ignored. Therefore, the third and fourth terms on the RHS are important and cannot be ignored, because the liquid fraction ( = 1) at the beginning of solidification. However, its magnitude will depend on the value of St. Furthermore, if = 0 (fully solid) at the end of solidification, the third and fourth terms on the RHS are ignored for all values of St, and the second term on the RHS becomes important. The increase in Pe is balanced by the turbulence energy diffusivity ratio * (last term on the RHS) in a way similar to the function of the µ * in the momentum equation.
The non-dimensional * , * and * turbulence equations (Equation (14) to Equation (16)) are controlled by the phenomena of the mean flow gradients which are represented by the momentum equations. Therefore, the effects of Re on the terms are similar to the previous discussion on the last term of Equation (12). The turbulent Grashof number ( ) tends to augment the strength of the third terms of Equations (14) and (15). The is reduced at high Re (high casting speed) where the flow is controlled by forced convection and increased at low Re (low casting speed) where the flow is dominated by natural convection. Therefore, the , Re and µ * are the controlling parameters that increase or decrease the contribution of turbulence energy production due to buoyancy. The last term on the RHS of Equations (14)-(16) vanishes at higher Re and Da. The second term of Equation (17) is reduced at a high turbulent Reynolds number Re and vice versa.

Description of the Mushy and Slurry Zone
The solidification of hot liquid metal in the HDC casting consists of different flow zones ( Figure 2): the bulk liquid zone ( = 0); the slurry zone ( = 0 to 0.27); the mushy zone ( = 0.27 to 0.99); and the solid zone ( = 1). The hybrid model was used to account for the transition from the slurry to the mushy zone [1,2].
The drag forces acting to dampen the fluid flow in the mushy region are embodied by the third term on the RHS of Equation (2), and the term µ can be represented as suggested in [14].
where = (19) = 180 and d is the secondary dendritic arm spacing (400 µm). The ϵ is assumed to be 0.001 in order to circumvent the numerical difficulties that can arise as a result of division by zero. The critical solid fraction used for the transition from the slurry to the mushy zone is 0.27 as suggested in [31]. The 0.27 is the solid fraction value above which a rigid solid structure is formed. The hybrid model removes the need for using separate equations each for the slurry and the mushy zone [5]. The relative viscosity is written as suggested by Nzebuka and Waheed [5]: = 0.27 The relationship between the liquid fraction and temperature is assumed to be linear from a metallurgical phase diagram and is given as follows: = (23) where and are the solidus and liquidus temperatures, respectively.

Specification of the System, Boundary Conditions, and Materials Properties
The schematic of the various flow zones of the computational domain and the BCs are given in Figure 2 and Table 2, respectively. The inlet location can be placed towards the upper or lower position rather than at the central position. The computation geometry ( Figure 2) and the choice of the inlet location are adopted from [1,2,18]. The inlet is 10 mm and the width of the casting domain is 80 mm. Therefore, to ensure that continuity is satisfied, the inlet velocity ( ) is equal to eight times the casting speed. The condition of = 0 at the domain outlet thermal boundary condition seems to be valid at lower casting speeds (low Pe) and longer domains. As the casting speed increases, the condition of zero temperature gradients may not hold since the hot liquid metal may approach the outlet, thereby creating a temperature gradient. However, for consistency of simulation settings, the condition of = 0 will be adopted for all cases investigated. This assumption appears to be appropriate concerning HDC casting since, for all cases simulated, the outlet temperature is well below the solidus temperature (as shown in the later section) of the alloy investigated and any variation in temperature towards the outlet will not affect the solidified casting.

Solution Methods
The governing equations represented by Equations (1)- (7) were solved using a CFD commercial software ANSYS FLUENT 2020 R2, Canonsburg, PA, USA [32]. In the commercial software, the gravity tab was left inactive since the density of liquid and solid are different [1]. Consequently, the fourth term (the Boussinesq approximation) on the RHS in Equation (2) for the y component of the momentum equation was implemented in FLU-ENT using a user-defined function (UDF) code. Additionally, further UDF codes were written for the third term on the RHS of Equations (4) and (5), the mushy region control term of Equation (19), the hybrid models of Equations (20)- (22), and the BCs for the primary and secondary cooling regions. The ANSYS FLUENT graphical user interface for user-defined functions was used to compile and activate all the UDFs in their respective section tabs. For instance, the UDFs for the Boussinesq approximation term were implemented in the cell-zone conditions source term tab and Equation (19) was implemented in the solidification and melting model mushy zone Table. The ( ) term in Equation (18) is inbuilt into the commercial software and it turns to zero at a fully liquid region which switches off Equation (18). As the solidification proceeds, decreases, of Equation (22) tends to 1, increases to a large value (because c is in the order of 10 ), and the A term of Equation (18) increases and dominates the other terms in Equation (2), thereby switching the flow from fully liquid to a mushy and solid zone.
The momentum equation is nonlinear, therefore, the coupling system was resolved using the SIMPLEC algorithm. The transient terms were discretized using the first-order implicit formulation. The calculation was computed until the solution converges at every time step. The thermophysical properties of the material used in the present work are given in Table 3. The adaptive meshing technique was used and it was found that above 16,605 nodes the grids were independent of the mesh size. Hence, for computational saving purposes, the current computation was run using 16,605 nodes.

Numerical Model Validation
The VDC casting had experimental and numerical work documented in the literature. The VDC and HDC have the same model equations. Therefore, the present work model was validated and compared with the experimental and numerical sump profile (shape) reported by Baserinia et al. [8]. They experimented at the Novelis Global Technology Center in Kingston, Ontario, Canada with VDC casting for an AA3003 aluminum ingot.
The initial temperature of the melt was 973 K, a casting speed of 2.33 mm/s, and cooling water flow rates of 1.8 L/s. These data were used in our numerical model and the solution was iterated to convergence. They [8] measured the shape and depth of the liquid sump using a melt poisoning technique with a 50Al-50Zn. Additionally, their [8] numerical solution utilized a 3D Eulerian multiphase model and they ran the equations using an ANSYS CFX commercial package to predict the shape and depth of the sump profile. Figure 3 displays the comparison of the experimentally measured and numerically predicted trace of the sump profile of [8] with the present numerical model predicted result. It is apparent from the figure (Figure 3) that the predicted sump profile for the present model is largely comparable in shape to the measured and predicted sump of Baserinia et al. [8].
The sump shape predicted using the current model tracks closely with that of Baserinia near the nose-shaped region of the sump depth. However, the present model profile widened slightly from the measured and predicted the sump profile of Baserinia in any other region. The discrepancy seen may be associated with the two-dimensional model considered in the current study [1]. Nonetheless, this is at no time a reason for concern since in the DC casting CFD modeling, sump depth (or length) is normally used to measure the possibility of melt break-out and occurrence of macrosegregation [1,12,33]. It should be noted that the sump depth predicted by the present model was 67.50 mm, while the predicted and measured sump depth by Baserinia et al. [8] were 68.60 and 66.52 mm, respectively. This suggests that the predicted sump depth agrees reasonably despite the 2D approximation of the domain. Therefore, it is assumed that the only effect that the 2D approximation may have on the HDC casting is on the sump shape due to asymmetric cooling conditions on billet surfaces.  Table 2. Figure 4 illustrates the velocity vectors placed together with the isotherms of various solidification zones. The fully liquid zone is before the 918 K (liquidus line) isotherm. Between the 918 and 905 K isotherms is the slurry zone. The mushy zone is between the 905 and 843 K (solidus line) isotherms. The fully solid zone is below the 843 K isotherm. The liquidus isotherm is asymmetric with an obvious bend while the solidus isotherm is symmetric about the middle of the casting. This indicates that the asymmetric nature of heat transfer is dominant in the fully liquid zone where gravity has a strong effect with negligible effects on the solid zone. In all cases investigated, the solidus isotherm is below 250 mm (total length is 300 mm) of the casting domain, this justifies why the condition of = 0 was imposed at the outlet. The reference velocity vector arrows placed above the figures were used to compare the magnitude of the velocity vectors inside the domain. The overall flow dynamics during the solidification of the alloy were majorly controlled by forced convection and natural convection. The forced convection is induced from the inlet region and is seen as clustered liquid jets deflected upward towards the top mold region. Natural convection is caused by buoyancy forces and it majorly influences the recirculating vortex seen within the fully liquid and the slurry zones. The figures show that the velocity vector arrows increased with the increase in the Pe. This indicates that the strength of forced convection is stronger than the other transport mechanicisms for the higher Pe.
The effect of increasing the Pe is also seen in the contour of the isotherms. The liquidus (918 K) and solidus (843 K) isotherms extend towards the outlet section of the domain with an increase in the Pe. This results in deeper melt penetration as a consequence of the deepening of the liquid sump length, the slurry, and mushy zones widen toward the central region for the higher Pe. Similar observations have been reported in [1,2] for HDC casting and in [9][10][11][12][13][14]33] for VDC casting at higher casting speeds. This reveals that heat extraction from the molds and solidification rate are lower for the higher Pe than the lower Pe.
The physical implication of the higher Pe is that at higher casting speeds, the production thruput increases and the cost of production lowers, but the possibilities of melt break-out, due to the deeper and stronger melt penetration beyond the mold region, increases. The predicted sump length, which corresponds to the center of the 843 K isotherm (see Figure 4) [18].
The streamlines of the velocity vectors of Figure 4 are presented in Figure 5 and they clearly describe the flow pattern more vividly. The turbulent viscosity ratio was also plotted together with the streamlines to estimate the strength and spatial variation of the turbulence field during the solidification of the Al-Cu alloy in the HDC system.
The flow pattern for all the Pe's is similar with major central and minor top clockwise and counter-clockwise recirculating cells, respectively. An exception in the flow pattern is that at a Pe of 59.11 and 82.75, there is a presence of another minor counter-clockwise recirculating cell seen towards the lower region of the large central clockwise recirculating cells. This particular recirculating cell feature is not present for Pes greater than 82.75 and can be attributed to the presence of a strong convective force of the upstream fluid. The strong jets of hot liquid metal at higher Pes tends to destabilize the recirculating tendency in that region. The figure also reveals that the number of the major central and top recirculating cells increases with an increase in the Pe due to strong forced convection. Further observation indicates that towards the center of the slurry zone (see Figure 4), the streamlines at Pes greater than 82.75 experience flow instability and tend to curve downward. Additionally, the central streamlines for Pes of 153.68 and 177.33 are not essentially straight due to the effect of deep melt penetration.
The turbulent viscosity ratio increases with an increase in the Pe as depicted in the color legend in Figure 5. The turbulent viscosity ratio measures the strength of turbulence in the flow domain [1,2]. The strength of the turbulence is strong towards the top mold due to the influence of turbulent buoyancy production [1]. The spatial variation of the strength of the turbulence field is apparent with an increase in the Pe. This indicates that an increase in the casting speed influences the strength and the spread of the turbulence field within and beyond the mold region.
Temperature distributions are presented in Figure 6 for the six values of Pe. In all the simulated cases, the effect of the forced convective heat transfer mechanism is evident. At a high casting speed (i.e., high Pe), the diffusion term in Equation (13) decreases relative to the advective term. Therefore, heat transfer by conduction decreases relative to convective heat transfer. It should be noted that in metal casting, volumetric heat extraction by conduction is much larger than volumetric heat extraction by convection since the thermal conductivity of the solid Al-Cu alloy is almost twice the liquid melt thermal conductivity. For the smallest casting speed (i.e., Pe = 59.11), the rate of solidification of the metal alloy is efficient since a large section of the domain is solid (depicted by the green-blue region in the figure). As the Pe increases, the liquid region (red-brown) advances spatially towards the mid-section of the domain. To further investigate the effect of the coupling between the momentum equation, energy, and turbulence equation, the x-component velocity, the temperature, and the turbulent viscosity ratio were plotted as a function of . The dimensionless is plotted from the lower to top cooling regions.
The L is the length of the sump and it was used to normalize the distance along the horizontal direction x. The first profile (x/L . ) is plotted in the fully liquid zone close to the inlet region, the second profile (x/L . ) is plotted to cover partly the solid, fully liquid, slurry, and mushy zone, and the third profile (x/L ) is plotted at the end of the iso − curve. The positions of these profiles are presented in Figure 7 which corresponds to a case of Pe of 130.04 with symmetric solidus isotherm (843 K).  The local velocity magnitude towards the top cooling regions for the /L . profile is slightly higher than the /L . profile for Pe values of 59.11 to 106.4. Additionally, the local velocity magnitude is slightly larger than the inlet velocity towards the top cooling regions for the profile of /L . and Pe values of 59.11 to 106.4. This may be due to the larger contribution of the buoyancy forces in the flow since the buoyancy term in Equation (12) is bigger at low Pe values. A similar observation in the increase of the flow velocity due to the buoyancy effect has been reported by Nzebuka et al. [12] in the VDC casting processes. However, the magnitude of their [12] predicted velocity is very much larger when compared to the casting speed since the flow situation considered is buoyancy-assisted flow. However, as the Pe increases more than 106.40, the local velocity magnitude for the profiles of the /L . and /L . becomes equal and approximately assumes the value of the inlet velocity towards the top cooling regions. This is expected since the forced convection dominates the natural convection, thereby decreasing the magnitude of the buoyancy term in Equation (12). Figure 9 presents the turbulent viscosity ratio as a function of Y. It is apparent from the figure that for Pe = 59.11, the turbulent viscosity ratio is lowest towards the top cooling regions for the /L . and /L . profiles when compared to the profiles with higher Pe values. The largest value of the turbulent viscosity ratio is less than 25 for the lowest Pe profiles. With an increase in the Pe, the turbulence viscosity ratio increases towards the top cooling regions to a maximum value of 60 for the /L . profile and Pe = 177.33. For the Pe = 59.11. 82.75 and 106.40 the maximum turbulence viscosity ratio for the /L . profile is greater than the /L . profile. However, at the Pe = 130.04, /L . and /L . profiles tend to equal each other towards the mid-section and top cooling regions of the profile. Then, for Pe = 153.68 and 177.33, the turbulent viscosity ratio is maximum towards the mid-section with /L . reporting slightly higher values than /L . . The turbulence viscosity is less than one for the /L profile for all the values of Pe simulated, and the flow laminarization occurs at this location. Additionally, in all the Pe values and profiles plotted, the turbulence viscosity ratio decreases close to the wall indicating a strong presence of viscous forces. The temperature drop during metal alloy solidification measures the extent of heat extraction in the casting domain. In the present work, the temperature drop was quantified with a dimensionless temperature variable given as follows, T * * = , and plotted in Figure 10. If T * * > 0, then heat removal is poor, and if T * * < 0, heat removal is efficient. Whatever the values of Pe, the T * * < 0 at the walls of the cooling regions. The T * * ≥ 0 for the profile of x/L . for all the values of Pe investigated. The T * * for the /L . profile towards the cooling regions is lower than the /L . profile of all the values of Pe simulated. Additionally, the profile of /L has a T * * curve lower than the other profiles and the relative position of the T * * curves for this profile are approximately constant for all the Pe values simulated. One notable observation in Figure  10 is that for the profile of /L . , the T * * becomes more negative (turns lower) at the wall cooling regions as the Pe increases up until Pe = 106.4. For a Pe greater than 130.04 to 177.33, there is no appreciable change in the T * * at the wall cooling regions. This indicates that during metal solidification in the HDC casting, there is a location where an increase in the casting speed does not affect the wall temperature drop.

Influence of Different Cooling Conditions on the Solid Fraction, Dimensionless Temperature and Velocity, Viscosity, and Sump Length
Further simulations were made for casting temperature ( = 973 ) and Pe = 106.40. The cooling conditions at the primary and secondary cooling regions were represented as Nusselt numbers (Nu = ), where h is the heat transfer coefficient and W is the width of the domain. The protocol established in this section will be used to obtain a limiting casting speed relation for HDC casting.
The Nu for the primary cooling region is represented as Nu and is fixed to a constant value of 10, while the Nu for the secondary cooling region is given as presented in Table 4. Variation in secondary cooling regions was chosen since the majority of heat extraction is in that region [34].   Figure 11 presents the solid fraction evolution and temperature (T * * ) as a function of the vertical distance Y from the center to the top secondary cooling regions for different Nus. The profile covers both the solid and fully liquid zone and the range of Nu values chosen was 10-50. This is the range we observe slight appreciable differences in the solid fraction evolution and the temperature drop. It is evident from the figure that an increase in the Nu increases the formation of solid fractions and decreases the temperature. This means that the rate of solidification increases at higher heat transfer coefficients in the secondary cooling regions. The effect of the Nu number increase extends up to the center of the casting at the fully liquid zone as shown in the zoomed-in portion of Figure 11b. It is obvious from the figure (Figure 11a and b) that the variation of solid fractions and T * * is not pronounced with Nu. Therefore, additional values of Nu were simulated up to 70 to investigate its effect on the sump length ( ) as shown in Table 5.
(a) (b) Figure 11. (a) Solid fraction and (b) Dimensionless temperature as a function of the dimensionless vertical distance at X = 1. X is the dimensionless horizontal distance.
It is apparent from the table that an increase in the Nu from 10 to 40 appreciatively decreases the sump length. Above the Nu value of 50, the increase in the Nu weakly affects the sump length. The increase in the sump length from the Nu value of 60 to 70 is approximately 1.1%, which can be taken to be a minor variation. 108.34 Figure 12 shows the behavior of the as a function of Nu. The curve tends to be asymptotic at Nu greater than 50. This corresponds to a limit where an increase in the cooling water flow rate has no appreciable effects on the sump length and the volumetric heat extraction as also reported in [33]. To obtain the sump length for any given values of casting speed and cooling conditions represented by the heat transfer coefficient or Nu, / was plotted as a function of the Nu (see Figure 12b). An exponential asymptotic nonlinear curve fitting can be written as follows: = Pe(1.025 + 0.863 * (0.923) (mm) (24) The relation in Equation (24) can be exploited to provide the correct casting condition for the HDC casting. If the is replaced by the entire length of the HDC casting L (mm) in Equation (18), the maximum inlet velocity and casting speed above which melt break-out can occur is obtained as follows: (m/s) = ( * ) (25) where = (1.025 + 0.863 * (0.923) (mm), is the width of the casing (m). Therefore, the maximum casting speed is written as: (m/s) = ( * ) (26) where = .
Equation (26) predicts the maximum casting speed for the HDC casting system above which melt break-out may occur. Equation (26) is valid for a contact mold heat transfer coefficient above 2000 W/m 2 K, which is common for the most practical situation. It should be noted that Equation (26) gathers information about the properties of the material, cooling conditions in the secondary cooling regions, and the geometry of the HDC casting domain. Therefore, Equation (26) may be applied to a wide range of HDC casting alloys and can be used to predict a wide-casting parameter operating window before melting break-out occurs. Further simulations indicate that an increase in the casting temperature from 943 K (equivalent to 10 K superheat) to 996 K (equivalent to 60 K superheat) only results in a sump length increase of 2.29%. This reveals that an increase in the has a trivial effect on the sump length in the HDC castings system, which is also true for the VDC system as reported in [10,12].

Conclusions
The present paper focuses on the computational fluid dynamics (CFD) prediction of the flow redistribution and thermal history in a horizontal direct-chill casting of an aluminum-copper alloy system. The governing solidification equations and the turbulence models were non-dimensionalized to understand the relevant contribution of each term in the solidification processes. With an increase in the Pe (equivalent to an increase in the casting speed), the liquidus (918 K) and solidus (843 K) isotherms extend toward the outlet section of the domain. This results in deeper melt penetration as a consequence of the deepening of the liquid sump length, and the slurry and mushy zones widen toward the central region for higher Pes.
An increase in the Nu increases the formation of solid fractions and decreases the temperature. This means that the rate of solidification increases at higher heat transfer coefficients in the secondary cooling regions. Additionally, an increase in the Nu from 10 to 40 appreciatively decreases the sump length. Above the Nu value of 50, the increase in the Nu weakly affects the sump length. The increase in the sump length from the Nu value of 60 to 70 is approximately 1.1%. A simple maximum casting speed as a function of the geometrical data and the alloy properties relation was obtained. This relation is useful in defining a wide range of casting speeds for different alloys and geometrical data.
The current work is useful to casting engineers who always rely on trial and error in choosing a new casting speed whenever a new alloy is to be produced. Hence, with the new information and the casting speed relations, it is possible and easy to predict the operating window over which melt break-out can occur during HDC. Further study may compare the sump shape predicted by the 2D simulation and that of the 3D one, investigate the sump shape effect on the surface and centerline macrosegregation by considering species transport, and use an experimental method to verify Equation (26) for the limit of casting speed which will lead to melt break-out. Funding: This research received no external funding Data Availability Statement: The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.