Abstract
This article demonstrates implementation of immersed boundary method in continuous casting simulation involving boundary movement. In this methodology, the immersed boundary method is coupled with the second-order accurate finite difference solution of unsteady three-dimensional heat conduction equation. The moving molten metal front is modelled using the immersed boundary method in a Cartesian mesh framework that provides simplicity in its implementation and reduces the computational time as compared to the adaptive mesh solutions. A parallel programming paradigm using message passing interface has been implemented to obtain enhanced computational efficiency. This study has focused on capturing moving boundary during continuous casting and predicts the temperature distribution and shell thickness under different cooling ambiences and casting function. Good agreements with published data and correlations are obtained through numerical analysis. Mould-region shell thickness agrees well with Chipman–Fondersmith correlations. A new correlation has been further proposed for the delay constant at different heat extraction rates. The effects of key parameters like casting speed, convection and radiation from the continuous casting are also quantified in attempt to avail the data for optimal design of continuous caster.
Keywords
Introduction
Thermal transport phenomenon plays an important role in many manufacturing processes, and several experimental/numerical studies have been aimed to understand the thermal physics involved in these processes.1,2 Continuous casting is one of the heavily practiced metallurgical processes where thermal transport phenomenon plays a significant role. It has an edge over other conventional casting practices due to several factors like simplicity, higher productivity, eco-friendliness, lower operating cost and wide applicabilities in product manufacturing at industrial scales. In this process, superheated liquid metal is fed into water-cooled mould along the gravity from which it is withdrawn continuously in order to have partially solidified billet, bloom or slab, which are further processed in rolling mills. Electromagnetic stirrer is also present to facilitate stirring process. It helps in achieving better homogeneity and properties of casting product.2–6
During continuous casting process, superheated molten metal flows from the tudnish to the primary zone through meniscus and then to the secondary zone in order to obtain fully solidified state under controlled cooling ambiences. Superheat, as well as the latent heat, is extracted from the liquid regime while sensible heat is extracted from the solid regime. The primary section in the casting consists of water-cooled mould. Sets of nozzle are provided in the secondary section to spray water-mist on the upcoming strand. After this, there is a subzone for radiation cooling. These three zones appear in longitudinal succession along the length of the casting strand. Electromagnetic stirrers are present at different locations along the length of the casting. Depending upon the location, the electromagnetic stirrers are classified as mould, strand or end-electromagnetic stirrers. It also helps in maintaining the hyperbolic shape of solidification front due to the presence of magnetic field.4,7 It provides intensive flow of molten metal inside the mould and thus assists in better exchange of heat between the surface and the interior of the casting product. This is due to the fact that the stirring is done using an inductor, which produces magnetic field and helps in setting eddy currents across the metal pool. In effect, Lorentz force is obtained, which helps in the movement of the liquid metal to the core of the strand. 5 The stirrer prevents sticking of casting strand to the mould surface (which may lead to breakout). Moffatt 4 showed that application of electromagnetic stirrer results in better surface finish of the casting product. Furthermore, it also helps in controlling the growth of dendritic microstructures during solidification. 8 Nozzles in the secondary zone spray air–water mist on the surface of the out-coming casting strand. In the radiation zone, ambient still air is the sole cooling medium. At the starting of the process, a dummy bar is horizontally held at the end to restrain molten metal from spilling. It is removed once complete solidification is achieved in its vicinity. A general schematic diagram shows the arrangements of these different zones along a continuous caster (Figure 1).

Schematic diagram of continuous casting process.
The casting process parameters are designed to obtain a sufficient solidified shell thickness at the primary zone, which should be able to deter the bursting phenomenon (a phenomenon in which the liquid molten metal comes out of surface below the exit of the mould). Occurrence of aforesaid phenomenon halts the continuity, and it becomes inevitable to restart the process from the beginning. Restarting the process is time consuming and results in wastage of molten raw materials, and thus increases the cost of production. Hence, the physics of metal flow in the primary zone has a very important role in the success of continuous casting process as well as in the quality of output products.
In the field of continuous casting, over several years, attempts have been made to develop this understanding in order to avoid surface defects, turnarounds and other defects. Rates of heat extraction from different sections of the cooling zones are of interest due to their influences in directional solidification and effects in the grain structures of the final product. Heat transfer rates can be obtained through analysis of the temperature profiles obtained for different arrangements of cooling section, dimension of cross section, rate of molten metal flow and intensive flow of metal due to electromagnetic stirrer. With the help of a priori information about the cooling parameters, the height of the mould can be determined which can sustain bursting (sufficient solidified shell thickness to balance ferrostatic pressure posed by molten metal). The other important parameters are length of secondary zone, radiation zone, sets of nozzle required in cooling zone, flow rate of water flow from nozzles and rate of molten metal flow in mould. Better understanding of the thermal physics of continuous casting process can help to build a pool of information on the effects of these parameters and finally to specify the optimal design of the process.
Several numerical models have been developed with respect to one-dimensional (1D),9–11 two-dimensional (2D)12–15 and three-dimensional (3D)16,17 steady and transient assumptions in order to predict the temperature distribution, heat fluxes and solidification patterns. Researchers have used the finite element method, 16 inverse finite element method, 18 finite difference method,10,12,13,15,16,19 finite volume method 14 and boundary element method 20 to solve numerically. The first few solidification models were based on 1D energy equation. Mizikar 9 utilized 1D model to predict the thermal field and shell growth. Brimacombe 10 predicted shell thickness and initial crack location using the finite difference scheme and further studied effects of reheating in continuous casting. Lally et al. 12 used 2D slices to model 3D steady-state temperature fields during continuous casting. Meng and Thomas 13 performed heat transfer simulation in the mould as well as in the spray zone. Their model coupled 1D transient calculation of heat conduction within the solidifying steel shell with 2D steady-state heat conduction through the mould wall. They also considered air gap between the mould and solidified shell due to shrinkage of casting strand, effect of oscillation marks and mass and momentum balances at phase interface. Boehmer et al. 19 studied the thermo-mechanical behaviour in continuous casting at high temperature and its effect on casting product using the finite difference method. They inferred that if effect of casting parameters is known a priori, then optimal casting condition could be addressed and hence able to avoid thermal cracking. Li and Thomas 11 proposed transient heat transfer model for single roll continuous strip casting process. Heat transfer through rotating wheel was also considered in their work. Tieu and Kim 16 performed 3D simulation to get thermal field and determined shell thickness for multi-zonal cooling ambiences. Janik and Dyja 17 used 3D energy equation to solve temperature field using the finite element method (ANSYS). Shell thickness and heat flux inside the mould were determined from the temperature distribution. Alizedeh et al. 14 also utilized temperature field to predict shell thickness, interfacial gap between shell and mould and liquid pool positioning. They used a coordinate transformation technique to model the billet movement.
As continuous casting involves change of phases, it is important to track the interface between the molten and solidified phase. Siegel 21 did extensive work on studying the interface behaviour during continuous casting and showed that its shape depends upon a single dimensionless parameter, which is function of casting velocity, the boundary conditions at the surface (cooling conditions) and width of the ingot. Malamataris and Papanastasiou 18 proposed an empirical correlation for the shape of solidification front in terms of process parameters. They extended inverse isotherm finite element method to consider distorted isotherms. They also analysed high-speed continuous casting utilizing energy equation. Majchrzak 20 utilized boundary element method to trace the solidification front and obtain thermal field during continuous casting in order to predict temperature pattern across the billet. Das and Sarkar 22 proposed a prediction model to predict the thermal profile during continuous casting using a coordinate transfer technique in a non-orthogonal control volume framework. This model involves tracking of the solid-liquid interface. A meshless method is also used to study the thermo-physics of pure metals during continuous casting method. 23 These approaches demand high computational time and complexity as it involves grid generation/mesh deformation at each time step to model the moving boundary.
Review of the literature (referred earlier in this section) reveals that for the study of thermal physics of continuous casting, different analysis approaches have been utilized. In broad way, these approaches could be categorized as follows:
On the basis of how energy equation is implemented in the considered domain: Single region: ingle energy transport equation is considered for the entire casting irrespective of phase; Multi-phase region: Energy transport equation is separately specified for each phase. Boundary conditions are imposed at interface through heat balance equation.
On the basis, how moving boundary problem (lower plate movement) is formulated: By transforming the governing equation in a moving reference frame; By simulating the boundary movement using dynamic grid generation.
Earlier works established that numerical simulations could be an effective tool to augment the detailed knowledge of different manufacturing processes (e.g. die-casting, investment casting, welding and friction stir welding (FSW)).2,24,25 Recently, Soe 26 demonstrated that commercial softwares can be useful for accurate simulations of manufacturing processes and identify the important input parameters for optimization studies. Song and Kovacevic 25 also showed that for a manufacturing process like FSW, 3D transient heat transfer models can well predict the temperature variation during welding. Furthermore, numerical simulation has an edge over experimental approaches as it saves the material, time, money, energy and helps in building a pool of data for different materials along with the different process parameters. But, most of time, these manufacturing processes are involved with complex and moving geometries. Hence, they pose substantial challenges in boundary treatment while solving by traditional body conformal grid approaches. In contrast, immersed boundary method (IBM) (a non-conforming boundary approach) 27 has an edge over the boundary conforming grid. Here, as it offers much simplicity and computational efficiency in handling the moving boundary cases. In addition, immersed boundary (IB) also favours parallel paradigm without an increase in complexity.
IB method was first coined by Peskin, 27 who used this to simulate cardio-vascular flow problems. Most important feature of this method is that the entire simulation can be performed in the background of a Cartesian grid, which does not aligned with the complex geometry of the moving/fixed boundary, while the replication of boundary condition is achieved by imposing appropriate forcing term in the governing equation. This method needs an efficient algorithm to find the near-boundary mesh points where the interpolated values are to be specified. Simplicity and computational efficiency of IB methods have been discussed in a number of earlier published reports.28–30 The advantages arise mainly due to the following facts: (a) costly body-fitted mesh generation can be avoided for complex geometry; (b) in body-fitted mesh approach using non-rectilinear meshes, calculation of metric quantities and flux variables are computationally involved, which is absolutely not needed in IB and (c) for moving geometry, re-meshing and mapping of older solutions to new meshes after each computational time step is not required. Computational efficiencies of IBM schemes were demonstrated by Verzicco et al. 31 They reported that simulation time was much less in IBM than a coarser body-fitted mesh simulation 32 for flow in a piston-cylinder arrangement with piston movement. Tyagi and Acharya 33 utilized IB method for large eddy simulation (LES) of turbulent flow in complex geometries. Their work also demonstrated different strategies for implementation of the forcing function. Mittal and Iaccarino 34 extensively studied different methods of IB implementation for immersed solid boundaries and shows its applicability over wide range of flow problems (such as flow past sphere and flapping wing). Heat transfer through complex and moving boundaries has been predicted by several researchers using IB methods. 35 Fadlun et al. 36 utilized IB method to study behaviour of fluid flow past complex geometry and reported the second-order accuracy in their results.
In this article, finite difference method for the 3D transient heat equation along with the IBM27,34 has been used to predict the temperature profile and solidification pattern under multi-zonal cooling ambiences for continuous casting of steel alloy. A single region approach is used, which considered different phases of the steel alloy associated with respective thermal properties. For the mixed phase, the property values are interpolated at the local temperature using the solidus and liquidus properties. The moving boundary is modelled using a second-order accurate IB interpolation in a fixed coordinate system. The proposed methodology is simple to implement and moreover can be easily translated to parallel programming paradigm. This method further claims its applicability due to higher computational efficiency.37,38 Solidification patterns are predicted for multi-cooling ambiences and effects of different casting parameters are studied.
Methodology
Problem formulation
The whole domain is considered as single region. 3D unsteady state heat conduction in Cartesian framework (equation (1)) is considered as governing equation for predicting the temperature profile of casting strand of alloys
Due to pouring and electromagnetic stirring effects, flow near the mould zone shows turbulent behaviour (Re: 104–105) leading to higher heat transfer. However, this article attempts to maintain the simplicity of the equation (1) without incorporating turbulence models into it and thus reduces the costs of computation. In earlier works, the effects of heat transfer augmentation due to turbulence have been mimicked by assigning a higher equivalent value of conductivity at the liquid and mushy zones.9,39,40 The effective thermal conductivity can be four to eight times higher than its thermal conductivity of solid phase.9,10 The idea of using higher equivalent k for liquid is first given by Mizikar 9 to accommodate the effect of fluid motion and turbulence heat transfer. Till date, the higher conductivity modelling has widely been applied by many authors who envisaged the continuous casting phenomena. Saraswat et al. 41 also showed that for a steel alloy, an effective higher thermal conductivity 5–15 times greater than the solid thermal conductivity predicts the solidification thickness within a range of 10% error as compared to k-epsilon model. However, it was established by Choudhary and Mazumdar 42 that the predicted shell thickness depends on the values of equivalent conductivity (k). They also showed that better agreement with experimental result is obtained when a conjugate heat and flow simulation is performed. It is to be noted that all equivalent k approaches have completely ignored the convective effect due to motion of the bulk of molten/solidified metal during the entire continuous casting process. In this case, a partial k approach is used, where the equivalent k is considered only for turbulence due to electromagnetic stirrer (predominant in the mould zone liquid and mushy phase) while heat transfer due to bulk convection of the metal (both liquid and solidified) is taken care by a convective term. Our proposed model has reduced the dependency on higher equivalent k (choosing an equivalent k value in the range of 5–15 times of solid phase conductivity) as used in previous studies.9,10,17,39,43 A lower value of equivalent k (∼2–3 times of solid phase conductivity) has been considered in order to model the augmented heat transfer rates due to pouring, stirring and movement effects along with the convective heat transfer term (second-term in the right-hand side (RHS) of equation (1)). This formulation avoided solving expensive fluid momentum equations for the turbulent electromagnetic stirred mould region and reduced computational expenses.
Latent heat
The latent heat of solidification can be considered by adding an enthalpy term in the governing equation (1) or by defining a specific scheme for increase in temperature for the mushy zone. The released fraction of latent heat in this regime is a function of the solid fraction. 44 Linear relationship can be used to calculate the equivalent specific heat during mushy zone. The equivalent specific heat for mushy zone is given by equation (2) as alloys have different liquidus and solidus temperature
Density
Generally, density is a function of temperature but the rate of change of density in relation with temperature for steel alloys is small due to lower values of thermal expansion coefficient suggesting a piece-wise continuous distribution over different phases. However, three distinctive densities have been considered respectively for the liquid, mushy and solid states.
Thermal conductivity
Thermal conductivity is modelled as linear function of temperature. In this case, a convection term in the direction of flow (casting speed, u) has been considered, whereas the effects of turbulent velocity and heat transfer in the horizontal plane are modelled using an equivalent higher thermal conductivity. The effect of turbulence and stirred whirls are dominant in the liquid and mushy phase and negligible in solid phase. Hence, higher value of equivalent thermal conductivity is considered in the near-mould regime where the effect of electromagnetic stirring is dominant.
Moving boundary modelling using IBM
During continuous casting, the molten metal flows as the lower boundary (or lower plate) moves and the casting grows in size. In order to utilize the simplicity of IBM, simulation is performed on the background of a fixed rectilinear mesh. For ease of implementation, the lower boundary is kept at a constant height and the upper boundary is moved in the opposite direction to produce the same effect. Initially the nodes outside the billet do not belong to the physical domain of the casting strand, and hence, they are discarded in the calculation. However, the location of the boundary changes with time as the upper meniscus moves. As a result few of the nodes, which were earlier outside the physical domain, are now inside the casting and therefore considered for calculations. Simultaneously, the near-boundary tagged points are also changed accordingly. The same is represented in Figure 2. The reconsideration of tagged points and solution points are performed at the beginning of each computational time step. The advantage of using IBM can be realized in using a fixed Cartesian mesh as the steps of reconstruction or deformation of the mesh can be avoided, although the volume of casting increases due to boundary movement. Only the treatment of different mesh points changes as per their location with respect to the moving boundary. The points outside the physical domain, that is, above the upper boundary, are discarded in calculation. Temperature field is interpolated at the tagged points to satisfy boundary conditions and at the solved or inner points, the finite difference form of energy equation is solved to predict the temperature.

Growth of casting strand with time and implementation of Immersed boundary method at these instants.
As the upper physical boundary (i.e. meniscus) is moving upwards and the casting is growing in size, the fixed rectilinear mesh framework results into three cases for any node in the geometry at a particular time instant: the node can be inside, on the surface and outside of the billet (Figure 2). So, first, the nodes which are inside the billet (solved nodes) close to boundary of billet (immersed nodes) and on the surface of billet are identified. This identification is done at every time step to capture this position of moving boundary front before solving the temperature field for the next time level. Temperature at the inside nodes is computed by solving the finite difference equations and the temperature at nodes close to boundary front (immersed point) is interpolated from the boundary temperature and inside point’s temperature values. The interpolation function used to reconstruct the solution at these nodes is expressed as equation (3)
where n is the distance of the immersed grid point from the moving surface. Mittal et al. 45 presented a similar approach for solving momentum equations with a moving boundary. The immersed surface is identified and a quadratic interpolation function is used along the normal direction to specify temperature at the immersed nodes. As per Figure 2, at blue nodes (immersed points), the temperature is calculated using the interpolation of red and green nodes. The red nodes are at the actual boundary where boundary condition (T = meniscus temperature) is known to us. Green nodes are the solved nodes, which are solved directly from the finite difference form of energy equation.
Discretization methodology
Finite difference method with Forward-Time Central-Space (FTCS) scheme is used to solve this problem numerically. Explicit time integration scheme is used and a von Neumann criterion is considered for stability. Ghost cell method is used to satisfy the Neumann boundary conditions (at fixed boundaries).
Generally, time increment, Δt, is calculated on the basis of von Neumann stability criteria. The first-order upwind scheme is used for convection term. In addition to these criteria, Δt satisfied the relation (4) arising from the stability of the IBM implementation
Parallelization of the solver
3D transient equation is computationally intensive and needs higher computational time. So in order to reduce the computation time, a domain decomposition technique has been used, which distributes the nodes in several processors and solves each sub-domain independently. Communications among the processors are achieved using message passing interface (MPI). Thus, parallel programming has enabled us to reduce the computational time by utilizing high-performance computing resources with multi-core architecture. Figure 3 shows the current implementation up to eight processors with an almost-linear scalability.

Scalability up to eight processors.
Geometry and boundary conditions
The billet is of rectangular shape and it facilitates two axes of symmetry, so only a quarter of billet has been solved with symmetric boundary conditions at the respective edges in order to reduce the simulation time. Schematic of the boundary conditions representing different cooling ambiences are depicted in Figure 4. Since the meniscus is moving, the boundary conditions at the surfaces of casting strand are also moving along with the meniscus movement to maintain the exact locations of different boundary conditions.

(a) Horizontal section of billet with boundary condition, (b) vertical section of billet with boundary conditions and (c) different cooling section for multi-zonal cooling ambience.
Primary zone (mould)
The heat flux extraction rate from surface of billet is expressed as equation (5)
The value of hm could be determined by an integral heat balance on the mould. It depends upon average heat flux to the mould water and dwell time. The value of hm could also be determined by Hill’s 46 correlation.
Secondary zone
Air/mist cooling
In this zone, a number of nozzles spray water on casting strand with different water flow rate. This is done to achieve controlled cooling. The boundary condition which replicates this phenomenon can be expressed as equation (6)
Heat transfer coefficient, h, can be expressed as equation (7) in relation to water flow rate 47
Radiation zone (air cooling)
In this zone, the casting strand surface is only in contact with still air, which has negligible convection effect. The boundary can be expressed as equation (8)
where qr is govern by Stefan–Boltzmann’s law and expressed as equation (9)
The emissivity of oxidized surface is generally taken as 0.8. 10
Material data
Three different grades of steel alloys have been used for validation cases namely, 0.3 wt% C steel, 42CrMo and 5G2A steel for the set-1, set-2 and set-3, respectively. 15G2A steel alloy is used in remaining analysis (i.e. set-4 to -11). The composition and material properties of the alloys are available in the cited literature.15,17,43,48
Results and discussion
Validation study
In order to validate this IBM model, three different benchmark cases are considered. The first two cases focus only on the mould region and the third case considers multi-zonal cooling ambiences. Attempts have been made to compare shell thickness and surface temperature predicted by this model with reported data. Shell thickness is determined as the thickness of the metal pool which has temperature less than or equal to solidus temperature.
Mould region
Case 1: Geometric properties and thermal boundary conditions for the first validation case are specified in set-1 (Tables 1–3). The effective thermal conductivity for the liquid and mushy phases at electromagnetically stirred mould region has been considered as keffective = mk with m = 2. The density is assumed to be of constant value of 7400 kg/m3 for all phases. Our predictions for solidified shell thickness are compared with radioactive tracer based measurements of Lait et al. 39 (Figure 5). The comparisons show good agreement only with little discrepancy at the near meniscus height. This is to be mentioned that a higher equivalent k-based model of Choudhary and Mazumdar 42 has shown disagreement with the experimental data for a similar case. They reported good agreement only when a conjugate heat transfer simulation is performed. Our computation shows that this IBM based methodology coupled with consideration of the convective heat transfer term and a lower effective k value gives good agreement avoiding the involved computations needed for a conjugate solver.
Data of continuous casting used in the present computations.
Length of different section (m).
WS: Water spray zone, RAD: Radiation zone
Water flow rate (m3/min)

Comparison of predicted shell thickness with experimental result.
Case 2: Geometric properties and thermal boundary conditions for the first case are specified in set-2 (Tables 1–3). Janik and Dyja 43 neglected the convective term due to vertical casting speed and considered a higher value of effective thermal conductivity (keffective = 7k) through all phases to encompass the stirring effects. In this study, the similar approach is tested and in addition to that the same problem is solved using the other approach which considers the convective term due to z-movement of molten metal along with the effective k (keffective = 2k, a lower magnitude) at the liquid and mushy phases. The densities considered for liquid, mushy and solid phases are 7000, 7200 and 7400 kg/m3, respectively.
Chipman and Fondersmith 49 showed that shell thickness is proportional to square root of solidification time for ingot casting. Savage and Pritchard 50 showed that for continuous casting, the solidification time may be estimated as, time = distance from meniscus/casting velocity. They also derived their own correlation for continuous casting. Our predictions are also compared with Chipman–Fondersmith and Savage–Pritchard correlations (Figure 6).

Comparison of predicted shell thickness with available results and Savage–Pritchard correlation.
Since the partial k approach (lower magnitude of equivalent k along with convective term due to movement of billet) results are best comparable with experimental as well as Chipman–Fondersmith and Savage–Pritchard correlation, this method has been used for all remaining computations presented in this article. Our prediction considering the z-directional convection agrees closely with both experimental results and established correlations, and hence, this approach can provide a very good model for continuous casting simulations.
To confirm grid-independence of the solutions, cases are run with two different grid sizes, 80 × 80 × 120 and 120 × 120 × 180. For both of the cases, shell thickness is predicted to be the same.
Multi-zonal cooling ambiences
For validation of multi-zonal ambience cooling, the case similar to the Qing et al. 15 case has been simulated. The partial k approach along with convection effect is used. The set-3 (Tables 1–3) is considered for this case. The densities considered for liquid, mushy and solid zones are 7000, 7200 and 7400 kg/m3, respectively. The effective thermal conductivity, keffective = mk, has been considered for liquid and mushy phases, where m equals 1–2, in different cooling ambiences (or zones). The value of m is equal to 2 for mould and segment-1 and 1 for other segments. The heat extraction in mould for this case has been considered as per Savage–Pritchard 50 heat flux relation. The shell thickness with respect to metallurgical length is plotted for casting speed (u) of 0.085 m/min (Figure 7). Our prediction of shell thickness with partial k approach agrees well with Qing et al.’s result (Figure 7).

Comparison of predicted shell thickness with Qing’s results under multi-cooling ambiences.
Effect of heat extraction rate in mould region
In order to quantify the effect of the rate of heat extraction on shell growth, simulations are performed with different magnitudes of heat extraction rate. Casting velocity is kept constant at 0.0267 m/s for all the cases. The parameters are shown in set-4 to -9 (Tables 1–3). The effective thermal conductivity, keffective = mk, with m = 2 has been considered. The respective heat extraction rates are calculated using

Predicted shell thickness for different heat extraction rates at casting speed of 0.0267 m/s.
As the shell thickness suggests a parabolic behaviour with distance from the mould (or distance travelled due to the casting speed), the shell thickness varies linearly with square root of time (
where A is considered as a constant depending on the material property and heat extraction rate; C is considered as delay constant, which is related to amount of superheat; t is time in minute and d is shell thickness in millimetre. The obtained values of A and C in Figure 9 are within the ranges reported by other authors for low carbon steel. 51 It is found that at lower heat extraction rates, the Chipman and Fondersmith relation does not hold. The relation changes to

Predicted shell thickness for different heat extraction rates on square root of time scale at casting speed of 0.0267 m/s.
Due to this, the shell-thickness profile with respect to
The effect of heat extraction rate slows down and finally ceases at higher values. In this article, heat extraction rate above 19.008 MW/m2 has no further effect on shell thickness, and hence, solidified thickness profile becomes constant for higher heat extraction rates and so the value of A and C in aforesaid correlation.
Effect of casting velocity in mould region
In order to quantify the effect of the casting velocity on the growth of shell thickness inside the mould, three different velocities (0.01335, 0.0267 and 0.0534 m/s) are studied as presented in Figure 10(a) and (b), keeping heat extraction equal to 2.376 MW/m2 in all cases. The other parameters are considered as shown in set-5 (Tables 1–3).

Predicted shell thickness with different casting velocities with heat extraction rate of 2.3764 MW/m2 (a) with distance from meniscus and (b) with square root of time.
The results in terms of solidified shell thickness are presented with respect to distance from meniscus (Figure 10(a)) and square root of time (Figure 10(b)). Increased casting speed leads to high value of time-delay constant, and as a consequence, a low shell thickness is achieved in the mould region (this may lead to bursting). This could be observed in Figure 10(a), which shows that solidified shell starts to grow at higher distances from meniscus as the velocity increases. On the vice-versa, time delay becomes small at lower velocities. Figure 10(b) shows that the effect of velocity on temporal growth of solidified shell inside the mould region is negligible, and hence, Chipman and Fondersmith correlation holds for different casting velocities. However, the effect of velocity could be seen in terms of shell thickness with respect to distance from meniscus inside the mould. The higher speed needs more metallurgical length to achieve same amount of solidified shell thickness as compared to lower casting speed (Figure 10(a)).
Correlation for delay constant
The relation between solidified shell thickness and solidification time can be expressed as sets of modified Chipman–Fondersmith relations as shown in the legend of Figure 9:
Or,
Constants A and C are functions of heat extraction rate and material properties. C denotes the time delay or minimum time needed to start solidification. This is also related with the minimum distance from the meniscus where solidified shell can start forming. For a specific material, information on C at one particular heat extraction rate (considered as reference case) can be utilized to predict values of C for other different heat extraction rates. Correlation for time-delay constant is obtained as
where C is the delay constant, Cref is the delay constant for the reference case and fqm is the heat extraction rate as fraction of the reference case.
For 15G2A grade steel (low carbon steel), the value of A ranges from 26 to 36 and becomes constant for higher heat extraction rates (qm ≥ 19 MW/m2). The correlation-fit data are shown along with the computational results in Figure 11.

Shell thickness correlation and computational results.
Simulation of the full-length billet considering multi-zonal cooling ambiences
Numerical simulation on the full-length billet geometry under different cooling ambiences is performed using IB method in multiple processors (11 processors). The partial k approach (lower values of equivalent k along with convective term due to casting speed) is used. The multi-zonal cooling ambience is chosen as shown in Figure 4(c). The different values of equivalent k are used in different segments, k values are decreased away from the meniscus as the stirring/turbulence effects decrease. The effective thermal conductivity, keffective = mk, has been considered for all temperature above solidus temperature (i.e. liquid and mushy phases), where m varies in the range of 1–3 at different cooling ambiences (or zones). The values of m are considered as m = 3 in the mould, m = 2 at segment-1 and m = 1 for other segments. The caster and casting parameters are chosen according to Tables 1–3 (set-10 to -11).
In Figure 12(a), the solidus line front and 70% solidification front are plotted against the distance from meniscus at two time different levels. The casting velocity is taken as 0.0111 m/s. It shows that once the full domain is solidified, the solidus front or 70% solidus front profile remains the same with respect of time, and the casting problem turns to a case of steady-state heat transfer.

Predicted shell thickness under multi-cooling ambiences at casting speed of 0.0111 m/s (a) with distance from meniscus (at different time levels) and (b) with square root of time.
In Figure 12(b), the shell thickness with respect to square root of time is plotted. It is observed that the shell thickness is proportional to
Figure 12(b) also shows that the 70% solidus line holds
Effect of casting speed on billet or multi-zonal cooling ambiences
The effects of casting velocity in multi-cooling ambiences are presented in Figure 13(a) and (b). The caster parameters are taken according to set-10 of Tables 1–3. The three different casting speeds considered are 0.0055, 0.0111 and 0.0267 m/s, respectively. In Figure 13(a), shell thickness is plotted against the distance from meniscus, and it shows that for higher speeds, longer distance from meniscus is needed for complete solidification, which is obvious as time for solidification decreases with increase in speed. As the casting speed increases to 0.0267 m/s, the metallurgical length of the entire domain (30.2 m) does not emerge to be sufficient to achieve full solidification. However, to show the temporal effect of velocity in shell growth, shell thickness is plotted against

Predicted shell thickness with different casting velocities under multi-cooling ambiences (a) with distance from meniscus and (b) with square root of time.
It is also noteworthy that in Figure 13(a), a very low casting speed shows solidification of the entire casting cross section at a distance of 5 m/min from the mould. This distance increases with increase in speed and no complete solidification is observed for a higher speed of 0.0267 m/min. This observation can be attributed to the fact that with increase of speed solidification, time decreases that supersedes the enhanced convection effects and the molten alloy does not get sufficient time for solidification of the full cross section through its entire length.
So from Figure 13(b), it can be observed that the
Effect of heat extraction rate on billet for multi-zonal cooling ambiences
The effects of heat extraction rate in multi-cooling ambiences are presented in Figure 14(a) and (b). The caster parameters are taken according to set-10 to -11 of Tables 1–3. The casting speed is kept equal to 0.0111 m/s for both the cases. In Figure 14(a), the growth of shell thickness with respect to metallurgical length is presented while Figure 14(b) depicts shell-thickness variation with respect to square root of time. The lower and higher heat extraction rates are provided according to Table 1 and 3 (Set-10 to 11). The case with higher heat extraction rate has double heat extraction rates in every respective section as compared to the lower heat extraction rate. The effect of heat extraction is significant in terms of both metallurgical length and time scale. Increase in the heat extraction rate favours the process of shell solidification. Figure 14(b) also depicts that the slope of solidified shell remains same in both of the cases; however, a constant offset is observed. Hence, it can be inferred that a change in heat extraction rates in multi-zonal cooling ambience introduces a time delay in the solidification process.

Predicted shell thickness for different with different heat extraction rates at casting speed of 0.0111 m/s (a) with distance from meniscus and (b) with square root of time.
The variation of shell thickness at different cooling ambiences can also be explained via temperature contours over an axial plane (Figure 15). The temperature contours for set-10 with casting speed of 0.0111 and 0.0267 m/s and set-11 with casting speed of 0.0111 m/s are shown for metallurgical length of 16 m. It can be seen that high temperature zone expands up to the bottom plane in set-10 with increased casting speed. The solidified region (with temperature less than solidus temperature, 1722 K) increases especially near the lower planes for set-11 with casting speed of 0.0111 m/s, which has double of heat extraction rate than set-10 with casting speed of 0.0111 m/s.

Temperature profile for different velocities and heat extraction rates.
Conclusion
This article demonstrates implementation of an efficient IBM based finite difference scheme for prediction of solidified shell thickness and temperature profile in a continuous caster. The proposed methodology shows good agreement with previously reported data as well as linear scalability over multiple processors. It is also shown that an increased effective conductivity coupled with convection term due to casting speed can better model the heat transfer augmentation due to electromagnetic stirrer. Solidification patterns are observed with respect to distance from the meniscus and time of solidification. Time-delay in solidification is found to be a function of heat extraction rate from the mould. A correlation has been proposed to determine the effect of heat extraction rate on growth of solidified shell in the mould region. It is shown that at a lower heat extraction rate,
Different cooling ambiences (as used in an industrial continuous caster) have been tested and their effect in the shell thickness has been studied. It is observed that higher heat extraction rate and lower casting speed favour better solidification. The 3D heat transfer and solidification patterns are visualized to explore thermal physics of continuous casting. Graphically, the metallurgical length could be determined with change in casting speed if the data for growth of shell thickness for one case is known.
The present results and the obtained correlations can be utilized to obtain a pool of information for different cooling parameters and casting speeds, which could be helpful while designing the caster for better productivity and to avoid bursting or breakout. The proposed IB model also offers potential for simulation of complex casting configuration with moving boundaries. Furthermore, this method can be used to study the heat transport phenomena in other manufacturing processes which involve complex geometries with fixed or moving surfaces (like investment casting, friction stud welding).
Footnotes
Appendix 1
Declaration of conflicting interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work has been financially supported by BRNS grant 2012/34/68/BRNS/2972. Computations are performed using HPC infrastructures at IIT Patna.
