Implementation and Validation of Casting Simulation Methodology for Diagnostics of Lamellar Graphite Iron

This paper describes and validates a methodology for implementation of full-scale sand-casting simulation in a general-purpose finite element software, including mold filling, heat transport, solidification kinetics, chemical microsegregation and prediction of microstructure and material properties. The solidification model, customized for gray cast iron, includes novel methods for handling interaction between parallel dendritic and eutectic solidification modes and its impact of their interaction on the final microstructure. The validation involves a previously published gray iron casting experiment and involves comparison of simulated and experimental cooling curves, microstructure parameters and tensile strength. We believe that this is valuable to researchers and engineers seeking to improve the state of the art of casting simulation tools.


Introduction
Casting simulation is a valuable tool used by product developers and foundry engineers for product design, problem solving and optimization of the casting process. Improvement of simulation capabilities often builds on published research. Yet, progress is to some degree held back by the black-box nature of many commercially available casting simulation software, making it challenging to identify flaws and remedies and test new ideas. For this reason, it is essential that methodology is available for how to test new models in an open casting simulation setting. COMSOL Multiphysics is a general-purpose finite element (FE) simulation software with fully coupled multiphysics modeling capabilities and possibility for users to extend the standard models with custom differential equations. 1 In this paper, we develop and implement a COMSOL-based methodology to simulate key aspects of an iron casting process, such as mold filling, heat conduction, solidification kinetics in relation to a ternary equilibrium phase diagram, microstructure and material properties. A previously published model for tensile strength in pearlitic gray iron 2 was complemented by a new model for solidification kinetics, aiming to eliminate reliance on experimental cooling curve data. The formulated model is a micro/macromodel suitable for gray cast irons, including novel features such as restriction of eutectic to the inside of dendrite envelopes and a method for predicting the dendritic structure observed in the as-cast material based on its evolution in contact with the liquid and gradual engulfment by the eutectic. The solidification model was first implemented using the numeric computing environment MATLAB in one dimension (1-D). 3 Once all the features had been implemented and proven stable, it was implemented in COMSOL Multiphysics for simulation in three dimensions (3-D). Agreement between the 1-D and 3-D implementations was verified by comparing temperature curves on a simple spherical geometry. The 3-D implementation was then applied to previously published gray iron casting experiment data 2 to validate its applicability to a complex casting scenario, involving mold filling. Comparison between simulation and experiment involved cooling curves, microstructure parameters and tensile strength We believe this work is valuable for researchers and engineers interested in advancing the capabilities of casting simulation tools, both with respect to the outlined methodology for casting simulation in a general finite element software and with respect to the novel features of the formulated solidification and microstructure model for gray iron.

Simulation Methodology for Casting Diagnostics
The proposed simulation methodology is designed for gravity casting of lamellar graphite iron, to enable prediction of both the evolving features of the solidification structure as well as the final microstructure features found in the cast material. The evolving features of the solidification structure are important for ease of material transport during solidification, while the final microstructure features are important for prediction of local material properties. This work focuses on prediction of tensile strength using a published model utilizing a microstructure parameter called the hydraulic diameter of interdendritic phases. 2,4 However, extension with additional microstructure parameters and material properties should be easy following the same methodology.
The fully user-controlled workflow is overviewed in Figure 1. It begins with importing the computer-aided design (CAD) geometry of the casting and the mold, and simulation of the mold filling, where solidification effects are assumed negligible due to relatively high pouring temperature. The temperature field in casting and mold at the end of the mold filling are then transferred onto the initial temperature field (T init ) for the simulation of solidification, microstructure and property formation. It is also possible to import the initial temperature field from other mold filling simulation tools; however, an additional importing interface would be required between the software tools. The initial temperature field is best validated using experimental cooling curves or verified using more advanced Finite Volume Method (FVM) CFD software. The solidification model described in the following section is a key component of the simulation methodology.

Solidification Model
A brief overview of the solidification model customized for gray cast iron is presented here. The numerical scheme is a micro/macromethod which can be characterized as a fixeddomain source-term method. 5 In essence, a mathematical model of growing microscopic solidification units defines the rate of phase transformations which latent heat release is placed in the source term of the heat transfer equation. The heat transfer equation is solved numerically over the macroscopic mesh. Besides defining the source term, the micro-model involves modeling of chemical microsegregation, phase composition and representations of the solidification microstructure which provides important information about the cast material. While much of the micro-model is adopted from predecessors, the two more novel features are the gradual capture of evolving dendritic structure into the eutectic and the restriction of eutectic to the inside of dendrite envelopes.
As the melt undercools below the liquidus temperature, the solidification model is activated. Phase transformations are driven by undercooling in relation to a description of the Fe-C-Si equilibrium phase diagram, obtained by regression to phase equilibria calculated using Thermo-Calc 2019a with the TCFE7 database. 6

Growth of Primary and Eutectic Envelopes
Phases considered in the solidification model are liquid (L), austenite (c) and graphite (G). During the primary stage of solidification, randomly distributed spherical envelopes grow and impinge in accordance with the Kolmogorov-Johnson-Mehl-Avrami (KJMA) model, 7 assuming a constant number density of nucleation sites N D . The volume fraction of primary envelopes is given by Eqn. (1).
While the argument of the natural exponential function (exp) is the so-called extended volume fraction including overlap of the randomly distributed spheres, g D is the corrected volume fraction without the overlap. Once the temperature falls below the eutectic temperature, a new set of spherical envelopes is introduced by Eqn. (2), growing and impinging in the same manner as the primary envelopes.
The rate of radial growth of both envelopes described by Eqns. (3) and (4) is according to the Oldfield model, 8 though with separate coefficients K D and K E and exponents n D and n E for primary and eutectic envelopes, respectively.
The eutectic envelopes are superposed on the primary envelopes. The portion of primary envelopes not intersecting with eutectic (D\E) contains a phase mixture of austenite and liquid (c ? L), while the intersection between the two (D \ E) contains the eutectic mixture of austenite and graphite (c ? G). In other words, free dendrites remain only in the portion of the primary envelopes not yet overgrown by eutectic. Likewise, eutectic exists only in eutectic envelopes which are inside primary envelopes. The interaction between the growing primary and eutectic domains is illustrated schematically in Figure 2. The restriction of eutectic to the inside of primary envelopes is motivated by observations that the austenite in eutectic cells tends to have the same crystallographic orientation as the surrounding dendrite grains, indicating that eutectic austenites are merely extensions of the earlier dendrites. 9,10 Phase Composition and Microsegregation The total volume fraction of austenite is thus described as Eqn. (5).
The factor g D 1 À g E ð Þrestricts the internal austenite g iA to the portion of the primary envelopes which are not occupied by eutectic envelopes. Likewise, the factor g D g E restricts eutectic austenite to the eutectic envelopes inside primary envelopes. The total volume fraction of graphite is described in the same manner as Eqn. (6).
The volume occupied neither by austenite nor graphite remains liquid. Both envelopes adopt the average chemical composition (w 0 C and w 0 Si ) from the parent liquid; however, solute partitions between the phases. The partition coefficient k Si is updated in accordance with the phase diagram description based on the chemical composition of the interdendritic liquid (w LA C , w iL Si ). Microsegregation of Si is modeled by Eqn. (7) according to the differential form of the Gulliver-Scheil equation. 11,12 ow iL Here, g LD refers to the volume fraction of liquid in primary envelopes, corresponding to the volume fraction of the envelope occupied by neither primary austenite nor eutectic, expressed as Eqn. (8).
Microsegregation of other substitutional alloying elements can be modeled in the same manner, assuming that partition coefficients and initial concentrations are known. It follows from microsegregation that eutectic austenite may have a different average concentration of Si than primary austenite. A method for calculating the average in the eutectic is treated in section Capture of parameters in the eutectic.
The volume fraction of austenite in the primary envelope is calculated using a mass balance of carbon, where concentrations of liquid and austenite are taken from the phase diagram description for the local temperature and concentration of Si in the liquid. Accounting for unequal densities of the two phases, the mass balance leads to expression Eqn. (9) The equation applies also below the eutectic temperature by extrapolation of w LA C and k C to lower temperatures. This implies that the interdendritic liquid remains in equilibrium with austenite below the eutectic temperature. This assumption has proven useful for modeling of ductile iron, where liquid is rarely in direct contact with graphite. 13,14 However, it seems to be a reasonable approximation also for hypoeutectic gray iron, considering that the liquid is also here predominantly in contact with austenite. Primary precipitation of graphite is outside the scope of the present paper. The mass fraction of graphite in the eutectic envelope can be calculated in a similar manner as Eqn. (10).
However, graphite is here assumed to be pure carbon and the concentration of C in austenite, w AE C , is driven toward equilibrium with graphite by diffusion of carbon toward the graphite in accordance with the differential Eqn. (11).
For simplicity, diffusion is assumed to occur in one direction and the diffusion length is assumed to be the mean dendrite modulus captured in the eutectic M (more details in sections Dendrite morphology during solidification and Capture of parameters in eutectic).
Modeling of the desaturation of austenite is a prerequisite for extensions of the simulation to include the solid-state transformations of austenite and associated microstructure features.

Dendrite Morphology During Solidification
Austenite in the primary envelope grows with dendritic morphology, which is known to coarsen while in contact with liquid. The scale of dendrite structure is here characterized by its modulus (M), defined as the ratio of its volume over its surface area. 15 Interpretation of M is facilitated by the fact that it approaches half the radius of a long cylinder, a geometry comparable to dendrite arms. But the modulus has the advantage of being defined for any shape. The modulus of dendrites in gray iron has been shown to evolve in proportion to the cubic root of solidification time. 15 Details about the new coarsening model for the modulus, represented here by Eqn. (12), will be published elsewhere.
oM ot ¼ f ðt; M; g iA ; g D Þ Eqn: 12 As indicated by arguments g iA and g D of the function, the model draws benefit from the present solidification model formulation.

Capture of Parameters in the Eutectic
The dendrite modulus M continues to evolve until the dendrite is captured in the eutectic envelope. Consequently, the gradual capture of dendrites by growing eutectic causes spatial variation of M at the scale of eutectic envelopes. For many purposes, the mean modulus M captured in the eutectic is a more practical parameter, because it represents better what would be measured on a cross section of the material. The parameter M can be calculated using Eqn. (13) which is solved along with the other differential equations.
MðtÞ À M t ð Þ g E Eqn: 13 Equation (13) can be interpreted as dilution of the mean modulus M, as M of free dendrites is captured at the boundary of the growing eutectic g E . The captured mean of other parameters, like g iA and w AE Si , is calculated in the same manner by substitution of corresponding parameters for M and M in Eqn. (13). The mean of parameters which have no interaction with the solver is best calculated in post processing.

Model for Tensile Strength
This work employs an empirical Eqn. (14), which relates the tensile strength of pearlitic gray cast iron with the morphology of the dendritic structure. 4 Note that the units of r UTS , D IP hyd and k t are here different from, 4 meaning that scaling of k t is necessary.
The dendrite morphology is here represented by the hydraulic diameter of interdendritic phases, which is a shape-independent measure of the scale of the interdendritic space, calculated as the ratio of its volume to its surface area. In essence, it is almost identical to M, the only difference being that it characterizes the space between dendrite arms rather than the arms themselves. While the correlation of the parameter with tensile strength is an empirical finding, it may play a mechanistic role as a measure of the space in which crack-initiating graphite flakes may exist. As this parameter is measured over an area of a cross section of the solid material, it corresponds to the mean of the parameter captured in the eutectic, D IP hyd . This parameter is easily calculated from g iA and M using Eqn. (15) due to its direct geometric relation with the two.

Implementation of a 1-D Solidification Model
The solidification model was first implemented into a 1-D heat conduction simulation in MATLAB using implicit control volume-based finite difference method (implicit CV-FDM) with spherical symmetry. The numerical treatment of the heat conduction calculation has been described in detail elsewhere. 16 The rate of release of latent heat of fusion was calculated by means of Eqn. (16) as the sum of contributions from the two solid phases.
The derivatives og A =ot and og G =ot were obtained by differentiation of Eqns. (9) and (10). The resulting expressions are extensive and are for this reason not shown here.
Including all released latent heat of fusion in the source term made the simulation prone to numerical instability. To counter this issue, the portion of the latent heat related to the regulation of austenite in the primary envelope g iA was excluded from the source and instead accounted for in an apparent heat capacity C Ã The derivative og iA =oT was obtained from Eqn. (9) using the slope of the liquidus surface in the phase diagram. The term including og iA =ot must consequently be excluded from og A =ot in Eqn. (16).
This would make the numerical scheme used for the 1-D implementation a hybrid of the source-term method and the apparent heat capacity method (also known as effective specific heat method). Note that this was not used in the latter 3-D implementation of the solidification model into the FE software.

Implementation of a Full-Scale Iron Casting Process Simulation
The model for the mold filling simulation overviewed in Figure 3 was realized in COMSOL using its CFD and heat transfer modules. Heat transfer equations were solved both for the mold material and for the liquid iron. The coupling between the two was obtained via the interfacial heat flux determined by iron-to-mold heat transfer coefficients. The molten iron flow was assumed laminar creeping flow, and it was modeled as a two-phase flow with the moving boundary between the liquid iron and the air evaluated by level set method, with Navier slip condition set at the casting walls. While these simplifications are sufficient to provide a reasonable initial temperature field T init , which is the scope in this context, more advanced modeling is required to predict flow induced casting defects such as gas entrapment. The velocity field u was frozen after the mold filling, such that the solidification model did not include any fluid dynamics. The temperature distribution (T init ) at the end of mold filling was transferred to the heat transfer equations both for the mold and metal domains. While simulation of convection in parallel with solidification was beyond the scope of this work, such an expansion of the simulation would improve the accuracy of the temperature field and account for the effect of macroscopic chemical segregation.
The model for heat transfer during solidification in the metal domain is overviewed in Figure 4  Implementation of the 3-D solidification model is rather straight forward, as the numerical treatment is to a high degree handled by the software. The solidification model was implemented using the Heat Transfer in Solids interface (COMSOL Heat Transfer module). Note that, contrary to the 1-D implementation, all latent heat of fusion was here introduced into the heat conduction equation as source.

Heat Flux on the Metal-Mold Interface
Since the solidification model includes only the solid domains (mold and metal), the heat flux at the metal-mold interfaces is set up differently than in the case of convective fluid-solid heat transfer. In case of metal casting, the temperature difference controls the heat flux at the shared boundary between the sand and the metal. In casting simulation software, the heat transfer coefficient between the metal and sand domains is specified as the function of metal temperature. This is achieved using the Thermal Contact feature. Assuming for simplicity one-dimensional heat flux, the Thermal Contact feature expresses heat fluxes q d and q u at the boundary, as Eqns. (18) and (19).
Eqn: 19 where h (W m -2 K -1 ) is the reciprocal of the contact resistance R eq (Km 2 W -1 ).       Subscripts d and u refer to the down-and upsides of the boundary. Which is the upside of the boundary is determined internally by the software when the geometry is set up, meaning the metal may end up on either side. To ensure that T metal points to the temperature on the metal side of boundaries, domain indices can be determined graphically using operators up(dom) and down(dom), where ''dom'' is the domain index. For example, application of operator up(dom) to the boundary in Figure 5 returns domain index 2, which in this example corresponds to the metal domain, meaning T metal = up(T) in Eqns. (18) and (19) on this boundary.

One-Dimensional Model Verification: Spherical Casting
Agreement between the 3-D and 1-D solidification model implementations was verified by comparison with cooling curves taken from two locations of a spherical geometry, as shown in Figure 6. Molten gray iron was solidified inside a thin-walled spherical shell, having the inner radius 22 mm. Utilizing symmetries of the sphere, the geometry was simplified to a thin slice to speed up the FE simulation. Temperature curves from the center and the inner surface of the sphere were compared, as shown in Figure 7. The curves are almost perfectly overlapping. The residual small difference is likely a consequence of the different numerical treatments in the MATLAB and COMSOL implementations.

Three-Dimensional Model Validation: Cylindrical Castings
To validate the simulation methodology for a more complex casting geometry, as well as to demonstrate how it may combine with a mold filling simulation, the implemented 3-D model was applied to a casting experiment. 2,4 The geometry (Figure 8), process parameters, probe locations and the experimental cooling curves were adopted from. 2 The meshed finite element geometry is shown in Figure 9, where a finer mesh is applied on the casting and insulation and chill cylinders. The temperature distribution on a cross section after the mold filling is shown in Figure 10. Figure 11 shows a comparison between the simulated and experimental cooling curves at the probe locations within insulated and sand-encapsulated cast cylinders, which both show good agreement. The chilled cylinder is excluded from the comparison based on its carbide content, given that the solidification model does not account for growth of the metastable carbide-austenite eutectic. Table 1 shows agreement between simulation and experiment with regard to microstructural parameters and the UTS for the two cooling rates. All the % errors fall within 10% of the corresponding measured quantity, which can be regarded good agreement given the complexity and multiphysical nature of the solved problem.
The agreement between simulation and experiment in Figure 11 and Table 1 must be interpreted in light of the adjustments of heat transfer coefficients and parameters of the solidification model.
Heat transfer coefficients between the metal and the other domains (insulation, sand, chill) were adjusted with a focus on agreement of cooling curves for temperatures outside the solidification range. Heat transfer coefficients constitute a major uncertainty in casting simulations, varying significantly in their implementation between different commercial simulation tools.
Parameters related to solidification kinetics were subsequently adjusted within reasonable limits such as to reproduce the observed recalescence with good agreement. While growth law parameters were global, the number of nucleation sites was fitted separately for each cylinder. The individual fit of the number of nucleation sites is motivated given that it is known to vary considerably depending on cooling conditions. Obviously, in order to make the presented simulation approach useful, it must be complemented with a nucleation law which predicts the number of nucleation sites based on the local cooling conditions and the metallurgy of the metal.
With the uncertainties coming with these adjustments in mind, the agreement between simulated and experimental cooling curves in Figure 11 and the parallel agreement of microstructure parameters in Table 1 demonstrates (a) that the presented simulation approach is feasible for testing of new models for prediction of solidification kinetics, microstructure and material properties and (b) that the introduced solidification model is capable of predicting realistic microstructure parameters and tensile strength given realistic temperature evolution and (c) that the solidification model is capable of generating realistic cooling curves.
However, more work is needed to determine appropriate solidification model parameters more rigorously and to introduce a responsive nucleation law.
The introduced model for solidification kinetics and prediction of microstructure parameters complements the previously validated tensile strength model in that it relies on input about chemical composition and parameters related to growth kinetics, rather than experimentally determined apparent specific heat capacity (C p *) used in its previous implementation in Flow-3D. 2

Conclusions
A methodology was presented and validated for simulation of key aspects of gravity casting process in a general FEA software, including mold filling, heat transfer in mold and metal, solidification and prediction of local microstructure and material properties. The methodology is valuable for researchers who seek to advance the capabilities of casting simulation and opens an alternative route for empowering users of casting simulation tools with the latest models. The implemented solidification model for gray cast iron includes unique features which are valuable for improved accuracy of solidification kinetics and prediction of microstructure and properties. The main contributions are highlighted below: • Casting process simulation model equations, structure and data flow were presented to enable an open iron casting simulation setting. • Practical guidelines were provided on how to implement heat transfer coefficients at boundaries between mold and metal in the FEA software. • The casting process simulation methodology was validated by application to a complex casting experiment from a previous publication, demonstrating that it is feasible for testing of new models for solidification, microstructure and properties. • The solidification model handles interaction between primary and eutectic solidification units in two novel ways: (1) eutectic only grows inside dendrite envelopes and (2) the as-cast dendritic structure is calculated by gradual capture of the evolving free dendrites into the growing eutectic. • Agreement between simulation and experiment demonstrated that the solidification model provides reasonable microstructure parameters and tensile strength given that cooling curves are reproduced with decent accuracy. Moreover, it demonstrated that the solidification model is capable of generating realistic recalescence. • Connection of the solidification model to a multicomponent (in this case ternary) phase diagram was found useful for estimation of the impact of changes in chemical composition of the alloy on cooling curves, microstructure and properties. • While the solidification model was designed for hypoeutectic gray irons, the phase diagram and growth laws are easily adapted to other hypoeutectic cast irons or other alloys.

List of Symbols Phases and Envelopes
g D Volume fraction of primary envelopes g E Volume fraction of eutectic envelopes g A Total volume fraction of austenite g G Total volume fraction of graphite g iA Internal volume fraction of austenite in the portion of primary envelopes not occupied by eutectic envelopes g iG Internal volume fraction of graphite in the intersection of eutectic and primary envelopes g LD Volume fraction of liquid in primary envelopes g iA Mean internal volume fraction of austenite captured in eutectic r D (m) r E (m) Specific heat capacity of molten iron during mold filling (J kg -1 K -1 ) C p * Apparent specific heat capacity of iron during solidification (J kg -1 K -1 )