A computer simulation of solidification taking into account the movement of the liquid phase

In the paper, we present the results of solidification simulation taking into account the movement of the liquid phase. The results are obtained from an author software which is implemented on the base of a stabilized finite elements method (Petrov-Galerkin formulation). Using that formulation the Navier-Stokes equation is solved together with the convection term (Boussinesq approximation). The Finite Element Method (FEM) formulation is responsible for solidification, approximating the solution of the heat conduction equation (with the internal heat source term responsible for the heat released during the phase transition). The movement of the liquid phase in a solidifying cast that is caused by convection can significantly affect the process of heat transfer from the casting to the mold, which in turn has an influence on the temperature distribution in the cast and may cause a change in the location of the defects. The presented results allow to assess under what conditions the effect of convection on the solidification process is significant.


Introduction
Solidification is still one of the most difficult process for modelling.Many physical processes that takes place from micro scale up to macro adds present significant difficulties, when someone tries to prepare suitable computational model.During years, many researchers came with models, that could give reasonable results [1,2].
However, many models had some simplifications and limits.Example of such limit was omitting forces of convection in liquid metal.This can lead to significant differences between temperature map from simulation and temperatures observed in real casting.Additional problem is that without convection it is impossible to model such things like distribution of species in casting or solidification shrinkage, which are important factors of castings quality [3].
However, increasing complexity of models, present problem in their implementation.Nowadays, complex models must be implemented with use of technologies like parallel computers [4], accelerated architectures like GPU [5], or using special organization of calculations [6].

Mathematical model
The governing equation for modelling solidification process is based on heat transfer equation with source term: where Boundary conditions, that are used in model with equation ( 1) are the Newton boundary condition on external sides of mold and contact condition for heat exchange between casting and mold.
With use of the apparent heat capacity formulation it is possible to obtain following form of equation ( 1): where c * is the approximation of the effective heat capacity.There are various methods for obtaining this approximation.Here, the Morgan method is used: where H is enthalpy [J] and t in upper script is time level.Liquid metal in this model is assumed to be Newtonian fluid.This allows to write Navier-Stokes set of equations as: (5 where p is pressure [Pa], μ is viscosity [kg/(m s)], f l is liquid fraction [-] and K ε is the permeability of the mushy zone [m].
This equation needs to be accompanied by appropriate set of initial and boundary conditions.Initial u is set as initial condition and between mold and casting the no-slip condition is used.
The right hand side part of equation ( 4) describes body forces that arose in liquid.The first part is connected with buoyance force that is approximated by Boussinesq formula: where , which in this case was temperature from initial conditions.
The last part of left hand side of equation ( 4) is drag force that appears in mushy zone from interaction between liquid and already solidified metal.This model assumes that solidified part is immovable [7].The permeability of the mushy zone is approximated by the Kozeny-Carman equation: where K 0 is secondary dendrite arm spacing [m] and f l is liquid fraction [-].Liquid fraction and solid fraction are connected with simple relation: where Equations ( 1), ( 4) and ( 5) are discretized using Finite Element Method, which leads to following set of equations: , (12) where superscript S is used to set apart matrices used in solidification equation and T, u, p are vectors of unknown temperature, velocity and pressure.Elements of matrices M s , M, N s , N, K s , K, G T , G and D can be calculated from the following formulas: where α, β are indices of single element matrix, a and b have range up to number of nodes in element, i, j, k have range up to number of dimensions, S are shape functions of finite elements and δ is Kronecker delta.
Equations in this form can be solved with carefully selected finite elements [8].In this work, approach that uses stabilized Finite Element is used [9], which allows to avoid limits imposed by the Ladyzhenskaya-Babuska-Breezi condition.
Special care is needed for treating drag force part in stabilization [10].This work uses approach, that determines values of stabilization coefficient not only by velocity of liquid, but also limits them proportionally to volume of liquid fraction.
In calculations, it is assumed that small time step is used, hence in determination of actual values of material properties, when temperature was needed, temperature from previous time step was used.This strategy allows to treat solidification equation as linear and even solve it independently from Navier-Stokes equations, which gives better overall performance [11].Moreover, such approach makes possible to use lumped mass matrix in solidification equation [11,12].
With those assumptions in mind and using theta scheme for time integration, the final form of equations solved in the model is: where matrices with SUPG and PSPG are terms supplied by stabilization, Δt is time step and θ is parameter determining type of time integration scheme (0 for Euler Backward and 0.5 for Cranck-Nicolson).
This model uses Newton-Raphson approach for linearization of Navier-Stokes equations.

Simulations
Model described in section 2 was implemented as an in-house C++ program.This program uses TalyFEM [13] and PETSc [14] libraries.TalyFEM library implements finite element formulas like calculations of shape functions or Gauss quadrature.PETSc library is used for linear algebra operations.
The results for two cases are presented.All simulations use linear triangle finite elements.First is small casting (square shape with 0.05 [m] side length), the second is larger one (0.5 [m] side length), see Fig. 1.Those two cases allow to asses importance of convection in solidification.

Fig. 1. View of domain shape used in simulations
The boundary conditions utilized following parameters: Newton boundary condition on all sides of mold had the environment temperature equal to 300 [K] and heat exchange coefficient was equal to 10 [W/(m K)].Continuity condition assumed value of 1000 [W/(m K)] for heat transfer through separation layer.
The material properties are summarised in  The first series of results present temperature map in four time moments (12.5 [s], 25.0 [s], 37.5 [s], 50.0 [s]) for small square.It is presented in Fig. 2. The effect of convection is only slightly visible in the beginning of simulation.With increasing time, effects of convection are less visible.However, in time moment 50 [s], location of hot spot can be viewed as minimally moved toward top.In this case symmetry of results is easily observable.
Corresponding pressure and velocity map is presented in Fig. 3.The pressure is presented as colour map, while velocity is presented as black vectors.Their size and density corresponds to the magnitude of velocity.The main observation here is that even in places with mushy zone there are no visible fluctuations of pressure, which shows effectiveness of used stabilization technique.After 25 [s] there is almost no liquid phase movement.
In larger case the effect of convection is better visible.It can be observed as fluctuations of temperatures presented in Fig. 4. In Fig. 5, we can see that during whole simulation time, liquid metal had nonzero velocity.More detailed comparison can be done with help of plots presented in Fig. 6 and Fig. 7, where also results from case without convection are used.Fig. 6 presents cooling curves taken in different places of casting (in centre, top, bottom and close to right wall).When there is no convection, curves from bottom, top and right spot are the same.However, with convection, it can be easily seen, that bottom cools faster then top spot.Fig. 7 shows results of the same analysis, but for the large case.Here differences are even more visible.Without convection, the results from places close to walls are the same.

Table 1 .
Material properties for casting

Table 2 .
Material properties for mold