A Numerical Solution Model for the Heat Transfer in Octagonal Billets

Abstract

In the quest for high-quality steel products, the need of cast billets with minimum surface and internal defects is of paramount importance. On the other hand, productivity is required to be as high as possible in order to reduce production cost. Different billet shapes have been applied with emphasis upon square, rectangular, and circular cross-sections. It is obvious that the best billet shape that minimizes surface and subsurface defects is the circular one. Nevertheless, this shape creates some problems with respect to handling and safety reasons. One recent attempt is to produce normal octagonal-shaped billets that appear to approach the circular shape albeit easier to handle. In this study, a numerical solution for the heat transfer during solidification in the continuous casting of octagonal billets has been carried out. The developed model deploys an implicit scheme in order to solve the differential equations of heat transfer under the appropriate boundary conditions in a section of an octagonal billet, assuming fully axisymmetric cooling of the bloom. The geometry of the octagonal billet plays an interesting role in the development of the heat transfer analysis. Based upon fundamental principles, a computer program has been developed for this purpose. Consequently, results from the numerical solution are presented and discussed.

1. Introduction

The mathematical analysis of heat transfer in a billet during the continuous casting process has been investigated from the early stages of the development of continuous casting machines. In this study, some fundamental works related to the numerical solution of heat transfer in solidifying billets are presented. A pioneer work is mentioned in [1] in which a mathematical model of heat transfer in continuously cast steel slabs is described. A computer program was developed to model the problem using a 1D transient conduction equation with the appropriate boundary conditions. Some years later, the design of the mold and spray sections of a continuous casting machine were examined in detail with the aid of a 2D heat-flow mathematical model [2]. Experimental data obtained from commercial casters were used for the validation of the deduced results. In that monumental work, the various aspects of heat transfer in the various sections of a casting machine were analyzed in detail. Later on, surface and subsurface defects in the cast billets became apparent, and research work included mechanical phenomena (like stress and strain) together with the heat transfer analysis; an early typical work on the subject was presented in [3]. It is worth mentioning that the caster automation includes most of the heat transfer formulation and analysis, albeit not published. Consequently, the present literature survey mostly presents indicative published works upon heat transfer analysis. In the 1990s, the continuing casting of steel started maturing, and this was reflected by two studies with practical implementations [4, 5]. In a review paper [6], the evolution of heat transfer and mechanical studies for the continuous casting were presented. Furthermore, an analysis on the ideal taper prediction for billet casting and a thorough analysis of thermomechanical behavior in billet casting were also presented [7, 8]. The salient features of heat transfer in the secondary cooling zone during the continuous casting of steel were examined in [9]. Primary cooling is considered the start-up of steel solidification process in a water-cooled copper mold, and secondary cooling zone is considered the region just after the mold until the end of the caster in which the billet is cooled by spray nozzles or air-mist systems. Radiation heat exchange between the billet surface and the environment plays also an important role in the secondary cooling zone. Fluid-flow phenomena in the mold captured the interest, and as of that, heat transfer studies were also reported on the subject later on; a typical work is given in [10]. In most similar works, general-purpose heat-transfer software has been deployed either as a direct prediction tool or as a verification one. A typical work is given in [11] in which a 1D heat-transfer simulation program (CON1D) was successfully verified against a 3D finite element model developed in Abaqus. Near-net-shape cast sections appeared in industrial practice in the early 2000, and therefore, their thermal-mechanical modeling during casting was developed as well [12]; in this type of analysis, the use of a package like Abaqus has been indispensable. In a relatively recent work [13], the 1D differential equation for heat transfer was solved in order to study the evolved microstructures during the solidification of round billets (so, the radial direction for heat transfer prevails). In the cases under study, the billet diameters were in the range from 210 up to 410 mm; the validation of their model was performed by deploying a semi-analytical model [14] for the prediction of the surface temperature of a cast round billet. The numerical model was used to calculate the local solidification time, the local gradient, and the local solidification time as a function of the distance from the round billet surface; furthermore, a simplified relation for the prediction of the columnar to equaxial transition was proposed. Deploying the finite volume method, a real-time 2D heat transfer and solidification model for continuous casting [15] were developed and tested online. The behavior was satisfactory with less than 10°C deviations on surface central temperatures. A dynamic heat transfer model was developed [16] in order to study the effect of the casting speed on the temperature field of continuously cast billets; various steel grades were taken under consideration. A real-time heat-transfer model and a heat-transfer coefficient identification method [17] were developed. The model validation was achieved by the measurement of surface temperatures by a CCD system that appeared to eliminate the impact of scales on the temperature measurement and keep the measured surface temperature fluctuation within the range of ±10°C. A finite volume method and the alternating direction implicit algorithm (ADI) were selected in order to develop a real-time heat transfer model [18] for the dynamic control of continuous cast billets. The system was applied online to control the electric current of the final electromagnetic stirring (FEMS) system. A simulation model of solidification and heat transfer of 150-mm-square billets was developed [19] for the continuous casting of grade 40 Cr. Nailing tests and temperature measurements were applied in order to fine-tune the model with maximum errors of less than 2%. In another study [20], different micro-segregation models coupled with fluid flow and heat transfer were run to study macrosegregation phenomena in a round billet with a 165-mm diameter. It was predicted that heavy centerline segregation occurs at the end of solidification when the central solid fraction reaches the value of 0.7. The need for high-speed casting under very controllable conditions has led the researchers to seek billet sections that can succeed in the demanding steelmaking environment. Recently, one of the leading manufacturers in the steelmaking sector, Danieli, has announced [21] the octagonal billet as a proven potential shape for successful high-speed casting. Consequently, the present author has taken the opportunity to study the heat transfer in a solidifying billet based on fundamental principles. The work has been carried out for a larger billet size than the one tested in [21] and for two special steel grades, the peritectic grade S355 J2 and the medium-high-grade 42CrMo4, that are both important for many end-user applications.

2. Formulation of the heat-transfer problem

The general three-dimensional (3D) heat-transfer equation that describes the temperature distribution inside a solidifying object is given by Eq. (1) according to [22]

$$\mathit{\rho c}\frac{\partial T}{\partial t}=\nabla \cdot k\nabla T+S$$ E1

The source term S may be considered [23] to be of the form:

$$S=S_C+S_{P}T$$ E2

Furthermore, T is the temperature, and ρ, c, and k are the density, heat capacity, and thermal conductivity, respectively. The 2D heat-transfer equation in Cartesian coordinates can be written as [22]

$$\mathit{\rho c}\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left(k\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(k\frac{\partial T}{\partial y}\right)+S$$ E3

The boundary conditions that apply in the octagonal billet case are very similar to the ones for the rectangular billet that have been presented in detail in [24] and will not be repeated here. Nevertheless, two important boundary conditions apply in the case under study, which follow:

• Due to symmetry, the upper diagonal side (Figure 1, segment OB) of the domain of interest is supposed to be adiabatic, as well as the lower side (Figure 1, segment OK) is. In this way, the following formulation holds:

$$\begin{array}{l}\frac{\partial T}{\partial s}=0,\text{where}s\text{is the vertical axis}\mathrm{on}\text{the line}\\ y=x\text{tan}\left(\pi /8\right)\end{array}$$ ₍₄₎ E4

• The primary (mold) and secondary cooling systems are the ones applied in Stomana; this is presented in detail in a previous publication [24] and will not be repeated here. However, due to the potential of the octagonal mold, a surplus of water was used in favor of an enhanced cooling behavior for the solidifying billet; this was expressed as an extra percentage of water flow deployed for the octagonal billet.

Figure 1.

The salient features of an octagon and its circumscribed circle.

3. Theoretical solution

There is no published work on the theoretical solution for the octagonal billet until now. However, a theoretical solution can be approximated by solving the DE of heat transfer on an equivalent circular geometry. According to Figure 1, the selected radius of that circle was R_e = 156.6 mm, as given by Eq. (5):

$$R_e=\frac{1+\text{cos}\left(\pi /8\right)}{2}R=0.962162.84=156.6\mathrm{mm}$$ E5

This radius is actually the mean value between the circumradius and in-radius of the octagon [25]. The circumradius R being equal to 162.84 mm makes the currently studied octagonal billet equivalent in size to a 250 × 300 mm² billet produced in the Stomana plant, at Pernik, Bulgaria. Now the original problem can be converted into a heat-transfer problem of cylindrical geometry, and the DE together with the boundary conditions is formulated according to Eq. (6):

$$\begin{array}{l}\frac{\partial T}{\partial t}=\alpha \left(\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}\right)\\ T=T_0,t=0;\frac{\partial T}{\partial r}=0,r=0;-k\frac{\partial T}{\partial r}=h\left(T-T_f\right),r=R_e\end{array}$$ E6

By variable transformation using Eq. (7)

$$\mathrm{\Theta }=\frac{T-T_f}{T_0-T_f}$$ E7

the boundary condition at r = R_e becomes:

$$\frac{\mathrm{\partial \Theta }}{\partial r}+\frac{h}{k}\mathrm{\Theta }=0,r=R_e$$ E8

In order to fulfill the boundary condition at r = R_e according to Eq. (8), Eq. (9) has to be solved for β:

$$\beta J_1\left(\beta \right)=A{J}_0\left(\beta \right),A=\frac{R_{e}h}{k}$$ E9

where J₀ and J₁ are the Bessel functions of zero and first order, respectively.

This theoretical solution is presented in the book of Carslaw and Jaeger [26] and is followed here. There are an infinite number of values for β, which solve Eq. (9); nevertheless, the first six roots [26] are enough for the computations. In this way, the derived solution arrives in the form of Eq. (10):

$$\begin{array}{l}\mathrm{\Theta }=\frac{T-T_f}{T_0-T_f}=\sum _{n=1}^{\infty }\text{exp}\left(-{{\beta }_n}^2\xi \right)\frac{2A{J}_0\left({\beta }_{n}r/R_e\right)}{\left({{\beta }_n}^2+A^2\right)J_0\left({\beta }_n\right)}\\ \xi =\frac{\alpha t}{{R_e}^2},A=\frac{R_{e}h}{k}\end{array}$$ E10

Eq. (10) was applied in order to compute the temperature distribution T(r,t) inside the equivalent to the octagonal geometry cylindrical billet; some results are depicted in Figure 2. In fact, the following parameters were used in these computations: thermal conductivity, k = 50 W/m/K; density, ρ = 7600 kg/m³; heat capacity, c = 400 J/kg/K; thermal diffusivity, α = k/(ρc) = 1.645 × 10⁻⁵ m²/s; heat transfer coefficient, h = 200 W/m/K; cooling water temperature, T_f = 30°C, and initial billet temperature, T₀ = 500°C.

Figure 2.

Theoretical predictions of the temperature distribution in various positions in a cylindrical billet.

4. Numerical solution

Figure 1 illustrates the important geometrical features of the octagonal-shaped billet. Considering that the flow of heat takes place only in the radial direction toward the center of the octagon, and certainly toward the center of its circumscribed circle, it appears that due to symmetry, only a small trigonal part becomes the important region (the shaded area in Figure 1) for the solution of the DE of heat transfer. Actually, when uniform cooling takes place around the solidifying billet, then symmetrical conditions prevail in the relevant heat transfer. Consequently, heat flow takes place only in the radial direction. The strongly implicit scheme as presented by Patankar [23] was deployed for the solution of the DE of heat transfer. Figure 3 illustrates the control-volume selection for a small number of nodal points.

Figure 3.

Presentation of control volumes and nodal points for a coarse grid.

The five-point formulation for the discretization equations was applied. Special attention was paid in the boundary condition formulation of the discretization equations, and a more detailed analysis on this is presented in Appendix I. One may realize that once the boundary conditions are properly written, then a finer grid will ultimately permit the convergence to the engineering solution. In Figure 3, only the lower trigonal part is necessary to be examined. The vertical side, MZ, is the boundary from which the main heat flow (cooling) takes place; the horizontal side OZ and the diagonal side OM, with an inclination of π/8, or 22.5°, are considered to be adiabatic to heat flow due to the aforementioned symmetry reasons. This study is part of a series of published works with respect to the numerical solution of the heat transfer equation in 2D (square and rectangular) and 3D domains [24, 27, 28, 29]. Although the core of the existing program remained intact, due to the nature of the present trigonal domain, the part of the program related to the boundary conditions had to be developed from scratch. The octagonal billet with R = 162.84 mm can be discretized with a fine grid of 200 × 200 nodal points. In this way, the space intervals (δx, δy) were about δx = R cos(π/8)/199 = 0.756 mm and δy = R sin(π/8)/199 = 0.313 mm, respectively. Actually, keeping the ratio, δy/δx = 0.313/0.756 = 0.414 = tan(π/8), exhibited a relatively good convergence; the time interval was Δt = 1 sec.

5. Running the computer program

The Gauss-Seidel algorithm was applied for the iterative solution of the matrix equations in the case of the 42CrMo4 grade. Over-relaxation was used for the fastest possible convergence, and the over-relaxation parameter used was ω = 1.870, which has exhibited good behavior for this type of studies [30]. The applied maximum error (tolerance) for the computed temperature at each nodal point was 5 × 10⁻³. In addition to this, the Strongly Implicit Procedure (SIP) [31] was deployed for the iterative solution of the matrix equations in the case of the S355 J2 grade. In this case, however, the same tolerance (5 × 10⁻³) was used for the maximum residual value and not the temperature difference at each nodal point. The C++ language was selected for the development of the program, and the deployed compiler was the Intel Parallel Studio XE 2015 Composer Edition for C++; the open source, cross-platform Code::Blocks, version 17.12 software was used as an integrated development environment for program edition and compiler invoking. The software was run in a DELL Alienware laptop with the Intel i7-6700HQ CPU (8 cores) at 2.6 GHz, 16GB RAM, running under a 64-bit Windows 10 Professional O/S.

6. Results and discussion

In order to validate the new developed program, a solution against the aforementioned theoretical one had to be sought. For this reason, an octagonal billet initially at 500°C, cooled afterwards by an hypothetical fluid at 30°C and with a heat transfer coefficient h = 200 W/m²/K, was simulated. The values for the properties of density, thermal conductivity, and heat capacity were the same with the ones used in the theoretical solution. In total, 38 sets of data were randomly selected at various positions inside the cylindrical billet and were compared with similar results deduced from the numerical solution described here in order to verify its validity. Figure 4 presents the compared results of these two sets of computed data. An analysis of variance (ANOVA) was performed for these two sets of results, and the following statistical results were derived: residual standard error, 3.133 on 36 degrees of freedom (DF); multiple-R squared, 0.9957; and F-statistic, 8333 on 1 and 35 DF, with a p-value <2.2 × 10⁻¹⁶. Consequently, a standard error of 3.1°C was found to represent the mean statistical error between the results obtained from the theoretical and current numerical solutions, respectively.

Figure 4.

Comparison of temperature values between the theoretical (for a cylindrical billet) and numerical (for an octagonal billet) solutions.

The advantage of studying the heat transfer in an octagonal billet marks the importance of such an experiment: the potential to test some results against theoretical ones gives the confidence about attaining the proper numerical solution; in fact, the octagonal geometry approximates the circular cross-section better than a square one for this comparison to be accomplished.

Table 1 presents the chemical analysis for the two selected grades under examination: 42CrMo4 and S355 J2. In the bloom caster of Stomana, in Pernik, the billet size of 250 × 300 mm² is normally used for the production of special steels. The equivalent billet to this size in octagonal shape has a circumradius of R = 162.84 mm, as aforementioned. An interesting idea is that putting the octagonal billet in practice, there is a potential to increase productivity by using more cooling water, an advantage that has been tested in practice [21] for normal rebar grades. Consequently, in this study a surplus of mold water by 12% for both of the two grades 42CrMo4 and S355 J2 was used compared to the normal water-cooling process applied in Stomana for these grades, respectively. Furthermore, a surplus of water in the air-mist region by 3% for the high-carbon grade 42CrMo4 and 4% for the peritectic grade S355 J2 was applied, respectively.

Grade	%C	%Si	%Mn	%P, S max	%Cr	%Mo	%Ni	%Cu	%V	%Al
1	0.42	0.22	0.78	0.015	0.98	0.20	0.07	0.22	—	0.03
2	0.17	0.21	1.06	0.015	0.10	0.02	0.09	0.22	0.03	0.03

Table 1.

Chemical analysis for grades 42CrMo4 (1) and S355 J2 (2).

Figure 5 illustrates the shell thickness (curve 3) increase as an octagonal billet of grade 42CrMo4 travels down Stomana’s continuous casting machine. One may notice that the solidification completes at about 32.4 m from the liquid-steel meniscus in the mold.

Figure 5.

Surface (1) and centerline (2) temperature values in an octagonal billet of the 42CrMo4 grade; thickness evolution (3) as the billet solidifies downwards the continuous caster.

At the same time, the centerline temperature (curve 2) drops appreciably at that point revealing the complete solidification at that point, as well; the surface temperature is also presented by curve 1. The casting speed was 0.70 m/min, and the superheat (SPH) was 30°C in the computation. The SPH is actually the excess temperature above liquidus temperature; the liquidus temperature is exclusively a function of the chemical composition of a steel grade.

Figure 6 depicts the temperature distribution in the domain of interest of a 42CrMo4-grade octagonal billet at a position of 19.2 m from the meniscus. Carbon steels normally solidify gradually from liquid to solid phase. The degree of solidification is generally described by a parameter that is called solid fraction (f_s) and represents the percentage of the solid phase in the mixture of solid and liquid phases. When f_s = 0, then we have 100% liquid phase, and the steel temperature at this point is the liquidus temperature, T_l. When f_s = 1, then we have 100% solid phase, and the steel temperature at that point is the solidus temperature, T_s. Although the liquidus temperature is always a function of the chemical composition of a steel grade, this is not the case with the solidus temperature. The solidus temperature is also affected by the local cooling rate at solidification (C_R), which is expressed as

$$C_R=\frac{T_P-{T_P}^0}{\mathrm{\Delta }t}$$ E11

Figure 6.

Temperature distribution in the selected part of a 42CrMo4-grade octagonal billet. Solid fraction values (f_s) are presented in the mushy zone.

Eq. (11) shows that the temperature difference between a temperature T_P and the previous one T_P⁰ at point P inside a solidifying billet within a time interval Δt, divided by this time interval, defines that local cooling rate. Consequently, the solidus temperature may be given by a formula of the type:

$$T_S=G_1\left(\mathit{\text{chemical analysis}},C_R\right)$$ E12

On the other hand, the solid fraction may be considered a function of the local cooling rate and temperature:

$$f_S=G_2\left(T,C_R\right)$$ E13

Therefore, during solidification of carbon steels, there is always a range between solidus and liquidus temperatures in which the solid fraction is in the range from 0 to 1; this physical range inside a solidifying steel body is called mushy zone, and the whole phenomenon is related to micro-segregation. The simple micro-segregation model [32] has been adopted in this work, and the analysis has been described in similar previous studies [24, 29, 33]. It appears that the larger the carbon content in a steel grade, the larger the mushy zone gets, and central porosity becomes inevitable in the final stages of solidification.

The local cooling rate distribution for the case of a 42CrMo4-grade at the caster position presented in Figure 6 is illustrated in Figure 7.

Figure 7.

Cooling rate distribution in the selected part of a 42CrMo4-grade octagonal billet.

A short analysis showed that for the data presented in Figure 7, the C_R (cooling rate in °C/s) has the following statistics: average value μ = 0.106, standard deviation σ = 0.128, minimum = 0.023, and maximum = 0.494. More than 99% of the C_R population is within μ + 3*σ = 0.490.

Figure 8 depicts the surface (1) and centerline (2) temperatures in an octagonal billet as the billet solidifies downstream the Stomana billet caster; the shell thickness (3) is presented as a function of the distance from the meniscus of liquid steel in the mold. It is worth mentioning that a S355 J2-billet cast at a speed of 0.85 m/min solidifies even faster than a similar 42CrMo4-grade billet at a lower speed. Figure 9 illustrates the temperature distribution in the critical geometrical region of a S355 J2-grade octagonal billet at a distance of 18.4 m from the meniscus, cast at a speed of 0.85 m/min and a SPH = 30°C.

Figure 8.

Surface (1) and centerline (2) temperature values in an octagonal billet of the S355 J2 grade; thickness evolution (3) as the billet solidifies downwards the continuous caster.

Figure 9.

Temperature distribution in the selected part of a S355 J2-grade octagonal billet. In addition to the solid fraction values (f_s) that are presented in the mushy zone, the percentages of δ- and γ- phases are shown, respectively.

Peritectic grades have the tendency to crystallize in a mixture of delta (δ) and gamma (γ) phases of iron in an intermediate temperature range, before continuing to a 100% γ-phase solidification. For this reason, the percentages of δ- and γ-phases are also presented in Figure 9 at the selected solid fractions.

All these calculations come also from the simple micro-segregation model [32] that is adopted in the developed program. Figure 10 depicts the local cooling rates at the same position from the meniscus as for the aforementioned temperature distribution presented in Figure 9.

Figure 10.

Local cooling rate distribution in the selected part of a S355 J2-grade octagonal billet.

For the cooling rate distribution (C_R in °C/s) presented in Figure 10, a short statistical analysis derived the following results: average value μ = 0.138, standard deviation σ = 0.156, minimum = 0.027, and maximum = 0.490. Almost 99% of the population is within μ + 3*σ = 0.606.

The minimum average number of iterations required for convergence by the Strongly Implicit Procedure (SIP) was attained at the value of 0.95 for the iteration parameter α as explained by Stone [31] in order to speed up the convergence process. Actually, this value (α = 0.95) was used in order to solve the derived system of matrix equations in the case of the S355 J2 grade. Figure 11 illustrates the effect of this iteration parameter on the ratio of the required average number of iterations to the minimum number of iterations required at α equal to 0.95. In general, the convergence behavior improved as the selected iteration parameter increased. It is worth noting that convergence succeeded within an extensive range of the iteration parameter although above the limiting value of α = 0.95, the solution procedure started to diverge. For verification purposes, two sets of temperature results were generated by the Gauss-Seidel and Strongly Implicit Procedure, respectively, taken after a long simulation time (at the same time instant t = 2460 sec, equivalent to a distance of 34.85 m from meniscus). Using R [34], a Pearson’s product–moment correlation test gave a coefficient of 0.9999344 for these two sets of results with a t-Student value equal to 17,455, 39,998 degrees of freedom and a p-value < 2.2 × 10⁻¹⁶.

Figure 11.

Effect of the iteration parameter on the convergence behavior of the strongly implicit procedure.

7. Conclusions

The 2D differential equation of heat transfer was numerically solved for the derivation of the temperature distribution inside an octagonal billet. The deduced model was applied for the casting case of octagonal billets for an equivalent shape to the size of 250 × 300 mm² billets that are normally cast at the Stomana billet caster. Furthermore, the simulated octagonal billets were considered to be from the two special grades of 42CrMo4 and S355 J2. Two different iterative methods were applied for the solution of the resulting system of equations, the Gauss-Seidel method and the strongly implicit procedure, both exhibiting satisfactory behavior. A surplus of cooling water especially at the mold may result in a potential productivity increase. Computational results from the simple micro-segregation model were also included in the present study.

Acknowledgments

The author would like to express his gratitude to the top management for the continuous support on this kind of studies.

A.

The applied solution is based on the implicit scheme as described by Patankar in 1980 [23]. Figure 3 shows the control volumes and the nodal points upon which the numerical solution was based; in the graph, only a very rough discretization was considered with 6 × 6 nodal points just for illustrative purposes. Normally, for nodal points in the interior of the geometry under investigation, the implicit scheme works on five adjacent points; in this way, for every internal point P, the four adjacent nodal points of interest are nominated as E (east), W (west), N (north), and S (south). The numerical discretization equation for the heat transfer equation at P is then given by Eq. (14)

$$a_P{T}_P=a_E{T}_E+a_W{T}_W+a_N{T}_N+a_S{T}_S+b_P$$ E14

where the temperature coefficients and the constant term are given by Eq. (15):

$$\begin{array}{l}{a}_E=\frac{k_e\mathrm{\Delta }y}{{\left(\mathit{\delta x}\right)}_e},a_W=\frac{k_w\mathrm{\Delta }y}{{\left(\mathit{\delta x}\right)}_w},a_N=\frac{k_n\mathrm{\Delta }x}{{\left(\mathit{\delta y}\right)}_n},a_S=\frac{k_s\mathrm{\Delta }x}{{\left(\mathit{\delta y}\right)}_s},\\ {a_P}^0=\frac{\mathit{\rho c}\mathrm{\Delta }x\mathrm{\Delta }y}{\mathrm{\Delta }t},b_P=S_C\mathrm{\Delta }x\mathrm{\Delta }y+{a_P}^0{T_P}^0,\\ a_P=a_E+a_W+a_N+a_S+{a_P}^0-S_P\mathrm{\Delta }x\mathrm{\Delta }y\end{array}$$ E15

The described formulation is proven to satisfy the energy balance within the control volume that is included by the sides “n,” “e,” “w,” and “s.” Extra care is required in the derivation of the discretization equations at the borders of the studied geometry. For example, Eq. (16) describes the formulation for the temperature of a typical point on the outer area of the octagonal billet that is cooled by the water-cooled copper mold; the heat transfer is driven by a lower water temperature (T_f), by a rate given by a heat transfer coefficient h:

$$\begin{array}{l}{a}_C{T}_C=a_G{T}_G+a_D{T}_D+a_F{T}_F+b_C\\ a_G=\frac{k_w\mathrm{\Delta }y}{{\left(\mathit{\delta x}\right)}_w},a_D=\frac{k_n\left(\mathrm{\Delta }x/2\right)}{{\left(\mathit{\delta y}\right)}_n},a_F=\frac{k_s\left(\mathrm{\Delta }x/2\right)}{{\left(\mathit{\delta y}\right)}_s}\\ {a_C}^0=\frac{\mathit{\rho c}\mathrm{\Delta }\mathit{vol}}{\mathrm{\Delta }t},b_C=S_C\mathrm{\Delta }\mathit{vol}+{a_C}^0{T_C}^0+h{T}_f\mathrm{\Delta }y\\ a_C=a_G+a_D+a_F+{a_C}^0-S_P\mathrm{\Delta }\mathit{vol}+h\mathrm{\Delta }y\\ \mathrm{\Delta }\mathit{vol}=\left(\mathrm{\Delta }x\mathrm{\Delta }y/2\right)\end{array}$$ E16

Similarly, another very interesting point of analysis on the derivation of the discretization equation is at a diagonal point A; the formulation is presented by Eq. (17); in this case, there is no heat transfer (due to symmetry) in the area above the diagonal (adiabatic formulation):

$$\begin{array}{l}{a}_A{T}_A=a_B{T}_B+a_N{T}_N+b_A\\ a_B=\frac{k_e\left(\mathrm{\Delta }y/2\right)}{{\left(\mathit{\delta x}\right)}_e},a_N=\frac{k_s\left(\mathrm{\Delta }x/2\right)}{{\left(\mathit{\delta y}\right)}_s},\\ {a_A}^0=\frac{\mathit{\rho c}\mathrm{\Delta }\mathit{vol}}{\mathrm{\Delta }t},b_A=S_C\mathrm{\Delta }\mathit{vol}+{a_A}^0{T_A}^0\\ a_A=a_B+a_N+{a_A}^0-S_P\mathrm{\Delta }\mathit{vol},\mathrm{\Delta }\mathit{vol}=\left(\mathrm{\Delta }x\mathrm{\Delta }y/4\right)\end{array}$$ E17

Nomenclature

c

heat capacity in J/kg/K

f_s

solid fraction

h

heat transfer coefficient in W/m²/K

k

thermal conductivity, in W/m/K

S, S_c

source term, constant term of S, in W/m³

S_p

source term compensating for the temperature dependence of S, in W/m³/K

SPH

superheat (casting temperature—liquidus temperature) in °C

t, Δt

time, time interval in sec

T

temperature in a position inside a billet, or at a nodal point, in °C

T_l, T_s

liquidus, solidus temperature, in °C

x, y, z

spatial coordinates, in m

δx, δy

distance between adjacent nodal points along the x- and y-axis, in m

Δx, Δy

control-volume dimensions along the x- and y-axis, in m

ρ

density, in kg/m³

ω

over-relaxation factor

Subscripts 0

referring to an initial value

f

fluid (referring to a cooling fluid)

Superscripts 0

referring to a value at a previous time step