Numerical model of solidification process of Fe-C alloy taking into account the phenomenon of shrinkage cavity formation

Mathematical and numerical models of Fe-C alloy solidification process based on the finite element method (FEM) are presented in the paper. The phenomenon of the defect called "shrinkage cavity" is introduced to the model. It is very important aspect of presented work which makes possible the prediction of the location and shape of mentioned defect depending on the geometry of the casting and cooling condition of the process. An original computer program using Visual C++ environment has been written in order to simulate described process. The results of computer calculations of the three-dimensional solidification problem of the casting with riser are presented and discussed in details.


Introduction
The alloy's solidification process is accompanied by additional phenomena which are important from a technological point of view. They can be responsible for defects in the final product. Defect of the casting is deviation of weight, size, internal structure, external appearance, physical or mechanical properties from requirements described in the norms for particular casting. Usually observed and one of the most important defects of the casting are macroscopic cavities appearing due to contraction of the material during phase change from liquid to solid state. These defects are called shrinkage cavities and they are often discussed in literature [1][2][3][4]. Full classification of the physical and mechanical imperfections appearing in castings contains more than eighty types.
Shrinkage of the material is observed during the three stages of the casting process. This phenomenon is already observed while cooling the alloy in the liquid state, above the liquidus temperature [4]. Due to the small value of the shrinkage coefficient during cooling of the liquid, the effect of this process is difficult to observe. It does not affect the formation of a shrinkage cavity. The most important from the point of view of this work is the loss of material during solidification. The contraction starts in liquidus temperature and it goes on until it reaches solidus temperature. It is observed because solid phase has greater density than the liquid. Depending on the process of solid phase growth, the resulting loss of the volume causes the formation of shrinkage cavity. The last stage of the shrinkage is observed during cooling process in the solid state. This phenomenon causes significant stresses during interactions between the casting and the mold and may be the reason for the formation of cracks.
Presented paper only considers the shrinkage of the casting during solidification. Numerical modelling of the process of shrinkage cavity formation is widely discussed in literature [5][6][7][8][9].

Mathematical description
Considered system consists of the casting which is the essential part and the riser. Risers are removed during final stage of the process. The goal is to obtain defect free casting and the riser containing shrinkage cavity. The mold is neglected in presented model but its influence on heat transport is modelled by appropriate boundary conditions. Considered system is composed of regions ΩS, ΩL, ΩS+L, ΩSC. (Fig. 1). ΩS is filled with solid phase while ΩL contains liquid. ΩS+L is the mixture of solid and liquid and is called "mushy zone". ΩSC is the region occupied by the shrinkage cavity. Air material properties are assigned to this region.
represents effective specific heat. Thermal conductivity and density in equation (1) are averaged according to following formulas: is the liquid (L), solid (S) or gaseous fraction (A).
Variables fs and fl are the functions of temperature: where TL, TS are the liquidus and solidus temperature respectively. The value of gaseous fraction fA can be 0 or 1 due to the assumption that gas cannot form a mixture with the solid and liquid phases. The value of effective specific heat raises significantly between liquidus and solidus temperatures due to releasing of the latent heat of solidification [10]. The equations presented below include a step change in the value of C within the region of solidification: : (4) where c [J•kg -1 •K -1 ] is the specific heat, L [J•kg -1 ] represents latent heat of solidification. In the region ΩSC parameter C = cA. Equation (1) describes initial-boundary problem, thus it must be supplemented by the appropriate conditions. Initial condition is as follows: The following boundary conditions are also used: where: n is the vector outward normal to the boundary Γext, α [J•s -1 •m -2 •K -1 ] denotes coefficient of the boundary heat transfer, T∞ [K] represents ambient temperature. The plane Γsym is thermally insulated because it is the plane of symmetry of solidifying volume and only half of it is taken into account.
According to the method of weighted residuals equation (1) is multiplied by the test function w and integrated over the volume Ω=ΩS ΩL ΩS+L ΩSC: Test function w depends on the chosen numerical formulation of the problem. The following relation is then used: to obtain following equation: To lower the order of the above differential equation Green-Ostrogradsky theorem is used: where Γ = Γext Γsym.
Volume Ω is subjected to spatial discretization to obtain the set of finite elements. Fournode tetrahedrons are chosen for this problem. Galerkin formulation is used to develop local equation for each finite element. Backward Euler method is used for time discretization procedure. Finally global FEM equation is obtained in the following form: where M represents the matrix of heat capacity, K is the matrix of thermal conductivity, B represents the vector of boundary conditions.

Algorithm of shrinkage cavity formation
It is assumed that the process of shrinkage cavity formation depends on gravity. It begins when the solid phase appearing in the analyzed volume. The process of defect formation uses nodes of finite element mesh. Initially each node contains particular volume of liquid material. The sum of nodal volumes equals the total volume of the casting. Then the following operations are done in the each time step of the analysis:  The current volume loss is calculated as the volumetric increment of the solid phase ΔVS [m -3 ] multiplied by the contraction coefficient Sh [-].  Starting from the highest node, the liquid phase is replaced by gas. A volume of gas equal to the calculated shrinkage replaces the liquid in the nodes.  If the shrinkage volume is less than the liquid phase at the node, it goes to the next step. The process ends when the entire casting is filled with solid phase. The assumptions made cause that the shrinkage cavity always begins to form in the upper part of the casting. Then, the expansion of the defect moves down the casting as the liquid level drops.
The computer simulation of the solidification process in the system containing casting with riser is performed. The dimensions of the system together with the boundary and initial conditions are presented in Fig. 2.

Fig. 2. Dimensions (in millimeters), boundary and initial conditions
The entire volume is initially filled with hot, liquid alloy of temperature T0 = 1800 K. The forced and convective motions of the fluid phase are neglected. The plane of symmetry is thermally insulated (qb=0 J•s -1 •m -2 ) while the heat transfer coefficient is defined on the top boundary (α=20 J•s -1 •m -2 •K -1 ) and on the others (α=200 J•s -1 •m -2 •K -1 ). The temperature outside the solidifying system is T∞=300 K. The material properties of Fe-C alloy used in the model are collected in Table 1. Time step Δt is constant and equals 0.005 s. The mesh is built of tetrahedrons and contains 112688 nodes.
Convective boundary conditions are responsible for heat transfer across the external boundaries of the casting and riser. In the first stage of the process the cooling of the liquid without phase change is noticed. There is no shrinkage cavity development during this early step of the process. When the temperature drops below TL the solidification begins and also the shrinkage cavity appears. with an air appears in the upper part of the riser. It evolves downwards following the upper surface of liquid region. The moment at the end of solidification is shown in Fig. 3c. It shows fully developed shrinkage cavity in the form of a funnel. The entire defect is located within the volume of the riser. It is important to choose the shape of the riser and the cooling condition in the appropriate way. It means that shrinkage cavity must not evolve outside the riser into the region of the casting. Such situation could make the final product unusable. Presented results of calculation obtained with the use of original computer program written in Visual C++ environment shows the possibilities of the model in predicting location and shape of the shrinkage cavities in the metal casting process. Appropriate design of the riser and cooling conditions assure the high probability of obtaining casting products of good quality. Accuracy of the obtained results could be higher if the mould is taken into account.