TranSIenT heaT TranSfer and SolIdIfIcaTIon ModelIng In dIrecT-chIll caSTIng uSIng a generalIzed fInITe dIfferenceS MeThod

The goal of this work is to achieve a novel solution of the transient heat transfer problem in the start-up phase of direct-chill casting processes using a Generalized Finite Differences Method. This formulation is applied in order to solve the heat transfer equation in strong form under a Lagrangian description. The boundary conditions incorporation is done in a simple and natural way. The meshfree nature of this approach allows to capture the motion and phase boundaries evolution without using remeshing approaches. The simplicity, efficiency and suitability of this numerical formulation is demonstrated by comparison of the obtained numerical results with the results already published by other researchers. This shows that our approach is promising for the numerical simulation of heat transfer problems during the start-up phase of direct-chill casting processes.


Introduction
Semi-continuous casting techniques are the most used processes for mass-production of aluminum alloys. One of these techniques is the so-called "Direct-Chill Casting" method (DCC) [1]. It begins with a step in which molten metal is poured into a water-cooled mould where heat from the nearby mould walls is extracted causing its solidification. Once the metal leaves the mould cooling zone, it is further cooled in a secondary cooling zone where water is sprayed on the casting surface. This cooling process continues until the casting reaches the possible maximum length [2]. One of the most important phases in the DCC process is the start-up phase which is the period from the operation start to the time when a steady state is achieved. This phase is highly important since many defects in the end product could originate during this step. In order to get homogeneous casting products, the control of process parameters as cooling rates, casting speed, inlet velocities, and the casting temperature is needed. Notwithstanding, it is hard to optimize and enhance these techniques with experimental studies because it is very tough to quantify temperatures, stresses fields, pressures, or velocities, especially within the mould region during the start-up phase where complex physical phenomena exist [3]. Numerical modelling and simulation is commonly used to improve and prevent the occurrence of casting defects because it produces a lot of detailed information which cannot be obtained with other kind of techniques [4].
The numerical modelling and simulation of physical phenomena involved in continuous casting (CC) and DCC processes are usually based on classical mesh -based approaches [5]. For example, several numerical algorithms based on Finite Element Method (FEM) have been developed for the modelling of three-dimensional fluid flow, heat transfer, and solidification in CC processes [6], the analysis of convective heat transfer in the molten metal and phase change in aluminum DCC [7], prediction of stresses, strains, mushy zone length, and heat transfer in DCC [8], the study of cooling and solidification of semi-continuous casting processes of copper [9], determination of two-fluid flow and the meniscus interface movement in an electromagnetic CC of steel [10], heat transfer study in the primary cooling zone in DCC using experimental measurements of the ingot, mould, and cooling water temperatures during casting [11], the simulation of mould cavity filling process, study of the influence of molten steel flow, and the pouring types on the solidification kinetics at the initial stage in CC [12].
Other numerical approaches based on the Finite differences Method (FDM) have been used to model the steady-state three-dimensional heat flow in CC of steel [13], the numerical prediction of 3D laminar or turbulent liquid flow, heat transfer, and macroscopic solidification in DCC of aluminum alloys to investigate the temperature distributions and solidification patterns in the mould and post mould regions [14], and the analysis of the steady-state heat transfer phenomena in CC [15]. Similarly, numerical procedures based on the Finite Volume Method (FVM) have been proposed for the prediction of turbulent flow and temperature in complex flows of molten steel in CC [16], and the modeling of heat transfer and solidification in CC, considering the establishment of experimental non-uniform and uniform water distributions for the cooling conditions [17].
Meshless or meshfree methods have been developed in recent decades as an alternative to the classical mesh-based techniques since they allow to overcome some of the drawbacks in such methods [18]. The advantages of meshfree methods over meshbased methods rely on the fact that they use a set of finite nodes scattered within a domain of interest as well as on its boundaries, without needing some information with respect to the connectivity and relations between particles such that they do not constitute a mesh of elements. This makes them very attractive in problems involving high deformations and discontinuities in the computational domain without employing remeshing techniques. As a result, it gives the freedom to remove or incorporate particles whenever and wherever required, allowing a simple development of adaptive strategies.
The Local Radial Basis Function Collocation Method (LRBFCM) is a strong form meshless method that has been used for the numerical simulation of CC process of steel considering turbulent fluid flow and solidification [2], the simulation of transient heat transfer of the start-up phase in DCC of aluminum alloys [19], and the simulation of DCC under the influence of a low-frequency electromagnetic fields [20]. Other strong form meshless method approach that has been tested for the modeling of heat transfer and solidification in CC processes in the primary and secondary cooling regions is the Oñate's Finite Point Method (FPM) [21].
Apart from the mentioned meshless methods given in strong form, there exist some other weak form approaches such as the element-free Galerkin method (EFG) that has been used for the simulation of transient heat transfer of the start-up phase in DCC of an aluminum alloy round billet considering the nonlinear aspects related to the material properties and boundary conditions [22]. Recently, another weak form meshless method that has been tested in the analysis of the solidified shell thickness, the mushy zone thickness, the metallurgical length, and residual stress in CC processes is a 3D thermo-elastoplastic Petrov-Galerkin (MLPG) method [23].
A truly meshfree Generalized Finite difference Method (GFDM) is the Finite Pointset Method (FPM) proposed by J. Kuhnert [24]. It has shown to be much better than other meshfree methods and the longestablished mesh-based methods for treating multiphase or free surface flows, fluid mechanics problems involving quickly evolving domains, and for problems including heat transfer or convective flows [25,26,27]. It is a Lagrangian method of strong-form that utilizes a weighted least-squares (WLSM) method to solve elliptic partial differential equations and to compute spatial derivatives [28]. It has numerous advantages when compared with other numerical techniques because it is capable of simply and naturally incorporating any form of boundary condition without requiring some stabilization or special treatments, and its implementation is straightforward [27]. Nevertheless, the transient problem arising in semi-continuous casting processes has not been tested yet. Consequently, the application of FPM to analyze the transient heat transfer and solidification problems in direct-chill casting is proposed in this work; to the author's knowledge this is the first time that this approach has been applied in order to solve this practical engineering industrial process in metal casting. With the purpose of getting some information on its performance, the numerical solution of the start-up phase in a DCC example is compared with numerical published data predicted by other researchers using different numerical methods. The structure of this article is as follows: Section 2 describes the partial differential equation to be considered. Section 3 describes the ideas behind FPM and the numerical scheme applied for the solution of the governing equations. The numerical examples and their results are presented in Section 4. Finally, in the last section some conclusions are given.

governing equations
The modeling of the transient heat transfer and solidification in DCC will be done by means of the heat transfer equation which in Lagrangian form is given by (1) where c is the specific heat, k is the conductivity, T is the temperature, ρ is the density and t is the time. The problem is defined with proper initial and where ∂Ω indicates the kind of boundary, ∂Ω k denotes an inflow boundary condition and ∂Ω l refers to a wall boundary condition, q is the heat flux density, n is the outward unitary normal vector on ∂Ω and T 0 is the initial temperature. Furthermore, q=0 for isolated boundaries or q=h c (T-T ∞ ) in convective boundaries, where T ∞ is a reference temperature and h c is the convective heat transfer coefficient.

numerical Scheme
In this section we present the ideas behind FPM, and the numerical scheme applied for the solution of this transient heat transfer problem in DCC. We start with a temporal discretization of Eq. (1) using an implicit Euler scheme, which leads to (4) where the superscripts indicate the level of time for T and ∆t is the time step. Eq. (4) can be rewritten in general form as (5) where A, B, C and D are defined as: A=ρc, B=∆t∇k, C=k∆t and D=ρcT n [26].

The Finite Pointset Method
In this subsection, the general ideas behind FPM are described, which is a GFDM that uses the WLSM. Following [27,28]: Let Ω be a material domain with boundary ∂Ω and suppose that a set of nodes r 1 , r 2 ,⋯, r N are distributed with corresponding function values f(r 1 ), f(r 2 ),⋯, f(r N ). Then, we are interested in computing the value of f at whatever location f(r) using the function values at node positions in the vicinity of r. A weight function w(r i -r) is proposed in order to determine the number of nodes and the vicinity of r whose form is (6) where r i is the position of the i -th node in the vicinity of r, α is a positive constant whit a value of 6.5 and h is the smoothing length. A Taylor's series approximation of f(r i ) around r is (7) where f k and f kl (f kl =f lk ) depict the set of first and second spatial derivatives at the node position r, r ki and r k are the k -th components of the positions r i and r, respectively, and e i is the truncation error. f k and f kl can be determined minimizing e i for the n p Taylor's series approximations of f(r i ) for to the n p nodes in the r vicinity. The resultant system of equations can be expressed in matrix form as with ∆r ki =r ki -r k , ∆r kli =(r ki -r k )(r li -r l ) and ∆r kki =0.5(r ki -r k ) 2 , where k≠l and k,l = 1,2,3.
The value of a is computed with WLSM minimizing the quadratic form (14) which reads Thereby, the values of the function and its derivatives at r are automatically calculated.

FPM form for general elliptic partial differential equations
General elliptic partial differential equations like Eq. (5) have been studied earlier in [25]. Thus, in this subsection we present the corresponding FPM discretization for these general equations [27].  (19) where f(r j ) represents the unknown function value at the node j and n(j) the number of j-th nodes in the vicinity. As Eq. (19) holds for j = 1,2,⋯ ,N, it forms a full sparse system of linear equations LT = P that can be solved with iterative procedures. Therefore, any kind of solidification and thermal phenomena governed by Eq. (1,5) can be computed with this procedure, just aggregating proper entries in the systems of equations [25,26].

numerical examples and results
With the goal of validating the suitability of this FPM approach to model the transient heat transfer problem in the start-up phase of the DCC process, the solution of a heat transfer and solidification general benchmark problem and the solution of a simplified model of the start of the DCC process in axisymmetry are reported and compared with the published numerical and theoretical data [19,22].

Solidification in an infinite corner
This transient heat transfer and solidification general benchmark problem corresponds to an infinite corner of liquid that starts to freeze under a Dirichlet boundary temperature condition. This example was selected since it allows to validate the proposed FPM formulation against the results of an analytically solvable two-dimensional solidification problem [29]. In this example, an infinite corner of liquid with a uniform initial temperature T 0 =273.45 K starts to freeze under a Dirichlet boundary temperature condition T D =272.15 K on the lower left corner, and in the rest of the edges a Neumann outflow boundary condition was used without taking into account the convective fluid flow generated during the cooling. The melting temperature of the liquid in the corner is T m =273. 15  The analytical solution for the solidification front location assuming constant thermal properties and densities is given by Stapor [29] and it reads The solidification front location at t = 0.5 s obtained by means of this FPM formulation is shown and compared with the previous analytical solution, and with the XFEM simulation results of Stapor [29] in Figure 1. In this case the interphase position calculated with FPM perfectly matches the analytical solution and the Stapor's XFEM simulation with a fine mesh. Further, it could be observed that the solidification front position predicted by FPM with a coarser point cloud is closer to the analytical solution than the corresponding numerical solution obtained with XFEM. These results show that this FPM formulation works very well for the computation of heat transfer and solidification processes.

Simplified model of the start-up phase in a DCC process
Once the FPM formulation was tested for the previous example of heat transfer and solidification, the next natural step is to model and obtain a numerical approximation of the thermal effects involved in the start-up phase of a DCC process. With this purpose in mind, the simplified model of the start of the DCC process in axisymmetry studied by Vertnik et al. in [19] was selected.
The initial configuration of this problem is a cylinder, 0≤z≤0.01m, 0≤r≤0.25m, whose initial temperature is T 0 . The boundary conditions on z=0 are of Dirichlet type with T|(z=0)=T 0 , and the boundaries at the moving bottom are isolated. The boundary conditions at the outer surface are of the Robin type with a reference temperature T ∞ .
In this study case the convective heat transfer coefficients are distributed as shown in Table 1 and the material properties are specified in Table 2. In this formulation, the effect of the latent heat was taken into account using the effective heat capacity method: c=c s f s +c L (1-f s )+h f (∂f s )∕∂T. The liquid fraction increases linearly with temperature between T L and T s . The thermal conductivity in the mushy zone vary linearly with temperature as k m =k s f s +k L (1-f s ), where f s is the solid fraction and it is defined as (21) The smoothing length used in this simulation was h = 0.00875 m with a time step ∆t =0.1 s. The solution of this example has been obtained with an initial discretized domain of 505 nodes with an average spacing of 0.0025 m. In this example, the domain growth is imposed moving the current point cloud according with the casting speed, and adding a new row of nodes in z=0 when a gap with the size of the initial mean spacing is formed.   The temperature profiles over the billet at different time steps predicted by this formulation of FPM are shown in Figure 2 where a comparison with the corresponding numerical results of Vertnik et al. [19] is shown. As it can be seen in this picture, the temperature distributions, and the liquidus and solidus isotherms predicted by FPM match very well with their numerical counterparts in [19].
Further, these figures depict smooth and physical temperature fields with a stable evolution through time. This indicates the FPM potential for the numerical simulation of this start-up DCC transient heat transfer problem since the accuracy of the solutions is appropriate, and the non-linear aspects related to material domain growth and the phase changes are well reproduced.
In Figure 3, the centerline, mid-radius, and surface temperatures at different time steps predicted by this formulation of FPM together with the corresponding numerical results of Vertnik et al. [19] are shown. Regarding the computed temperature along these lines, the graphs are in a very good agreement. However, minor differences in the temperatures up to around 7 K can be observed in some points which are directly attributed to the differences in the numerical approaches. Further, in the FPM solution the domain growth was imposed moving the whole current point cloud according with the casting speed, whilst in [19] the domain growth was imposed moving only the boundary points at the end of the billet according with the casting speed, and adding a new row of points between the boundary and inner ones when it is needed. Nevertheless, these results show the effectiveness of this approach to model the transient thermal behavior and the evolution of the solidified shell thickness in the start-up phase in DCC.

conclusions
Once the current formulation of FPM was implemented and tested in a simplified model of the start-up phase in a DCC process, the following could be concluded: 1. This approach could be used properly to simulate this kind of transient heat transfer and solidification problems in DCC. 2. This approach has shown an excellent behavior for the numerical simulation of these transient problems.
3. This is the first time, to the author's knowledge, that the discussed version of FPM for transient heat transfer and solidification in DCC has been successfully tested for this industrial process.
4. The range of application of FPM has been extended in this work in the context of DCC.
5. Since this formulation is a truly meshless method, there is no need to keep a regular node distribution in order to obtain good numerical approximating solutions, not even to compute any numerical quadrature.
6. It can be used in the future for the analysis of more complex problems in the start-up phase in DCC as the coupling of transient heat transfer and fluid flow.
7. FPM is promising since it is a feasible and a much simpler approach for implementation; it is able to handle high deformations in the domain, and finally, it is capable of naturally incorporating any form of boundary condition without requiring some stabilization or special treatments.