Density Variations during Solidification of Grey Cast Iron

As part of moving towards a sustainable production of diesel engines for heavy vehicle applications, the ability to predict casting defects has become ever so important. In order to model the solidification process for cast components correctly, it is of essence to know how the material will actually behave. To produce sound castings, often of complex geometry, the industry relies on various simulation software for the prediction and avoidance of defects. Thermophysical properties, such as density, play an important part in these simulations.Previous measurements of how the volume of liquid grey iron changes with temperature has been made with a conventional dilatometer. Measurements have also been made in the austenitic range, then on iron-carbon-silicon alloys with a carbon content lower than 1.5 wt%. Based on these measurements the density variations during solidification were calculated. The scope for this paper is to model the volume changes during solidification with the control volume finite difference method, using data from the density measurements.


Introduction
Solidification of cast metals is to a large extent a question of heat transfer. Thus, the software used in casting simulation is most often based on the heat transfer laws. These laws consist of heat balances and, within them, relations of the so called thermophysical properties. One such property is density, which plays an important part for the change of heat content per time unit. Density variations during solidification is also regarded as one of the factors influencing the formation of shrinkage porosity. [1] However, data on liquid density variation of cast iron is scarce and quite widespread over time. Measurements of liquid density have also been carried out with a number of different methods. [2] The measurement methods chosen in this work has been a conventional pushrod dilatometer from Netszch, used for 1D-measurements, and an instrument developed to measure the volume changes in all three directions.
With the increasing use of simulation software to predict casting outcomes it has become more and more important to base the calculations on reliable data, to obtain reliable results, since the accuracy of these data determines the accuracy of the whole simulation [3,4]. A thorough knowledge of density variations will enable simulation of heat conduction, solidification, elasticplastic deformation and fluid flow, which in turn will improve product quality [5].
The aim of the present investigation has therefore been to build a model in order to simulate volume change based on heat distribution in cylindrical and spherical geometries, and to compare these results to data from measurements performed on such bodies. The purpose being to enable further modelling of the volume change during the solidification process in the mushy zone.

Theoretical Background
Heat transfer can occur in three ways: conduction, convection and radiation. [6] Since the focus of this work will be on heat transfer by conduction, the other modes, although of importance in casting simulation, will not be further discussed.
Heat transfer by conduction. Conduction is best described as a heat flow from the high temperature region to the low temperature region. The heat flow per unit area is proportional to the temperature gradient, with the thermal conductivity as proportionality constant. [6] The defining equation for thermal conductivity is expressed as also known as Fourier's law. However, heat will not flow unhindered, therefore thermal resistances must be introduced. If the rate of heat transfer is considered as a flow and the factors thermal conductivity, thickness of the material and the area through which the heat transfer occur, together will form a resistance to this flow, the heat flux can be described as: which can be regarded as an analogy to Ohm's law.
For the use of numerical methods, such as the Control Volume Finite Difference Method (CV-FDM), control volumes are created, over which the heat balance is written.
Heat distribution in cylindrical and spherical samples. For the purpose of simulation of heat distribution in cylindrical and spherical samples of grey cast iron from the dilatation experiments described above, polar cylindrical and spherical co-ordinates would present an advantageous way of reaching the solution. Equations describing the heat transfer in these geometries, using only the radius as descriptive space parameter, will be developed further below. For a more detailed description, please turn to refs. [6] and [7] .
Modelling. Numerical modelling of casting processes often employs the Control Volume Finite Difference Method [8], also known as the Finite Volume Method (FVM), or the Finite Element Method (FEM) [9][10][11]. (There are also models based on the Phase Field, Cellular Automaton and CALPHAD [12] methods, which could be used in combination with the earlier stated methods.)

Method
Cylindrical polar coordinates. The control volumes in the cylindrical polar co-ordinate form would be cylinder caps of thickness Δr i and the distance from the center of the cylinder to point "i" would be r i , ( fig. 1). The heat fluxes over the control volumes are expressed as in eq. 2 and the change of heat content per time unit will be the sum of the heat fluxes and the heat generation in the control volume. Change of heat content per time unit can also be expressed as

Science and Processing of Cast Iron XI
Combining these expressions and discretizing the time derivative will result in eq. 4 ( ) In the case of cylindrical co-ordinates, the size of the two resistance terms in the control volume is different. If the height of the cylinder is Δz and the resistances are inserted, eq. 4 will be expressed as: The transmission resistance terms on the interface are however only used in case of heat transfer between different materials.
Spherical polar co-ordinates. The control volume in the spherical system will be spherical shells, however, the thickness of the shell will be Δr i and the distance from the center of the sphere to point "i" will be r i , just as in the cylindrical system.
The numerical definition of heat conduction in spherical co-ordinates will be the same as in the cylindrical case, i.e. eq. 4. However, the volume will be different, as expressed by eq. 6, as well as the resistance terms. [7] This will yield a different governing equation: The implicit method. There are different ways to solve the governing equation, however, the implicit method has turned out to be stable in most cases. In order to solve the equation numerically it has to be expressed in a form that will suit the application. [6] The capacity and conductivity functions are therefore introduced.
The capacity function, cylindrical co-ordinates, is stated in eq. 8.
The conductivity function, cylindrical co-ordinates, is stated in eq. 9.
The capacity function, spherical co-ordinates, is stated in eq. 10.
The conductivity function, spherical co-ordinates, is stated in eq. 11.
"Time marching." Supposing the temperature field at time t is known, the temperatures at time t+Δt could be solved by the governing equation together with the boundary conditions formulated for points 1 and n. [6] The equation for the numerical solution by the implicit method will have three unknown temperatures to solve, here put on the left-hand side of the equation: in case of cylindrical co-ordinates: = 2 ∆ , (13) and in case of spherical co-ordinates: = 4 .
The solution will be reached by collecting the equations in an equation system, written in matrix form. For the nodes 2 to n-1, eq. 12 is applied and for the outer nodes equations formulated with the boundary conditions will be used. The equation system will be tridiagonal, best solved by the tridiagonal matrix algorithm procedure.
Both heat capacity and heat conductivity are temperature dependent, since the material parameters, among them the density, are temperature dependent. They should therefore change with each temperature.
Enthalpy change. During solidification, release of a certain amount of heat will occur. The way to model this and to ensure that all of the latent heat is released, is with the help of a push-back algorithm [6], ensuring that the temperature step is not to large at the beginning and the end of solidification (and thereby omit the increased portion of the c p value, or using a too large c p , respectively).
The amount of heat released could be written: By differentiating eq. 15, rewriting and inserting it into the heat conduction equation we get the expression: This yields that in the solidification interval the c p -value must be adjusted:

Experimental Work
Volume change and density measurements. For the 1D-measurements a conventional push-rod dilatometer from Netzsch (DIL 402C) was used. Measurement method and results have been described in depth elsewhere [13] and will here only be explained in brief. Grey iron alloys of varying carbon contents were cast and specimens in the form of small cylinders (weighing approximately 2.8 g) were machined. These samples were aimed for testing in the liquid state, in a range from approximately 100 °C above liquidus temperature down to liquidus temperature. The samples were first encased in an aluminum oxide cylinder with a piston at each end, transferring the movements of the sample to the push rod. Data could thus be logged for length change at the corresponding temperature. In order to measure the contraction behavior of austenite, a series of steel alloys with varying carbon content were cast, and likewise machined into cylinders. These samples were tested in the temperature range where the material would be austenitic. Testing of these samples was performed without the aluminum oxide container. The calculations of the density have been based on the assumption of material isotropy. Thus, the axial expansion could be set as proportional to the expansion of the diameter. The density of a cylindrical sample (at an arbitrary temperature) can be written as: The density calculations for the liquid samples build on the same assumption, however, in the case of radial expansion, these samples were confined by the container. The diameter was therefore set to the inner diameter of the aluminum oxide container at the same temperature. It may here also be pointed out that because of the position of the thermocouple (near the outer wall of the aluminum 158 Science and Processing of Cast Iron XI oxide container) the sample temperature will deviate from the temperature logged by the instrument. The temperature difference depends largely on the thermal properties of the aluminum oxide.
For the 3D-measurements, a spherical vessel filled with liquid grey cast iron was used. [14] The volume change of the sphere was measured at three points at the surface and an average of the change in radial direction was calculated at each timestep.
The chemical composition of the alloys used in the simulation can be found in table 1. Simulation. Two different scripts, written and executed in MATLAB and based on the described calculation method, were used, one simulating the cylindrical case and the other simulating the spherical case. Both, however, were structured in the same way, calculating the temperature distribution by equation 12 and taking the enthalpy change into account. In the case of the cylindrical sample, a uniform mesh of Δr = 0.325 mm was used. The mesh contained 10 metal cells, six cells in the mold and two for the surrounding atmosphere. In order to simulate the temperature distribution, chosen values were applied. These can be found in table 2. For the simulation of the spherical sample a mesh of Δr =0.5 mm was applied. The mesh contained 33 metal cells, one for the mold and two for the surrounding air. Simulating the cylindrical case, the experimentally measured density data in the liquid and the austenitic ranges were read into the program and the temperature in the first surrounding-air node was put equal to the temperature measured in the experiment. The volume of the simulated cylindrical sample was calculated by: at each timestep, based on the temperature in the metal cells, the weight of the sample (measured after the experiment) and the density at the temperature in question. The simulated volume was then compared with the volume curve from the experiment. For the spherical samples the value for the weight was chosen to 145 g. The calculation of the volume in the liquid and austenitic ranges was performed in the same way as for the cylindrical samples, apart from that the temperature field was assigned to the mold cell.

Results
The results of the simulation of the cylindrical samples are presented in figures 2, 3 and 4. Since the result from the dilatometer experiments on this particular alloy showed two different behaviors, both types has been used in the comparison (samples C325-7 and C325-20).   The simulated volume and the measured volume for the spherical sample is presented in figures 5 and 6. The volume has been calculated based on the data from the same two density measurements as in figure 2, 3 and 4.

Discussion
The discrepancy between the simulated volume and the volume curve from the measurement in pictures 3 to 7 could be found in the measurement of the temperature in the Netzsch dilatometer. Since the thermocouple is placed outside of the container and the measurements were performed on a cooling liquid, the temperature in the liquid metal would be higher than the thermocouple is logging. The magnitude of the temperature difference should therefore be investigated and taken into account in the further modelling. The other modes of heat transfer have not been taken into account, since it was not the scope at this stage, but these may also have an impact on the results.
In the austenitic range, in all of the pictures, the simulated curve follows the contraction of austenite with a carbon content of 0.59 wt% carbon. Segregation effects has thus not been taken into account. The graphite fraction, which in the actual measurements has got an influence, has likewise not been modelled. Further modelling should therefore include these areas. Based on this work, the volume changes in the mushy zone will be modelled, thereby obtaining a deeper understanding of the mechanisms behind shrinkage porosity.

Conclusion
It is possible to model and simulate volume change based on temperature distribution and measured density data with good results.
It is likely that the addition of a graphite fraction to the curve in the austenitic range will position the simulated curve closer to the measured one.
The present investigation constitutes a good basis for further modelling of volume changes during the earlier stages of solidification.