# Computational Evolving Technique for Casting Process of Alloys

Authors: Amir M. Horr
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The challenging task of bringing together the advanced computational models (with high accuracy) with reasonable computational time for the practical simulation of industrial process applications has promoted the introduction of innovative numerical methods in recent decades. The time and efforts associated with the accurate numerical simulations of manufacturing processes and the sophisticated multiphysical and multiscale nature of these processes have historically been challenging for mainstream industrial numerical tools. In particular, the numerical simulations of industrial continuous and semicontinuous casting processes for light metal alloys have broadly been reinvigorated to investigate the optimization of casting processes. The development of advanced numerical techniques (e.g., multiscale/physical, finite zoning, and evolving domain techniques) for industrial process simulations including the transient melt flow, heat transfer, and evolution of stress/strain and damage during continuous casting processes have endeavored many new opportunities. However, smarter and broader improvements are needed to capture the underlying physical/chemical phenomena including multiscale/physical transient fluid-thermal-mechanical coupling and dynamic heat-transfer changes during these processes. Within this framework, the cooling system including its fluid flow and its characteristic heat transfer has to be modelled. In the research work herein, numerical studies of a novel transient evolving technique including the thermal-mechanical phenomena and Heat Transfer Coefficient (HTC) estimation using empirical and reverse analyses are presented. The phase change modeling during casting process including liquid/solid interface and also the implementation of dynamic HTC curves are also considered. One of the main contributions of this paper is to show the applicability and reliability of the newly developed evolving numerical simulation approach for in-depth investigations of continuous casting processes.


---

## Body

## 1. Introduction

The speed and the quality of industrial production processes are totally depending on efficient and sound material processes based on rigorous scientific and engineering knowledge, experiences, and numerical simulations. For many light metal alloy industries, it is required to establish a framework for continuous innovation in design and productions and a sound and comprehensive strategy for confronting crucial challenges and competitions. The combination of state-of-art innovative mathematical, physical, and phenomenological models (with accurate predictive powers) and sound parallel experimental and numerical simulation studies to explore what technical solutions are viable for these processes is not yet exploited to the extend it could. Today’s integrated computer-based simulation tools which consist of advanced/practical multiphysical/scale models (continuous and discrete), data-driven schemes (Artificial Intelligence, AI, Machine learning, ML), and pioneering integrated process simulation tools would reflect the essence of whole engineering and industrial experiences and knowledge. These new industrial design and production tools can be developed in a systematic way to handle vast related processes/databases, lightweight industrial experiences/knowledge-bases, and information [1]. One of the most important ingredients of these integrated systems is numerical simulation tools where properties and the performance of production chain processes based on conventional and/or advanced techniques are implemented.The combination of physical, mathematical, and computational topics to approach the fundamental question of optimum process parameters and conditions has been promoted in recent years. Instead of having single numerical domain with its associated grid (mesh), the emerging idea includes a new approach to represent various complexities (e.g., change of phase, micro cracks, etc.) within limited dynamic zones in the computational domain that works for both implementing advanced mathematical models (e.g., solidification, plasticity, etc.) and also computational time. This provides the basis for multiple numerical domains (zones), running in parallel (with multiple instances of solvers) and their communication using advanced interface routines (using mapping, condensation, etc., techniques). The interaction/communication of these evolving zones within the simulation system can be managed smartly using advanced controllers (AI and ML technologies) and optimum design can be carried out using the combination of analytical and smart data-driven schemes.One of the main concerns for these new-generation techniques is how to represent changes of the dynamic system (e.g., change of phase, cracking, nonlinearities, etc.) to understand the consequences of solution and optimization techniques on the process. The overarching research work herein is a step towards the founding of the procedure to combine complexity-driven solution techniques (e.g., nonlinearity within dynamic zones) with efficient interfacing and computationally efficient solution schemes. The technique relies on separate and parallel efforts devoted to predicting the response of complex system (e.g., solidification and heat energy evolution during casting process). The research work also combines different aspects of thermal-fluid-mechanical interactions of industrial processes to predict optimum parameters and avoid damage and failure through the implementation of multiscale (e.g., micro-macro) mechanical models at selected dynamic zones.As part of the integrated simulation approach for multiphysical, multiscale simulation, the development of cooling zones at boundaries, dynamic mesh modules along with parallel processing, and computational efficiency have also been promoted herein. The purpose of the research work is to present steps towards an integrated numerical method, which would help improve the simulation of dynamic material processes and their optimization for the future manufacturing.


## 2. Conventional Modelling Concepts

The life-long performance of engineering components along with their costs, agility, and versatility requires sound and practical processes within manufacturing chain. To help optimizing these processes for different light metal alloys, various numerical techniques/tools have been developed and employed in the last sixty years. The early numerical approaches developed during second half of the 20th century (e.g., Finite Element, FE, etc.) have employed simple discretization and solver techniques for static and dynamic processes including fixed and adaptive meshing [2, 3]. Additionally, different types of 2D and 3D elements along with variety of shape functions have been developed to discretize numerical domains for simulation of industrial processes including casting, extrusion, rolling, etc. as static, pseudostatic and/or full dynamic processes (e.g., element activation, and deactivation). The discretized domain is then solved at sequential time steps (or frequency steps) using conventional time/frequency integration techniques. The time-history results are then assembled using deterministic and/or statistical schemes to study the dynamic behavior of industrial processes under real-life scenarios [2, 3].For industrial casting applications, the fixed Lagrangian and/or Eulerian [4] discretized techniques have traditionally been used where the variation of material confined by the numerical grid is solved at different points in time (time steps). In the conventional pure Lagrangian and also updated Lagrangian approaches (with adaptive remeshing), the process of large deformation (or flow) of a body is resembled as a material flow where each material particle (cell) transports its own properties (such as microstructure, temperature, density, etc.). For the simulation of dynamic systems, as the deformation/flow front advances, its properties may change in time. However, their inability of following large distortions of the computational grid without recourse to frequent remeshing (e.g., adaptive remeshing) schemes [4] is its main shortcoming (it is a computationally expensive task). For the Lagrangian approach, let us consider two overlapping domains as a physical material domain, RX, and a spatial domain, Rx (spatiotemporal conditions are only considered at discrete time steps). The mapping function φ for these two overlapping domains can be defined as [5](1)φ:RX×t0,tT→Rx×t0,tT(2)X,t→φX,t=x,twhere t0,tT are time at original and updated (terminal) configurations and mapping at discrete time points can be allowed as(3)x=xX,tt=twhich states that the spatial coordinates depend and also strictly track the material particle at discrete time points (no convection). The definitions of time at both domains are real and no hypothetical time definition/interpretations are carried out in a single numerical Lagrangian domain.Alternatively, in the Eulerian fixed and moving mesh approaches, rather than tracking each material cell (e.g., at deformation front), the evolution of the material flow properties at every point in space (vortices) can be recorded at sequential time steps. This means that the material flow properties at a specified location within the grid depend on its location and time (advection). In this approach, the conservation of material is satisfied at spatial and time points and the computational mesh is mostly fixed where material moves with respect to the grid. The large distortions (and melt flow) of the material motion can be handled with relative ease using convective feature of the method where relative motion of deforming front in respect to the discretized grid is calculated [6–10]. Figure 1 shows the Lagrangian and Eulerian deformation concepts along a with typical temperature field for a 2D Eulerian simulation of material flow in horizontal casting.Figure 1
(a) Reference grid and deformation in (b) Lagrangian and (c) Eulerian schemes, (d) temperature field for 2D Eulerian material flow in horizontal casting simulation.
(a) (b) (c) (d)Many of the mainstream industrial casting simulation tools have combined the benefits of both techniques in the way that the simulations of melt flow into the mold (filling simulation) are carried out using Eulerian technique, while the thermal-mechanical simulations for the solidified billet are performed using Lagrangian fixed-meshed technique [5]. The Arbitrary Lagrangian Eulerian (ALE), Mixed Lagrangian Eulerian (MiLE), and other similar hybrid methods have been used in these tools where benefits of both numerical schemes are combined [9, 10]. For the ALE approach the material motion/deformation and the spatial configuration are referenced to the third overlapping domain through using mapping functions. The general practice for the discretized grid is to have a mapping function for the material configuration, spatial coordinates, and the third domain defining the grid points. Three different mapping transformations are defined to overlap the numerical domain where the motion of the grid points can be distinguished from the material motion. If the mapping from the third referential domain to the spatial domain is defined by ϕ, then [5](4)ϕ:Rχ×t0,tT→Rx×t0,tT(5)χ,t→ϕχ,t=x,twhere Rχ is referential domain. The grid motion (e.g., casting speed or billet speed) can then be defined as(6)vgχ,t=∂x∂tχwhere both material deformation/motion and grid motion are calculated with respect to real time during the simulations. Figure 2 shows the time snap shots of Eulerian mold filling simulation along with subsequent Lagrangian thermal-mechanical simulations for an aluminum semicontinuous industrial casting [11].Figure 2
(a) Eulerian filling simulation of mold and (b) subsequent Lagrangian extension of the cast billet with their temperature contours for double-symmetric quarter model.
(a) (b)The computational concepts for dynamic material processes have been going through fundamental changes since the start of 21st century by introduction of dynamic discretization and multiscaling concepts [12]. The novel ideas of evolution of material properties, updating domain’s discretization/boundaries and also computational partitioning during the numerical simulations have been developed [12, 13]. Furthermore, the introduction of material evolutions and their associated physical/mathematical models at different scales (macro, meso, micro, etc.) has helped to overcome numerical problems related to resolution and accuracy of dynamic processes [12]. For the material continuous casting case, the transient nature of the numerical domain and its continuous expansion throughout the process has conventionally been handled using either a birth and death approach or the splitting of element layers at the stationary-moving interface [14]. However, since the change of geometry, mesh, boundary, and energy source within dynamical domain (e.g., continuous casting simulation) can have profound effects on the computational accuracy/time, new schemes had to be developed for higher performance, more generality, and broader versatility.In the proposed dynamic multiscale and multiresolution approach herein, new subdomains/zones are generated at specific time points to resemble the extension of the cast billet. These newly generated zones are overlapping/attached to the main domain at discrete time points where comprehensive mapping procedure (including grid overlapping, nodal condensation, propagation interface, etc.) are applied to reassemble the domain to handle the new material and energy input (new melt and solidified melt material). The added dynamic subdomains/zones can be stationary or moving depending on requirement of the numerical model for the industrial process. As the numerical solution is performed on a full parallel-processing machine, multiple instances of the same solver (or alternatively different solvers) can independently be employed for the base domain and attached/overlapping zones. The generated subdomains (which can be at different length/time scales) would gradually integrate into the main domain (and become a part of main domain) as the evolution phases within the subdomains/zones come to the end (e.g., solidification and nonlinearities).For the interfacing between the main domain and subdomains (or between subdomains) and to overlap grids with multiresolutions and/or multiscale natures (e.g., overlapping meshes at different length scales) some mathematical representations have already been developed. The Partition of Unity Finite Element Method (PUFEM), Overlapping Sphere Elements (OSE), and few other techniques have already been proposed to handle the grid-inconsistency [15]. For overlapping zones within the numerical domain, let us assume two interfacing grids with overlapping at the boundaries of Ω1 and Ω2 numerical domains. For linear finite element grids with base functions vii∈I and wj∈J, the overlapping set can be written as [15](7)∅1vi,∅2wji∈I,j∈JIf the overlapping grid zone is represented with a linear independent relation as(8)∅1∑iαivi+∅2∑jβjwj=0and if it is assumed that ∅2=1-∅1 then it can be shown that(9)∅1∑iαivi-∑jβjwj=-∑jβjwj→∑iαivi-∑jβjwj=Const.and the stiffness matrix for the interfacing/overlapping grids can be assembled for firstly nonoverlapping elements (as usual) and, secondly, overlapping elements with multiple entries. Figure 3 shows the typical partial overlapping zonal grids which can represent different part of a numerical domain. More comprehensive discussions about the multiresolution interfacing/overlapping grids can be found in the literature [15–17].Figure 3
Overlapping zonal representation for nonmatching FE grids.


## 3. Finite Zonal Concept and Interfacing

The investigation of dynamic numerical simulation for industrial continuous casting processes where geometry, mesh, and boundaries of the domain can evolve during the solution has been considered in the last twenty years and different evolution schemes have been proposed [11]. The zonal evolution and smart subdomain technique have been proposed recently and the application of their techniques for the industrial continuous casting process has been investigated [11]. The introduction of dynamic boundaries, evolution zones, and physical/mathematical sink-sourcing energy concepts can facilitate steps towards the formulation of the next generation of numerical techniques for practical and industrial applications.


### 3.1. Multiresolution Evolution Scheme

The multiscale and multiresolution prerequisites of many industrial process simulation platforms have prompted the use of domain evolution concepts during numerical simulations. The research work herein briefly describes some aspects of basic concepts for innovative evolving domain technique for the full-parallel processing. The newly proposed approach would treat the evolving/progressing cast billet during the continuous casting simulation as a dynamic zone. The zone can move in predefined or calculated manner which would be attached to the main domain through overlapping boundary concept. The generated zones can gradually join the main domain (and become a part of main domain) or stay independent for further processing as flagged evolution phases within the zones (e.g., deformation, solidification, etc.) come to the end.For the continuous casting simulation (vertical or horizontal), a directional boundary evolution scheme has been implemented using extensive in-house coding. The technique can be summarized in the following steps.(i)
The initial starting conditions are modelled using the results of either Thermal-Computational Fluid Dynamics (CFD) melt delivery simulations or filling-thermal assumptions.(ii)
The initial geometry of the main domain (i.e., mold) and its discretization are defined and initial system matrices are assembled.(iii)
The initial time steps/iterations are solved (using a thermal-mechanical solver) and based on the casting speed the billet domain coordinates are updated.(iv)
The geometry and discretized grids are modified and adapted by adding a single/multiple layer mesh zone based on the evolution of the domain (updated spatial coordinates and shrinkage and thermal deformation changes).(v)
The system matrices are modified/appended and the input energy is distributed amongst the newly generated subdomain.(vi)
Time-history results (from previous time steps) are mapped to the newly generated subdomain.(vii)
The energy balance (thermal energy) is achieved through the energy sink/source concept and the previous converged solution is used as a first step for a newly updated domain.(viii)
The simulation scheme continues with the new geometry/mesh till the next evolution step is triggered.The lower-scale material evolution phenomena (e.g., microstructure formation during solidification) and its coupling with macro-scale properties can also be considered using numerically efficient overlapped grid technique (not presented in this paper). Figure4 shows snap shots of a casting numerical simulation of simple double-symmetric rectangular cross-section billet using dynamic evolution technique.Figure 4
Dynamic evolving technique for double-symmetric rectangular cross-section billet using directional evolving scheme.


### 3.2. Simulation of Lower-Scale Phenomena

The solidification of melt during the casting process can be described in different levels of scales/details. Since most of the conventional casting simulations are carried out at macro-scale, many of the conventional numerical solidification modules are based on the macro-kinematics of material behavior. For these macro-scale simulations, there are different approaches for defining the phase transition. Most common is the application of the Scheil-Gulliver equation [18, 19] where a nonequilibrium approach is formulated for solute redistribution during solidification (no solute would diffuse back into the solid). The solid fraction in this method is defined by local current, solidus, and liquidus temperatures as well as the partition coefficient. The Scheil-Gulliver equation in its simple and standard form can be written as [18, 19](10)CL=C0flk-1CS=kC01-fSk-1where CL,CS, and C0 are liquid, solid, and initial compositions and fS,fS, and k are solid and liquid mass fractions and coefficient k=CS/CL (which can be determined from phase diagram), respectively. The formulation is valid for a local equilibrium state between the solid and liquid interface and for moderate-to-low cooling rates. The numerical methodology for the simulation of solidification front during casting process where phase changes and different phases exist (liquid, mushy, and solids) can be challenging. Different numerical methodologies have been developed [11] to handle the phenomena within numerical domains, namely, conduction-based, conduction-convection based, and enthalpy-based methods. In the conduction-based method, the energy transfer during solidification is calculated based on conduction heat transfer and to account for the convection heat transfer effects (melt flow), an artificially high fluid conductivity is defined. While in the conduction-convection method the melt convection is taken into account explicitly with mass transfer and thermal convections are formulated using governing differential equations for the melt flow. The third alternative way is to use the enthalpy concept. Enthalpy is a thermodynamic quantity defining total heat energy of a system which can be combined with other numerical features to account for change of phase during solidification. In the mainstream industrial software tools for casting simulations, the enthalpy method is generally combined with Volume of Fluid (VOF) and domain porosity features (to damp the fluid momentum when it turns into solid). The VOF method can trace the solidification front within the numerical domain while the domain porosity would enable the smooth change of phase from liquid to solid by damping the momentum equation using momentum sink terms. The introduction of momentum sink terms in the energy equation might introduce more numerical challenges since the convergence of the system becomes more difficult. Consider the mass continuity and momentum equation for the melt flow inside the mold as [4](11)∇.U=0,∂ρU∂t+∇.ρU⊗U=-∇p+∇.μ∇Uwhere ρ,μ,p, and U are fluid density, viscosity, pressure, and velocity, respectively. Momentum sink term is traditionally inserted into melt momentum equation (for temperatures between liquidus and solidus) to damp the melt velocity due to solidification phenomena. A simple and practical approach for momentum damping can be defined using the Kozeny, Carman approach, and its modification by Blake [20–23]. Let us assume the basic equation for fluid flow in resistive media as(12)∆pL=α1-kβk3μvd2where ∆p/L, k, μ, d, and v are pressure gradient, domain porosity (liquid volume fraction), melt viscosity, primary channel width, and melt superficial velocity. Variables α, β can be defined for material processes during solidification based on the alloy and process properties. Two different types of momentum sink terms can be defined for the melt solidification, namely, viscous and quadratic terms. The basic viscous sink term for fluid flow in a porous media can be defined as [22](13)LV=μKUwhereK is permeability of the resistive domain. Voller et al. have modified [24–26] the fluid viscous loss model and the permeability of the domain as(14)Lq=-μK.u-uswhere us is the solidified slab speed (can be assumed as casting speed) andK is permeability which can be calculated based on the resistance against the melt flow inside the dendrite channels during solidification. Ignoring the remelting and the solute concentration effects (temperature exchange at interface, etc.), Gu et al. [27] have proposed the following equation for the permeability:(15)K=Cfl3fs2λ12where fl,fs, and λ1 are volume fractions of liquid and solid phases (at the solidification front) and primary dendrite arm spacing. The constant C can be assumed as C=6×10−4 [27]. The secondary arm spacing and its resistance effect has not been explicitly defined in (15). The broader equation for estimation of permeability has been derived by author using the combination of physical and analytical approach as [28](16)CLV=C∗λ21-fTαfTβfT=C1∗1-T-Ts∆Tγ∆T=TL-Tswhere C, C1 are constant coefficients, λ is the primary dendrite arm spacing, andα,β, γ are power coefficient which can be calibrated for the alloy processes (for aluminum alloysα= 3, β=2, and γ =1). Coefficients C and C1 can be assumed as C= 0.0006 (similar to Gu et. al) and C1=1 - 1.2 [28]. A more detailed version of the simplified loss term can also be developed [28] using primary and secondary dendrite arm spacing and melt isolation phenomena. Similarly, based on physical concepts of energy damping and a mathematical formulation of melt flow kinetic energy, a quadratic sink term can also be derived as [28](17)CLQ=C.β.1-T-TETE-TSαTE<T<TLα=1β+C.1-β.1-T-TSTE-TSαTS<T<TEα=2where TE is a temperature which denotes the starting point for dendrites bridging phase during solidification. For aluminum alloys, the simple assumption can be that TE=TL-TS/2 and β=γ.TE-TS/TL-TS where γ is a coefficient which depends on material microstructure development during solidification (when growth and progression of dendrites starts to create significant resistance against the melt flow near the interface). For aluminum alloys where equiaxed grain growth are perceived, it can be assumed as γ = 0.6. The coefficientC can be calculated/measured and calibrated from the total momentum of the melt within the mold prior to the solidification phenomena. More comprehensive discussions about these mathematical derivations of momentum sink terms can be found in [28].


### 3.3. Simulation of Cooling and Heat Transfer

Using water spray cooling in light metal alloy productions (e.g., continuous and semicontinuous casting processes) is popular and the numerical modelling of these processes is essential in developing a deep understanding of the thermal evolution and temperature fields (including thermal stresses) during industrial processes. For the casting process, the melt material flowing into the mold can be considered as a viscous liquid which confines in a tiny solidified shell formed during initial cooling (during convective air cooling). In many casting processes, the cast billet is then sprayed (or submerged) with a water flow to extract more thermal energy from inside the billet. During this secondary water cooling, the billet undergoes large temperature fluctuations, especially near its surface. Here it is essential to investigate and control the thermal event to minimize the occurrence of defects on the surface or subsurface (hot tearing, cold cracking, etc.). These thermal fluctuations arise due to primary air cooling at the top (thermal convection and radiation), water spray cooling (including shrinkage and thermal conduction throughout the body of the billet), and possible secondary air cooling.Numerical modelling/estimation of water spray cooling effects is one of the most challenging parts of the cooling system design. It would give rise to the estimation of thermal evolution and thermal fluctuations including its associated stresses. Many thermal evolution studies have been carried out in last two decades to develop an estimation technique for the dynamic HTC on the cast billet surface considering different water boiling regimes (nucleate, unstable, and stable filming regimes) and also Leidenfrost effects. One of the simplest models presented in the literature to start modelling of the water spray systems is to define a spray wet area as [29] (18)A=fx,y,zQVswhere fx,y,z, Vs, and Q are coordinate distribution function, surface flow rate, and water flow rate, respectively. Different parameters would affect the thermal energy dissipation during water spray cooling, namely, temperature difference between billet surface and water, water flow rate, angle of impingement, water quality, surface roughness, and some other minor parameters. The boiling regime of the water on the surface of the billet at high temperature has profound effects on the heat dissipation and the complex interface between alloy/steam/water would determine the real cooling HTC. Many researchers have tried to carry out the inverse heat transfer calculation using available experimental data and fitting techniques [30–32]. However, in recent years the combination of empirical and numerical techniques has been used to model the surface boiling regimes and HTC calculations [29]. To start with the physical-mathematical study of HTC calculation, let us start with the simplified governing equations for the water spray cooling system which can be written in a simple form assuming the incompressibility as [31](19)∇.U=0∂U∂t=-∇Pρ+1Re∇2Uwhere U,P,ρ, and Re are velocity vector, pressure, density, and Reynolds number,respectively. After some mathematical manipulation, for the case of continuous liquid phase (water) and dispersive gases (vapor and air) using the three-phase continuity equation, the phasekcan be written as [33](20)∂αk∂t+∇αkU→k=Γki-Γikρkand the total interfacial force in the momentum equation can be elaborated to find the bubble generation and movement relative to the billet surface (see Figure 5). For the water spray cooling system during casting process, the total heat flux can be defined as [33](21)qT=q1+q2+q3where q1,q2, and q3 are single-phase convective, evaporation, and quenching heat flux densities, respectively. Different models for the range of bubbling schemes in water flow boiling regimes have been developed in the past decades. In many of these models the condition that a vapor bubble would grow if the surrounding liquid temperature is equal to or more than the saturation temperature of the vapor inside the bubble is satisfied. In the research work herein, attempt has been made to formulate new mathematical models for the HTC at different boiling regimes during water spray (and submerge) cooling. For the initial air cooling during casting, conventional HTC approximation technique which includes convection and radiation terms can be adopted. If the temperature of billet surface is assumed asTb and air bulk (free-stream) temperature isTa, the thermal resistance (or inversely conductivity) at the interfaceAT can be defined as [34](22)R=1hr+hcATwherehr andhc are, respectively, the radiation and convection HTCs, calculated as follows:(23)hr=εσT4per unit areawhere σ is a constant (Stefan–Boltzmann constant) which can be assumed as 5.672×10−8 (W/m2 K4),T is the absolute temperature in Kelvin and ε is the billet surface emissivity. The convective HTC (hc) for air cooling can also be calculated using [34](24)hc=kNuλwhereNu is Nusselt number, k is the air thermal conductivity, and λ is the characteristic length. The Nusselt number (the ratio of conductive thermal resistance to the convective thermal resistance of the air) can be calculated as (25)Nu=CGrPrnGr=gβTb-Taν2λ3Pr=νkcpwhereCandnare constants andPr is the Prandtl number (defined as the ratio of momentum diffusivity to thermal diffusivity). Gr is the Grashof number (defined as the ratio of the buoyancy to viscous force acting on air) for vertical flat plates, β is the coefficient of thermal expansion (approximately 1/Ta, for air), ν is the kinematic viscosity, and g is gravitational acceleration. Figure 5 shows the schematic bubble generation and separation event during surface cooling and temperature dependent radiation HTC for aluminum alloy.Figure 5
(a) Bubble generation and separation mechanisms from hot surface and (b) radiation of heat from aluminum billet at different temperatures.
(a) (b)The water cooling HTC values change rapidly during the casting startup condition where continuous cooling regime adds to the modelling complications. The HTCs generally rise rapidly with the billet surface temperature until a critical temperature is reached above which unstable and film boiling regimes, with their lower associated HTCs’ starts. The estimation of critical heat flux (CHF) temperature point and its associated heat flux can help to understand and optimize the billet surface cooling rates and thermal-mechanical stress/strain conditions. The variation of HTC against the billet surface temperature can be divided into four distinctive zones, namely,(i)
below water saturation temperature,(ii)
between water saturation and CHF temperature (nucleate boiling regime),(iii)
from CHF temperature to Leidenfrost temperature (unstable boiling regime),(iv)
above Leidenfrost temperature (film boiling regime).Some aspects of these boiling regimes have already been investigated [34] where mathematical formulation and also empirical estimation of heat transfer have been developed. For the work herein, different sets of mathematical works and experimental verification for the heat transfer during cooling process have been carried out to optimize the HTC during casting [35]. For the preboiling regime a simplified approach can be employed (the preboiling has minor thermal energy exchange effect) and HTC can be defined as [35](26)HTCpre=αQLbβ7e2∗Tb-1.67e3whereQ,Lb, Tb are water flow rate, billet perimeter, and billet surface temperature, respectively. The parametersα, βthat can be defined using the impingement conditions are billet surface quality (for aluminum billetsα, β can be assumed as α = 1.35 and β = 0.3). For the more challenging nucleate boiling regimes, a more detailed formulation has been sought [35] and HTC can be defined as(27)HTCnu=μHfgTsatcfαCpTbHfgασμCpKgρw-ρvβwhere μ, Hfg, cf, Cp, K, g, ρw, ρv, and σ are coolant dynamic viscosity, latent heat of evaporation, flow heat parameter, specific heat, conductivity, gravitational acceleration, density, and vapor density, respectively. Tsat and Tb are defined as coolant saturation and billet surface temperature values where spatiotemporal distribution of temperatures on billet surface can be calculated using numerical domain. α and β are again cooling system parameters (these parameters can be assumed as α =0.8 and β =-0.515 for aluminum alloys).The third boiling regime during cooling process is the unstable film boiling where the interphase boundary is unstable and vapor film can be broken periodically for water to make contact with billet surface. The estimation of HTC for unstable boiling zone has been a challenge for both scientist and engineers trying to overcome the full dynamic and oscillatory nature of the boiling event. As the billet surface temperature decreases during cooling process, the duration of vapor near the surface reduces, while the number and duration of liquid-wall contact increase. Such sporadic on/off contacts of liquid water with the billet surface result in an abrupt discharge of vapor from the surface, causing splashing and turbulence (and driving away vapor from the wall). The HTC for this zone can be estimated using the fraction of billet surface wetted by liquid water during the cooling process. The proposed empirical formulation for this zone can be defined as [35] (28)q=qwAc+qv1-AcHTCtra=HTCnu1-Tb-TCHFTbγwhere TCHF is water critical heat flux temperature (can be assumed as TCHF=185°C for aluminum alloys). The final boiling regime during cooling process is the stable film boiling zone where a vapor blanket would be covering the surface of the billet. Generally, the estimated HTC values by the water spray is accredited to the intensified forced convection in the liquid film caused by the large flow momentum in the gravitational direction and consequent improvements in nucleate boiling. It also promotes the transient conduction which occurs on the heated surface due to the quick refresh of the liquid film. The proposed HTC estimation formula for cooling zone IV (film boiling regime) during cooling process can be proposed as [35](29)
H
T
C
f
=
H
T
C
t
r
a
1
-
T
b
-
T
L
e
T
b
H
1
+
H
2
ƙ
(30)
H
1
=
γ
T
b
×
k
w
3
g
ρ
w
H
f
g
ρ
w
-
ρ
v
1
+
c
p
T
b
/
2
H
f
g
μ
w
T
b
σ
/
g
ρ
w
-
ρ
v
β
α
(31)
H
2
=
5.4
×
10
4
T
b
1
+
0.055
T
b
v
s
κwhere TLe,Tb,μw,Hfg,Cp,σ,ρw,Kwρv, are water Leidenfrost temperature (≈400°C), billet surface temperatures, dynamic viscosity, latent heat of evaporation, specific heat, surface tension, density, and thermal conductivity, respectively. ρv is the vapor density and α,β,γ,κ, and ƙ are film cooling system parameters (these parameters can be assumed as α =0.4, β =0.5, γ=0.5, κ=0.6, and ƙ=0.15 for aluminum alloys). The film HTC calculation is valid for surface temperature (Tb) between 400°C and 550°C. More detailed discussion about numerical modelling of boiling and bubble generation can be found in [35].


## 3.1. Multiresolution Evolution Scheme

The multiscale and multiresolution prerequisites of many industrial process simulation platforms have prompted the use of domain evolution concepts during numerical simulations. The research work herein briefly describes some aspects of basic concepts for innovative evolving domain technique for the full-parallel processing. The newly proposed approach would treat the evolving/progressing cast billet during the continuous casting simulation as a dynamic zone. The zone can move in predefined or calculated manner which would be attached to the main domain through overlapping boundary concept. The generated zones can gradually join the main domain (and become a part of main domain) or stay independent for further processing as flagged evolution phases within the zones (e.g., deformation, solidification, etc.) come to the end.For the continuous casting simulation (vertical or horizontal), a directional boundary evolution scheme has been implemented using extensive in-house coding. The technique can be summarized in the following steps.(i)
The initial starting conditions are modelled using the results of either Thermal-Computational Fluid Dynamics (CFD) melt delivery simulations or filling-thermal assumptions.(ii)
The initial geometry of the main domain (i.e., mold) and its discretization are defined and initial system matrices are assembled.(iii)
The initial time steps/iterations are solved (using a thermal-mechanical solver) and based on the casting speed the billet domain coordinates are updated.(iv)
The geometry and discretized grids are modified and adapted by adding a single/multiple layer mesh zone based on the evolution of the domain (updated spatial coordinates and shrinkage and thermal deformation changes).(v)
The system matrices are modified/appended and the input energy is distributed amongst the newly generated subdomain.(vi)
Time-history results (from previous time steps) are mapped to the newly generated subdomain.(vii)
The energy balance (thermal energy) is achieved through the energy sink/source concept and the previous converged solution is used as a first step for a newly updated domain.(viii)
The simulation scheme continues with the new geometry/mesh till the next evolution step is triggered.The lower-scale material evolution phenomena (e.g., microstructure formation during solidification) and its coupling with macro-scale properties can also be considered using numerically efficient overlapped grid technique (not presented in this paper). Figure4 shows snap shots of a casting numerical simulation of simple double-symmetric rectangular cross-section billet using dynamic evolution technique.Figure 4
Dynamic evolving technique for double-symmetric rectangular cross-section billet using directional evolving scheme.


## 3.2. Simulation of Lower-Scale Phenomena

The solidification of melt during the casting process can be described in different levels of scales/details. Since most of the conventional casting simulations are carried out at macro-scale, many of the conventional numerical solidification modules are based on the macro-kinematics of material behavior. For these macro-scale simulations, there are different approaches for defining the phase transition. Most common is the application of the Scheil-Gulliver equation [18, 19] where a nonequilibrium approach is formulated for solute redistribution during solidification (no solute would diffuse back into the solid). The solid fraction in this method is defined by local current, solidus, and liquidus temperatures as well as the partition coefficient. The Scheil-Gulliver equation in its simple and standard form can be written as [18, 19](10)CL=C0flk-1CS=kC01-fSk-1where CL,CS, and C0 are liquid, solid, and initial compositions and fS,fS, and k are solid and liquid mass fractions and coefficient k=CS/CL (which can be determined from phase diagram), respectively. The formulation is valid for a local equilibrium state between the solid and liquid interface and for moderate-to-low cooling rates. The numerical methodology for the simulation of solidification front during casting process where phase changes and different phases exist (liquid, mushy, and solids) can be challenging. Different numerical methodologies have been developed [11] to handle the phenomena within numerical domains, namely, conduction-based, conduction-convection based, and enthalpy-based methods. In the conduction-based method, the energy transfer during solidification is calculated based on conduction heat transfer and to account for the convection heat transfer effects (melt flow), an artificially high fluid conductivity is defined. While in the conduction-convection method the melt convection is taken into account explicitly with mass transfer and thermal convections are formulated using governing differential equations for the melt flow. The third alternative way is to use the enthalpy concept. Enthalpy is a thermodynamic quantity defining total heat energy of a system which can be combined with other numerical features to account for change of phase during solidification. In the mainstream industrial software tools for casting simulations, the enthalpy method is generally combined with Volume of Fluid (VOF) and domain porosity features (to damp the fluid momentum when it turns into solid). The VOF method can trace the solidification front within the numerical domain while the domain porosity would enable the smooth change of phase from liquid to solid by damping the momentum equation using momentum sink terms. The introduction of momentum sink terms in the energy equation might introduce more numerical challenges since the convergence of the system becomes more difficult. Consider the mass continuity and momentum equation for the melt flow inside the mold as [4](11)∇.U=0,∂ρU∂t+∇.ρU⊗U=-∇p+∇.μ∇Uwhere ρ,μ,p, and U are fluid density, viscosity, pressure, and velocity, respectively. Momentum sink term is traditionally inserted into melt momentum equation (for temperatures between liquidus and solidus) to damp the melt velocity due to solidification phenomena. A simple and practical approach for momentum damping can be defined using the Kozeny, Carman approach, and its modification by Blake [20–23]. Let us assume the basic equation for fluid flow in resistive media as(12)∆pL=α1-kβk3μvd2where ∆p/L, k, μ, d, and v are pressure gradient, domain porosity (liquid volume fraction), melt viscosity, primary channel width, and melt superficial velocity. Variables α, β can be defined for material processes during solidification based on the alloy and process properties. Two different types of momentum sink terms can be defined for the melt solidification, namely, viscous and quadratic terms. The basic viscous sink term for fluid flow in a porous media can be defined as [22](13)LV=μKUwhereK is permeability of the resistive domain. Voller et al. have modified [24–26] the fluid viscous loss model and the permeability of the domain as(14)Lq=-μK.u-uswhere us is the solidified slab speed (can be assumed as casting speed) andK is permeability which can be calculated based on the resistance against the melt flow inside the dendrite channels during solidification. Ignoring the remelting and the solute concentration effects (temperature exchange at interface, etc.), Gu et al. [27] have proposed the following equation for the permeability:(15)K=Cfl3fs2λ12where fl,fs, and λ1 are volume fractions of liquid and solid phases (at the solidification front) and primary dendrite arm spacing. The constant C can be assumed as C=6×10−4 [27]. The secondary arm spacing and its resistance effect has not been explicitly defined in (15). The broader equation for estimation of permeability has been derived by author using the combination of physical and analytical approach as [28](16)CLV=C∗λ21-fTαfTβfT=C1∗1-T-Ts∆Tγ∆T=TL-Tswhere C, C1 are constant coefficients, λ is the primary dendrite arm spacing, andα,β, γ are power coefficient which can be calibrated for the alloy processes (for aluminum alloysα= 3, β=2, and γ =1). Coefficients C and C1 can be assumed as C= 0.0006 (similar to Gu et. al) and C1=1 - 1.2 [28]. A more detailed version of the simplified loss term can also be developed [28] using primary and secondary dendrite arm spacing and melt isolation phenomena. Similarly, based on physical concepts of energy damping and a mathematical formulation of melt flow kinetic energy, a quadratic sink term can also be derived as [28](17)CLQ=C.β.1-T-TETE-TSαTE<T<TLα=1β+C.1-β.1-T-TSTE-TSαTS<T<TEα=2where TE is a temperature which denotes the starting point for dendrites bridging phase during solidification. For aluminum alloys, the simple assumption can be that TE=TL-TS/2 and β=γ.TE-TS/TL-TS where γ is a coefficient which depends on material microstructure development during solidification (when growth and progression of dendrites starts to create significant resistance against the melt flow near the interface). For aluminum alloys where equiaxed grain growth are perceived, it can be assumed as γ = 0.6. The coefficientC can be calculated/measured and calibrated from the total momentum of the melt within the mold prior to the solidification phenomena. More comprehensive discussions about these mathematical derivations of momentum sink terms can be found in [28].


## 3.3. Simulation of Cooling and Heat Transfer

Using water spray cooling in light metal alloy productions (e.g., continuous and semicontinuous casting processes) is popular and the numerical modelling of these processes is essential in developing a deep understanding of the thermal evolution and temperature fields (including thermal stresses) during industrial processes. For the casting process, the melt material flowing into the mold can be considered as a viscous liquid which confines in a tiny solidified shell formed during initial cooling (during convective air cooling). In many casting processes, the cast billet is then sprayed (or submerged) with a water flow to extract more thermal energy from inside the billet. During this secondary water cooling, the billet undergoes large temperature fluctuations, especially near its surface. Here it is essential to investigate and control the thermal event to minimize the occurrence of defects on the surface or subsurface (hot tearing, cold cracking, etc.). These thermal fluctuations arise due to primary air cooling at the top (thermal convection and radiation), water spray cooling (including shrinkage and thermal conduction throughout the body of the billet), and possible secondary air cooling.Numerical modelling/estimation of water spray cooling effects is one of the most challenging parts of the cooling system design. It would give rise to the estimation of thermal evolution and thermal fluctuations including its associated stresses. Many thermal evolution studies have been carried out in last two decades to develop an estimation technique for the dynamic HTC on the cast billet surface considering different water boiling regimes (nucleate, unstable, and stable filming regimes) and also Leidenfrost effects. One of the simplest models presented in the literature to start modelling of the water spray systems is to define a spray wet area as [29] (18)A=fx,y,zQVswhere fx,y,z, Vs, and Q are coordinate distribution function, surface flow rate, and water flow rate, respectively. Different parameters would affect the thermal energy dissipation during water spray cooling, namely, temperature difference between billet surface and water, water flow rate, angle of impingement, water quality, surface roughness, and some other minor parameters. The boiling regime of the water on the surface of the billet at high temperature has profound effects on the heat dissipation and the complex interface between alloy/steam/water would determine the real cooling HTC. Many researchers have tried to carry out the inverse heat transfer calculation using available experimental data and fitting techniques [30–32]. However, in recent years the combination of empirical and numerical techniques has been used to model the surface boiling regimes and HTC calculations [29]. To start with the physical-mathematical study of HTC calculation, let us start with the simplified governing equations for the water spray cooling system which can be written in a simple form assuming the incompressibility as [31](19)∇.U=0∂U∂t=-∇Pρ+1Re∇2Uwhere U,P,ρ, and Re are velocity vector, pressure, density, and Reynolds number,respectively. After some mathematical manipulation, for the case of continuous liquid phase (water) and dispersive gases (vapor and air) using the three-phase continuity equation, the phasekcan be written as [33](20)∂αk∂t+∇αkU→k=Γki-Γikρkand the total interfacial force in the momentum equation can be elaborated to find the bubble generation and movement relative to the billet surface (see Figure 5). For the water spray cooling system during casting process, the total heat flux can be defined as [33](21)qT=q1+q2+q3where q1,q2, and q3 are single-phase convective, evaporation, and quenching heat flux densities, respectively. Different models for the range of bubbling schemes in water flow boiling regimes have been developed in the past decades. In many of these models the condition that a vapor bubble would grow if the surrounding liquid temperature is equal to or more than the saturation temperature of the vapor inside the bubble is satisfied. In the research work herein, attempt has been made to formulate new mathematical models for the HTC at different boiling regimes during water spray (and submerge) cooling. For the initial air cooling during casting, conventional HTC approximation technique which includes convection and radiation terms can be adopted. If the temperature of billet surface is assumed asTb and air bulk (free-stream) temperature isTa, the thermal resistance (or inversely conductivity) at the interfaceAT can be defined as [34](22)R=1hr+hcATwherehr andhc are, respectively, the radiation and convection HTCs, calculated as follows:(23)hr=εσT4per unit areawhere σ is a constant (Stefan–Boltzmann constant) which can be assumed as 5.672×10−8 (W/m2 K4),T is the absolute temperature in Kelvin and ε is the billet surface emissivity. The convective HTC (hc) for air cooling can also be calculated using [34](24)hc=kNuλwhereNu is Nusselt number, k is the air thermal conductivity, and λ is the characteristic length. The Nusselt number (the ratio of conductive thermal resistance to the convective thermal resistance of the air) can be calculated as (25)Nu=CGrPrnGr=gβTb-Taν2λ3Pr=νkcpwhereCandnare constants andPr is the Prandtl number (defined as the ratio of momentum diffusivity to thermal diffusivity). Gr is the Grashof number (defined as the ratio of the buoyancy to viscous force acting on air) for vertical flat plates, β is the coefficient of thermal expansion (approximately 1/Ta, for air), ν is the kinematic viscosity, and g is gravitational acceleration. Figure 5 shows the schematic bubble generation and separation event during surface cooling and temperature dependent radiation HTC for aluminum alloy.Figure 5
(a) Bubble generation and separation mechanisms from hot surface and (b) radiation of heat from aluminum billet at different temperatures.
(a) (b)The water cooling HTC values change rapidly during the casting startup condition where continuous cooling regime adds to the modelling complications. The HTCs generally rise rapidly with the billet surface temperature until a critical temperature is reached above which unstable and film boiling regimes, with their lower associated HTCs’ starts. The estimation of critical heat flux (CHF) temperature point and its associated heat flux can help to understand and optimize the billet surface cooling rates and thermal-mechanical stress/strain conditions. The variation of HTC against the billet surface temperature can be divided into four distinctive zones, namely,(i)
below water saturation temperature,(ii)
between water saturation and CHF temperature (nucleate boiling regime),(iii)
from CHF temperature to Leidenfrost temperature (unstable boiling regime),(iv)
above Leidenfrost temperature (film boiling regime).Some aspects of these boiling regimes have already been investigated [34] where mathematical formulation and also empirical estimation of heat transfer have been developed. For the work herein, different sets of mathematical works and experimental verification for the heat transfer during cooling process have been carried out to optimize the HTC during casting [35]. For the preboiling regime a simplified approach can be employed (the preboiling has minor thermal energy exchange effect) and HTC can be defined as [35](26)HTCpre=αQLbβ7e2∗Tb-1.67e3whereQ,Lb, Tb are water flow rate, billet perimeter, and billet surface temperature, respectively. The parametersα, βthat can be defined using the impingement conditions are billet surface quality (for aluminum billetsα, β can be assumed as α = 1.35 and β = 0.3). For the more challenging nucleate boiling regimes, a more detailed formulation has been sought [35] and HTC can be defined as(27)HTCnu=μHfgTsatcfαCpTbHfgασμCpKgρw-ρvβwhere μ, Hfg, cf, Cp, K, g, ρw, ρv, and σ are coolant dynamic viscosity, latent heat of evaporation, flow heat parameter, specific heat, conductivity, gravitational acceleration, density, and vapor density, respectively. Tsat and Tb are defined as coolant saturation and billet surface temperature values where spatiotemporal distribution of temperatures on billet surface can be calculated using numerical domain. α and β are again cooling system parameters (these parameters can be assumed as α =0.8 and β =-0.515 for aluminum alloys).The third boiling regime during cooling process is the unstable film boiling where the interphase boundary is unstable and vapor film can be broken periodically for water to make contact with billet surface. The estimation of HTC for unstable boiling zone has been a challenge for both scientist and engineers trying to overcome the full dynamic and oscillatory nature of the boiling event. As the billet surface temperature decreases during cooling process, the duration of vapor near the surface reduces, while the number and duration of liquid-wall contact increase. Such sporadic on/off contacts of liquid water with the billet surface result in an abrupt discharge of vapor from the surface, causing splashing and turbulence (and driving away vapor from the wall). The HTC for this zone can be estimated using the fraction of billet surface wetted by liquid water during the cooling process. The proposed empirical formulation for this zone can be defined as [35] (28)q=qwAc+qv1-AcHTCtra=HTCnu1-Tb-TCHFTbγwhere TCHF is water critical heat flux temperature (can be assumed as TCHF=185°C for aluminum alloys). The final boiling regime during cooling process is the stable film boiling zone where a vapor blanket would be covering the surface of the billet. Generally, the estimated HTC values by the water spray is accredited to the intensified forced convection in the liquid film caused by the large flow momentum in the gravitational direction and consequent improvements in nucleate boiling. It also promotes the transient conduction which occurs on the heated surface due to the quick refresh of the liquid film. The proposed HTC estimation formula for cooling zone IV (film boiling regime) during cooling process can be proposed as [35](29)
H
T
C
f
=
H
T
C
t
r
a
1
-
T
b
-
T
L
e
T
b
H
1
+
H
2
ƙ
(30)
H
1
=
γ
T
b
×
k
w
3
g
ρ
w
H
f
g
ρ
w
-
ρ
v
1
+
c
p
T
b
/
2
H
f
g
μ
w
T
b
σ
/
g
ρ
w
-
ρ
v
β
α
(31)
H
2
=
5.4
×
10
4
T
b
1
+
0.055
T
b
v
s
κwhere TLe,Tb,μw,Hfg,Cp,σ,ρw,Kwρv, are water Leidenfrost temperature (≈400°C), billet surface temperatures, dynamic viscosity, latent heat of evaporation, specific heat, surface tension, density, and thermal conductivity, respectively. ρv is the vapor density and α,β,γ,κ, and ƙ are film cooling system parameters (these parameters can be assumed as α =0.4, β =0.5, γ=0.5, κ=0.6, and ƙ=0.15 for aluminum alloys). The film HTC calculation is valid for surface temperature (Tb) between 400°C and 550°C. More detailed discussion about numerical modelling of boiling and bubble generation can be found in [35].


## 4. Case Studies

The first part of the casting simulation technique is the thermal-fluid mold filling and its material thermal evolution using CFD technique. CFD technique has extensively been used in the last couple of decades for fluid and fluid-structure engineering problems. Understanding the motion of fluids, its thermal conduction, convection/radiation, and also its interactions with solids is crucial for the design of engineering systems in many branches of engineering. To develop a new numerical industrial casting simulation tool, two comprehensive sets of parallel simulation-experimental studies have been planned, namely, thermal-fluid CFD simulation for mold filling/solidification and evolving thermal-mechanical simulation for stress/strain states and cracking. The most challenging parts of thermal-fluid CFD applications for the industrial process engineering and in particular casting processes are the numerical solidification and cooling simulations.


### 4.1. Thermal-Fluid Casting Simulation

As described in the previous sections, for the proposed casting simulation tool, it is essential to model the change of phase from fluid to mushy and solid state under the thermal gradient caused by cooling process. In the research work herein, a new solidification module based on momentum sink terms described in the previous section has been implemented and coded in-house. The thermal gradient and cooling effects during solidification phenomena have also been implemented using an in-house code which incorporate the developed HTC functions for different stages of water boiling regimes.For the experimental studies, two sets of casting trials have been selected, namely, horizontal and vertical direct chill casting procedures. These casting trials are performed on semi-industrial scale where aluminum billets are cast using direct chill method (air and water cooling). Figure6 shows the vertical direct chill casting trials in which rectangular and round cross-section billets were cast. Two sets of parallel numerical simulations have also been carried out to examine and calibrate the developed mathematical models for solidification and HTC estimations. For the vertical experimental casting trial, two series of Thermo-Couples (TC) are installed across the cross-section at different elevation to measure temperature for both “start of cast” and “during cast” conditions (thermal transient and steady state cases). These two TC series are installed on an especially designed supporting “frame system” (see Figure 6) which is fixed on the bottom tray during casting.Figure 6
(a) Vertical direct chill casting trial; (b) plan and side view of TC arrangement for temperature measurement; (c) thermal-fluid simulation; (d) measured temperature on billet surface; and (e) measured temperature in mold.
(a) (b) (c) (d) (e)A second experimental-simulation case study has also been carried out using a horizontal continuous casting technology to study the continuous casting process of lightweight alloys (see Figure7). The fluid-thermal conditions of the process have been simulated with user-defined modules developed for the CFD solver (2D axisymmetric and 3D half models). The casting process cooling rates and change of phase (solidification phenomena) within the mold have also been modelled using developed mathematical functions. Figure 7 shows casting of round billet using the horizontal casting machine for light weight alloys and its numerical 2D axisymmetric and 3D half models.Figure 7
(a) Horizontal direct chill casting trial, (b) 2D axisymmetric CFD model, and (c) temperature results for CFD half model.
(a) (b) (c)


### 4.2. Thermal-Mechanical Casting Simulation

The industrial casting process techniques for metals and alloys, including continuous and semicontinuous processes have intrinsically a dynamic nature, where billets are produced by solidification of melt material at “casting-speed” rate. The optimization of these casting technologies and the issues of product final mechanical properties, production speed, and energy/cost saving have been major topics for many years. The process of developing a numerical dynamic and flexible tool to model the thermal-mechanical evolution of the billets during a continuous casting process is among the most important ambitions for researchers and engineers in this field. As presented in the previous sections, following vast expansion of computational resources and also introduction of efficient parallelization schemes, a new concept of evolving domain has been developed herein to address the transient casting conditions. The combination of multidomain concept, energy input/output balancing (for melt thermal energy content), and advanced history mapping (using overlapping scheme) has been implemented [11].The dynamic mesh procedure with its system-matrix reassembly/appending capability would enable the development of a fully transient scheme to introduce/append new subdomains (melt entry) based on material and energy balance concepts. The thermal-mechanical casting simulation framework herein has integrated the following modules for industrial applications:(i)
evolving domain module which can handle different engineering dynamic processes with internal interfacing/moving parts using developed cooling, solidification, overlapping, and mapping techniques;(ii)
a comprehensive set of experimental mechanical tests at room and elevated temperatures for calibration of numerical damage models (including tension tests, etc.);(iii)
a multiresolution zoning technique [11] for simultaneous modelling of phenomena at different length scale during solidification (macro-meso simulation of grain growth using Lattice-Boltzmann and Cellular Automata which have not been presented in this paper; please refer to [36])Figure8 shows the experimental-simulation trials for thermal-mechanical numerical simulation using tension testing machine.Figure 8
(a) Aluminum sample geometry and its fracture simulations. (b) Tension test machine. (c) Measured and simulated true stress-strain curves for aluminum samples.
(a) (b) (c)


## 4.1. Thermal-Fluid Casting Simulation

As described in the previous sections, for the proposed casting simulation tool, it is essential to model the change of phase from fluid to mushy and solid state under the thermal gradient caused by cooling process. In the research work herein, a new solidification module based on momentum sink terms described in the previous section has been implemented and coded in-house. The thermal gradient and cooling effects during solidification phenomena have also been implemented using an in-house code which incorporate the developed HTC functions for different stages of water boiling regimes.For the experimental studies, two sets of casting trials have been selected, namely, horizontal and vertical direct chill casting procedures. These casting trials are performed on semi-industrial scale where aluminum billets are cast using direct chill method (air and water cooling). Figure6 shows the vertical direct chill casting trials in which rectangular and round cross-section billets were cast. Two sets of parallel numerical simulations have also been carried out to examine and calibrate the developed mathematical models for solidification and HTC estimations. For the vertical experimental casting trial, two series of Thermo-Couples (TC) are installed across the cross-section at different elevation to measure temperature for both “start of cast” and “during cast” conditions (thermal transient and steady state cases). These two TC series are installed on an especially designed supporting “frame system” (see Figure 6) which is fixed on the bottom tray during casting.Figure 6
(a) Vertical direct chill casting trial; (b) plan and side view of TC arrangement for temperature measurement; (c) thermal-fluid simulation; (d) measured temperature on billet surface; and (e) measured temperature in mold.
(a) (b) (c) (d) (e)A second experimental-simulation case study has also been carried out using a horizontal continuous casting technology to study the continuous casting process of lightweight alloys (see Figure7). The fluid-thermal conditions of the process have been simulated with user-defined modules developed for the CFD solver (2D axisymmetric and 3D half models). The casting process cooling rates and change of phase (solidification phenomena) within the mold have also been modelled using developed mathematical functions. Figure 7 shows casting of round billet using the horizontal casting machine for light weight alloys and its numerical 2D axisymmetric and 3D half models.Figure 7
(a) Horizontal direct chill casting trial, (b) 2D axisymmetric CFD model, and (c) temperature results for CFD half model.
(a) (b) (c)


## 4.2. Thermal-Mechanical Casting Simulation

The industrial casting process techniques for metals and alloys, including continuous and semicontinuous processes have intrinsically a dynamic nature, where billets are produced by solidification of melt material at “casting-speed” rate. The optimization of these casting technologies and the issues of product final mechanical properties, production speed, and energy/cost saving have been major topics for many years. The process of developing a numerical dynamic and flexible tool to model the thermal-mechanical evolution of the billets during a continuous casting process is among the most important ambitions for researchers and engineers in this field. As presented in the previous sections, following vast expansion of computational resources and also introduction of efficient parallelization schemes, a new concept of evolving domain has been developed herein to address the transient casting conditions. The combination of multidomain concept, energy input/output balancing (for melt thermal energy content), and advanced history mapping (using overlapping scheme) has been implemented [11].The dynamic mesh procedure with its system-matrix reassembly/appending capability would enable the development of a fully transient scheme to introduce/append new subdomains (melt entry) based on material and energy balance concepts. The thermal-mechanical casting simulation framework herein has integrated the following modules for industrial applications:(i)
evolving domain module which can handle different engineering dynamic processes with internal interfacing/moving parts using developed cooling, solidification, overlapping, and mapping techniques;(ii)
a comprehensive set of experimental mechanical tests at room and elevated temperatures for calibration of numerical damage models (including tension tests, etc.);(iii)
a multiresolution zoning technique [11] for simultaneous modelling of phenomena at different length scale during solidification (macro-meso simulation of grain growth using Lattice-Boltzmann and Cellular Automata which have not been presented in this paper; please refer to [36])Figure8 shows the experimental-simulation trials for thermal-mechanical numerical simulation using tension testing machine.Figure 8
(a) Aluminum sample geometry and its fracture simulations. (b) Tension test machine. (c) Measured and simulated true stress-strain curves for aluminum samples.
(a) (b) (c)


## 5. Discussion and Concluding Remarks

The speed and the quality of the industrial production process are depending on sound design and optimisation procedures which can be numerically simulated using rigorous scientific/engineering virtual tools. For many industrial material processes and their numerical simulation tools it is required to establish a framework for continuous innovation and sound predictive apparatuses for confronting accumulative challenges and competitive market. The combination of state-of-art innovative numerical approaches (with flexibility, agility, and predictive power) and thorough calibration and verification schemes to develop and explore what technical solutions are viable can be promoted for future of process engineering. The integrated framework for casting simulation developed herein consists of practical evolution scheme, pioneering analytical approaches (for solidification and cooling), and multiresolution grid technology (only briefly discussed in this paper; please refer to [36]) which would reflect the essence of new engineering and industrial virtual tools. The generic simulation framework would help to improve industrial material casting processes and to tailor future research works on optimized process parameters for various industrial applications using its developed modules. A series of carefully designed experimental trials have also been planned and carried out which helped to verify and calibrate numerical models. The idea herein was to develop a framework for an efficient numerical technique and integrated it into mainstream industrial simulation tools to help(i)
optimise and tailor casting production processes for reduction of costs and increase part quality and properties based on integration of practical numerical predictive modules;(ii)
develop innovative and reduced-risk of design/production strategies based on affordable advanced material process models;(iii)
amplify the use of computer-based virtual modelling and simulation systems for material process industries by developing affordable advanced numerical tool;(iv)
improve the confidence and reliance of industrial use of numerical tools and their predictive powersAlthough conventional numerical tools and software have been used for the simulation of industrial processes with relative success, their limited applications, large computational CPU timing, accuracy issues, and their lack of modelling for vital lower-scale phenomena are well understood. In addition, the single-domain and single-solver simulation framework of software would limit the application of them for multiphysical and multiscale process simulation.In the first part of this paper, the analytical and theoretical basis of the evolving domain technique and its overlapping and mapping methods were briefly presented. It was shown that these methods can be combined with new mathematical and analytical models to enable the making of accurate prediction for material evolution during industrial casting processes. The approach is new, as the conventional industrial numerical tools appear to lack accuracy and modelling power for underlying material phenomena. The parallel experimental-simulation trials have also been shown in the second part of the paper where the calibration and benchmarking of the numerical method have been studied. It can be concluded that the numerical simulation framework founded on combination of mathematical-analytical schemes with conventional solver technology can be implemented and employed in an upfront manner. The solidification and cooling models presented herein allow for the continuous transition of phases from melt to solid state.As far as numerical stability and efficiency are concerned, it can be assessed that the assemblage of global matrices has to be repeated for evolving numerical domain (starting with small matrices and appending at discrete steps) in comparison with the single assemblage procedure for the conventional methods (with large full matrices). However, further numerical benchmarking of these two methods has proved that the proposed method outperforms the conventional method for industrial-scale process simulations (results not reported in this paper; it would a subject of next submission). Furthermore, the proposed method is capable of dealing with material evolution at lower scales (e.g., solidification phenomena) through multiple overlapping and mapping schemes in a straightforward manner without the need for sequential simulation steps.



---

*Source: 6164092-2019-05-08.xml*