In my work, I focus on the numerical simulation of shrinkage defects in ductile iron castings. Ductile iron castings solidify mainly in a pasty manner, and their graphite expansion modifies the feeding behavior of liquid metal in a way that is not captured well by conventional models. I therefore developed a macro-micro simulation framework that couples microstructural evolution, temperature field calculation, isolated feeding zone identification, graphite expansion pressure, and flow pressure loss. The purpose is to predict the location and severity of shrinkage porosity and shrinkage cavities in ductile iron castings more reliably, and to support process optimization without relying only on trial-and-error pouring experiments.
I begin with the thermal problem because the temperature field controls the solid fraction, the formation of isolated liquid regions, the release of latent heat, and the final distribution of defects. In my model, the solid fraction is not assumed to vary linearly with temperature. Instead, I calculate the nucleation and growth of graphite nodules and austenite, and then use a Johnson-Mehl type relation to connect the microscopic phase evolution to the macroscopic solid fraction. This gives a more realistic description of latent heat release in ductile iron castings.
The following table summarizes the main symbols I use throughout my simulation framework.
| Symbol | Meaning | Typical unit |
|---|---|---|
| \(CE\) | Carbon equivalent of ductile iron | % |
| \(T_E\) | Eutectic temperature | °C |
| \(C_{TE}\) | Carbon content in austenite at eutectic temperature | % |
| \(f_s\) | Solid fraction | – |
| \(f_{gr}\) | Graphite fraction | – |
| \(f_{\gamma}\) | Austenite fraction | – |
| \(N_{gr}\) | Graphite nodule nucleation density | m\(^{-3}\) |
| \(R_{gr}\) | Graphite nodule radius | m |
| \(R_{\gamma}\) | Austenite shell radius | m |
| \(\Delta V_{iso}\) | Net volume change of an isolated feeding zone | m\(^3\) |
| \(P_e\) | Graphite expansion pressure | Pa |
| \(\Delta P_s\) | Pressure loss through a porous region | Pa |
| \(G\) | Feeding potential | Pa |
I treat the carbon equivalent as a first-order indicator of whether the alloy is hypoeutectic or hypereutectic. For ductile iron castings, the carbon equivalent is calculated from the major alloying elements:
$$CE = C + 0.33Si + 0.33P + 0.4S – 0.03Mn$$
where the elemental symbols denote mass fractions in percent. When \(CE \gt 4.26\), graphite is the leading phase during the initial stage of solidification. When \(CE \lt 4.26\), austenite is the leading phase. This distinction matters because the sequence of phase precipitation controls the local expansion and contraction behavior in ductile iron castings.
I calculate the eutectic temperature and the carbon content of austenite at the eutectic temperature as follows:
$$T_E = 1154.6 + 6.5Si$$
$$C_{TE} = 2.1 – 0.216Si$$
For the primary stage, I use a simple limiting fraction to represent the amount of primary austenite or primary graphite that can form before the eutectic stage begins. For a hypoeutectic alloy,
$$f_{au,pr} = \frac{4.26 – CE}{4.26 – C_{TE}}, \quad CE \lt 4.26$$
For a hypereutectic alloy,
$$f_{gr,pr} = \frac{CE – 4.26}{100 – 4.26}, \quad CE \gt 4.26$$
During the eutectic stage, I do not use a linear rule. I calculate graphite nucleation and growth, austenite shell growth, and the resulting solid fraction. The graphite nucleation rate is described by
$$\frac{\partial N_{gr}}{\partial t} = (1 – f_s) b \Delta T \exp\left(-\frac{c}{\Delta T}\right)$$
where \(b\) and \(c\) are nucleation parameters that depend on composition and melt treatment. In my simulations, I use representative values such as \(b = 4.0 \times 10^{13}\) nuclei per cubic meter per degree Celsius per second and \(c = 340\,^\circ\mathrm{C}\). The undercooling is \(\Delta T\).
Graphite growth is controlled by carbon diffusion through the austenite shell. I use
$$\frac{dR_{gr}}{dt} = \frac{(C_{\gamma/l} – C_{\gamma/gr}) \rho_{\gamma}}{(C_{gr} – C_{\gamma/gr}) \rho_{gr}} \frac{D_c^{\gamma}}{R_{gr}} \left[1 – \frac{1}{R_{\gamma}/R_{gr}}\right]$$
and the austenite shell growth is
$$\frac{dR_{\gamma}}{dt} = \frac{D_c^{\gamma} (C_{\gamma/l} – C_{\gamma/gr})}{R_{\gamma} (-1 + R_{\gamma}/R_{gr})(C_{l/\gamma} – C_{\gamma/l})}$$
The interfacial carbon concentrations are evaluated from local temperature and silicon content:
$$C_{l/\gamma} = \frac{1}{97.3}(1569 – T – 24.32Si)$$
$$C_{\gamma/l} = \frac{1}{177.9}(1528.4 – T – 32Si)$$
$$C_{\gamma/gr} = \frac{(T – 1154.6 – 6.5Si)(1.5 – 0.216Si)}{354.6 + 6.5Si} + 2.1 – 0.216Si$$
To connect the micro-scale phase evolution to the macro-scale solid fraction, I use a Johnson-Mehl type relation:
$$\frac{\partial f_s}{\partial t} = (1 – f_s) 4\pi N R^2 \frac{\partial R}{\partial t}$$
This equation allows me to obtain the solid fraction increment from the nucleation density and the growth rate of the eutectic cells. In ductile iron castings, this is important because the latent heat release is not uniform; it is concentrated in the eutectic stage, which produces the recalescence and the temperature plateau observed in real cooling curves.
I also track the mass fraction of graphite within the eutectic cell:
$$f_{gr,eu} = \frac{\rho_{gr} R_{gr}^3}{\rho_{\gamma}(R^3 – R_{gr}^3) + \rho_{gr} R_{gr}^3}$$
$$f_{\gamma,eu} = f_s – f_{gr}$$
For latent heat, I use an equivalent specific heat during the primary stage and a source term during the eutectic stage. The equivalent specific heat is
$$c_{eq} = c + \frac{L_{pre}}{\Delta T}$$
and the eutectic heat source is
$$Q = \rho L \frac{d f_s}{d t}$$
The transient heat conduction equation I solve for the temperature field in ductile iron castings is
$$\frac{\partial T}{\partial t} = \frac{k}{\rho c} \left(\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2}\right) + \frac{Q}{\rho c}$$
where \(k\) is thermal conductivity, \(\rho\) is density, and \(c\) is specific heat. I discretize this equation on an orthogonal hexahedral grid using the finite difference method. For each cell, heat exchange occurs with six neighbors. The explicit update can be written as
$$T^{n+1} = T^n + \sum_j A_j (T_j^n – T^n) + B$$
where \(j\) denotes a neighboring cell, \(A_j\) contains the thermal properties, time step, and grid spacing, and \(B\) contains the heat source contribution. To keep the explicit scheme stable, I constrain the time step by a Courant-Friedrichs-Lewy condition:
$$dt \lt \alpha \frac{\rho c \, dx \, dx}{k}$$
where \(\alpha\) is a safety factor less than one. At the domain boundary, I apply a fixed temperature for the outer air layer and a Robin condition for interfaces between different materials. Instead of specifying a highly uncertain interface heat transfer coefficient directly, I use an equivalent thermal conductivity between neighboring cells:
$$k_{eq} = \frac{1}{d_i/k_i + 1/h + d_j/k_j}$$
where \(d_i\) and \(d_j\) are the distances from the cell centers to the interface, \(k_i\) and \(k_j\) are the thermal conductivities, and \(h\) is the interface heat transfer coefficient. This reduces the sensitivity of the simulation to the interface parameter and improves the robustness of the temperature field calculation for ductile iron castings.
I validated the temperature field model with a sample cup experiment. A K-type thermocouple was placed at the center of a small cup-shaped casting, and the measured cooling curve was compared with the simulated curve. The alloy was a ductile iron with the composition shown below.
| Element | C | Si | Mn | S | P | Ti | Al | Mg | Ce | Fe |
|---|---|---|---|---|---|---|---|---|---|---|
| Mass fraction (%) | 3.64 | 2.35 | 0.124 | 0.0163 | 0.0355 | 0.0205 | 0.0202 | 0.0569 | 0.0107 | Balance |
The thermophysical properties used in the initial simulation are listed in the next table.
| Property | Value | Unit |
|---|---|---|
| Density of casting | 6.81521 | g·cm\(^{-3}\) |
| Thermal conductivity of casting | 0.076316 | cal·(cm·s·°C)\(^{-1}\) |
| Specific heat of casting | 0.191537 | cal·(g·°C)\(^{-1}\) |
| Liquidus temperature | 1171.2 | °C |
| Solidus temperature | 1140 | °C |
| Latent heat | 61.2383 | cal·g\(^{-1}\) |
| Initial casting temperature | 1300 | °C |
| Initial mold temperature | 20 | °C |
The measured cooling curve shows a rapid decrease, then a slower decrease, and finally another rapid decrease. The slow stage corresponds to the eutectic reaction, where graphite and austenite release latent heat. The simulated solid fraction also exhibits a nonlinear increase with time. At approximately \(1141\,^\circ\mathrm{C}\), the solid fraction increases while the temperature rises slightly, which reproduces the recalescence phenomenon found in ductile iron castings.
When I first used the initial latent heat value, the simulated temperature fell too quickly after about 175 s, and the temperature plateau was shorter than the measured one. I therefore adjusted the latent heat to \(80\,\mathrm{cal\cdot g^{-1}}\). With this adjustment, the correlation coefficient between the simulated and measured temperature curves increased from 0.93372 to 0.99357. This validation confirms that the microstructurally informed temperature field is sufficiently accurate for subsequent shrinkage prediction in ductile iron castings.
After obtaining a reliable temperature field, I construct the shrinkage defect prediction model. The first step is to identify isolated feeding zones. During solidification, the solid fraction increases. When the solid fraction exceeds a critical value \(f_{cs}\), the dendritic network becomes coherent, and liquid flow is strongly restricted. Regions with \(f_s \lt f_{cs}\) can still feed each other, while regions with \(f_s \ge f_{cs}\) are treated as solid or porous. The calculation domain is divided into three types of regions: liquid, mushy, and solid. An isolated feeding zone consists of liquid and mushy cells that are connected to each other but separated from other feeding zones by solid or porous cells.
I use a breadth-first search algorithm to identify isolated feeding zones. The algorithm is summarized in the following table.
| Step | Operation |
|---|---|
| 1 | Initialize a flag array. Set all cells to \(-1\), meaning unvisited. |
| 2 | Scan cells with flag \(-1\). If \(f_s \lt f_{cs}\), mark the cell as a new zone seed with flag \(n+1\). Otherwise set flag to 0. |
| 3 | Create a coordinate queue and push the seed. Pop a cell, examine its six neighbors. If a neighbor has \(f_s \lt f_{cs}\), assign the same flag and push it into the queue. Otherwise set its flag to 0. |
| 4 | Repeat until the queue is empty. Then return to step 2 to search for the next isolated feeding zone. |
| 5 | Stop when all cells have been visited. |
This algorithm has \(O(n)\) time complexity and \(O(n)\) space complexity in the worst case, where \(n\) is the number of cells. It provides the spatial distribution of isolated feeding zones, which is the basis for volume change calculation and cross-zone feeding.
For each cell \(i\) in an isolated feeding zone, I calculate the volume change during a time step. The volume change includes liquid shrinkage, graphite expansion, and austenite contraction:
$$\Delta V_{e,i} = \left(\alpha_l \Delta T_i (1 – f_{s,i}) – \alpha_G \Delta f_{G,i} + \alpha_A \Delta f_{A,i}\right)(1 – \phi_i) V_i$$
Here, \(\alpha_l\) is the liquid shrinkage coefficient, \(\alpha_G\) is the graphite expansion coefficient, \(\alpha_A\) is the austenite contraction coefficient, \(\Delta T_i\) is the temperature change, \(\Delta f_{G,i}\) is the graphite fraction change, \(\Delta f_{A,i}\) is the austenite fraction change, \(\phi_i\) is the porosity, and \(V_i\) is the cell volume. If \(\Delta V_{e,i}\) is positive, contraction dominates, and porosity increases. If \(\Delta V_{e,i}\) is negative, expansion dominates, and existing porosity can be reduced. If the expansion exceeds the local porosity, the excess liquid volume that can feed other cells is
$$V_{p,i} = \min\left(\left|\frac{\Delta V_{e,i}}{V_i}\right| – \phi_i, V_{l,i}\right)$$
The net volume change of an isolated zone is then
$$\Delta V_{iso} = \sum_i \phi_i V_i – \sum_j V_{p,j}$$
When \(\Delta V_{iso} \le 0\), the zone is an expansion zone. Its available expansion can feed other zones. When \(\Delta V_{iso} \gt 0\), the zone is a shrinkage zone. It has a net contraction that cannot be compensated internally and may become a shrinkage cavity or shrinkage porosity.
For shrinkage zones, I allocate the net contraction to cells according to a pressure field. I initialize the static pressure as
$$P_{init} = -\rho g \Delta z$$
where \(\Delta z\) is the vertical distance from a reference liquid surface. For open zones connected to air, I also account for atmospheric pressure. I then search from air cells into the alloy to mark cells connected to the atmosphere. For each open zone, the pressure is corrected by subtracting atmospheric pressure. Finally, the shrinkage volume is assigned to the cells with the lowest pressure first. If the shrinkage volume is small, it is distributed uniformly among the highest-priority cells. If it is larger than the liquid volume in those cells, those cells are fully converted to porosity and the remaining shrinkage is passed to the next priority level. If the shrinkage is intermediate, I sort the cells by liquid volume and allocate dynamically so that no cell receives more shrinkage than its liquid volume.
The more distinctive part of my model is the cross-zone feeding method among multiple isolated feeding zones in ductile iron castings. In ductile iron, graphite expansion can produce a significant internal pressure. An expansion zone can push liquid metal through narrow mushy channels or residual pores into a neighboring shrinkage zone. I treat the porous region between zones as a porous medium. The expansion pressure is calculated from the elastic response of the solid skeleton:
$$P_e = K \frac{-\Delta V_{iso}}{\Omega_{iso}}$$
where \(\Omega_{iso}\) is the total volume of the expansion zone and \(K\) is the bulk modulus:
$$K = \frac{E}{3(1 – 2\nu)}$$
Here, \(E\) is Young’s modulus and \(\nu\) is Poisson’s ratio. This pressure is the driving force for cross-zone feeding.
The flow resistance in the porous region is described by Darcy’s law:
$$\nabla P_s = -\frac{\mu}{K_p} v_l$$
where \(\mu\) is the dynamic viscosity of the liquid metal, \(K_p\) is the permeability, and \(v_l\) is the feeding velocity. I use the Kozeny-Carman relation for permeability:
$$K_p = K_0 \frac{\phi^3}{(1-\phi)^2}$$
with
$$K_0 = \frac{d^2}{180}$$
The secondary dendrite arm spacing \(d\) is calculated from the local solidification time:
$$d = \alpha t_f^n$$
where \(\alpha\) and \(n\) are constants. In my simulations, I use \(\alpha = 7.5 \times 10^{-6}\) and \(n = 0.33\). The porosity in the porous region is approximated by the sum of the local liquid fraction and porosity, \(\phi = V_l + \phi_i\). The feeding velocity is estimated from the solidification shrinkage rate:
$$v_l = -\frac{\rho_s – \rho_l}{\rho_l} \frac{H}{R} (V_l + \phi_i)$$
where \(H\) is the magnitude of the solidification rate and \(R\) is the magnitude of the temperature gradient. Combining these relations gives the pressure loss through a porous cell:
$$\Delta P_s = \frac{\rho_s – \rho_l}{\rho_l} \frac{H}{R} \frac{180}{(\alpha t_f^n)^2} \frac{(1 – V_l – \phi_i)^2 \mu}{(V_l + \phi_i)^2} dx$$
This equation shows that when the liquid or pore volume in a cell is very small, the pressure loss becomes extremely large. This means that a nearly solid cell blocks cross-zone feeding. When the liquid or pore volume is larger, the pressure loss is smaller, and liquid can pass more easily. This is consistent with the physical behavior of ductile iron castings, where residual liquid channels and small pores can still allow feeding under graphite expansion pressure.
To determine whether an expansion zone can feed a shrinkage zone, I calculate the minimum pressure loss along all possible paths. I represent each cell as a node and the pressure loss between neighboring cells as an edge weight. Then I apply Dijkstra’s algorithm to find the minimum pressure loss from the expansion zone to all reachable cells. I define the feeding potential as
$$G = P_e – \sum_i P_{s,i}$$
If \(G \gt 0\), the expansion pressure is large enough to overcome the pressure loss along the path, and cross-zone feeding can occur. The higher the feeding potential, the stronger the feeding ability. I also group shrinkage zones by feeding potential. If two shrinkage zones have feeding potentials that differ by less than a small threshold \(\gamma\), they are treated as the same priority level. In my implementation, \(\gamma\) is taken as the static pressure of one cell.
The cross-zone feeding algorithm is summarized below.
| Step | Operation |
|---|---|
| 1 | Assign each cell a flag indicating its isolated feeding zone. Non-zone cells have flag 0. |
| 2 | Initialize pressure loss. Set pressure loss to zero inside isolated zones. For other cells, calculate \(\Delta P_s\) from the Darcy model. |
| 3 | For each expansion zone, set the boundary cell feeding potential to \(P_e\). Insert these cells into a priority queue ordered by feeding potential. |
| 4 | Pop the cell with the highest feeding potential. If it belongs to a shrinkage zone and has not been recorded, record that zone and its total shrinkage volume. |
| 5 | For each neighboring cell, calculate \(G_{new} = G_{current} – \Delta P_s\). If \(G_{new} \gt 0\), insert the neighbor into the priority queue. |
| 6 | Continue until the priority queue is empty. This gives all shrinkage zones reachable from the expansion zone. |
| 7 | Group the reachable shrinkage zones by feeding potential. Allocate the expansion volume first to the highest-priority group. If the group has enough shrinkage volume, distribute the expansion proportionally. If not, fill the group and pass the remaining expansion volume to the next group. |
| 8 | Repeat for all expansion zones, then allocate the remaining shrinkage in shrinkage zones using the pressure field method. |
This algorithm allows me to simulate cross-zone feeding in ductile iron castings. An expansion zone in an early solidification stage can feed a shrinkage zone in a late solidification stage, even if the two zones are separated by a mushy or porous region. This is important because conventional models that assume isolated zones do not interact tend to overpredict shrinkage defects in ductile iron castings.
I implemented the entire model in C++ and integrated it with a casting simulation platform for mesh generation and visualization. The workflow consists of pre-processing, calculation, and post-processing. The pre-processor generates a uniform hexahedral mesh. The solver reads the mesh, assigns material properties, initializes temperatures, calculates the temperature field and solid fraction, identifies isolated feeding zones, computes volume changes, applies the cross-zone feeding algorithm, and writes porosity results. The post-processor displays temperature, liquid fraction, graphite fraction, and shrinkage porosity.
To verify the cross-zone feeding effect, I designed an irregular casting with two cylinders of different diameters connected by a narrow passage. The left cylinder is larger and cools more slowly, while the right cylinder is smaller and cools faster. During solidification, the narrow passage freezes first and separates the two cylinders into two isolated feeding zones. The right cylinder reaches the shrinkage stage earlier, while the left cylinder is still in the expansion stage because of graphite precipitation. With the cross-zone feeding model, liquid metal from the left cylinder can flow through the residual porous passage into the right cylinder, compensating part of its shrinkage. The simulation shows fewer defects in the right cylinder and more defects in the left cylinder compared with the case without cross-zone feeding. This result supports the physical consistency of the model for ductile iron castings.
I then applied the model to a drive-wheel-type ductile iron casting. The mesh size was 4 mm, and the total number of cells was 5,872,388, with 347,186 alloy cells. The alloy was QT600-3, and the initial casting temperature was 1480 °C. The mold was initially at 20 °C. The temperature field shows that the outer regions cool first, while the region connected to the runner and riser remains hotter and solidifies later. The cooling curves at several monitoring points show a clear temperature plateau, which is caused by the latent heat of the eutectic reaction. The graphite fraction distribution shows that graphite precipitates preferentially in regions with higher cooling rates. The liquid fraction distribution shows that the initially connected liquid region breaks into multiple isolated zones as solidification progresses. The predicted shrinkage porosity is dispersed, which is typical of ductile iron castings because of the nonuniform graphite expansion and the pasty solidification mode. The actual cut casting shows defect positions that agree with the predicted major defect regions.
I also applied the model to a bracket-type ductile iron casting with nonuniform wall thickness. The mesh size was 4 mm, and the total number of cells was 5,274,318, with 304,145 alloy cells. The alloy was QT700-2, and the initial casting temperature was 1370 °C. The temperature curves again show nonlinear behavior, with slow cooling during the eutectic stage. The predicted shrinkage defects are small and dispersed, which matches the pasty solidification behavior of ductile iron castings. After pouring and sectioning the actual casting, I compared the predicted shrinkage area with the actual shrinkage area. The predicted shrinkage area was 0.16 cm² in a section area of 17.04 cm², giving a predicted area fraction of 0.939%. The actual shrinkage area was 648 pixels in a section area of 64,917 pixels, giving an area fraction of 0.998%. The prediction accuracy, defined from the area fractions, was 94.09%. This result confirms that my model can quantitatively predict shrinkage defects in ductile iron castings.
The following table summarizes the simulation cases and validation results.
| Case | Alloy | Mesh size (mm) | Total cells | Alloy cells | Initial temperature (°C) | Validation result |
|---|---|---|---|---|---|---|
| Sample cup | QT400 | Uniform | – | – | 1300 | Temperature curve correlation 0.99357 |
| Irregular casting | QT600-3 | 3 | 482,400 | 73,477 | 1400 | Cross-zone feeding effect verified |
| Drive-wheel casting | QT600-3 | 4 | 5,872,388 | 347,186 | 1480 | Defect distribution agrees with cut casting |
| Bracket casting | QT700-2 | 4 | 5,274,318 | 304,145 | 1370 | Accuracy 94.09% |
My model has several innovations. First, it couples microstructural evolution with the macroscopic temperature field. The solid fraction is calculated from graphite nucleation, graphite growth, austenite shell growth, and a Johnson-Mehl relation. This produces a more realistic latent heat release and cooling curve for ductile iron castings. Second, it treats graphite expansion as a pressure source. The expansion pressure is calculated from the net volume change and the elastic bulk modulus of the solid skeleton. Third, it models the porous region as a flow resistance network using Darcy’s law and the Kozeny-Carman relation. Fourth, it uses Dijkstra’s algorithm to find the minimum pressure loss path between expansion zones and shrinkage zones. Fifth, it introduces a feeding potential and a priority grouping strategy to allocate expansion volume among multiple isolated feeding zones. This makes it possible to simulate cross-zone feeding in ductile iron castings, which is usually neglected in conventional shrinkage models.
There are also limitations. The thermophysical properties are treated as constants in the current implementation. In reality, density, thermal conductivity, and specific heat vary with temperature. Future work should include temperature-dependent properties obtained from measurements or inverse analysis. Gas precipitation is also not considered. In ductile iron castings, gas pores can form and interact with shrinkage porosity. This interaction should be included in future models. Finally, the graphite expansion is treated as an instantaneous source in each time step. The dynamic accumulation of expansion over time could be modeled more rigorously so that expansion generated at one time step can influence feeding at later time steps.
Despite these limitations, the model provides a practical and physically consistent method for predicting shrinkage defects in ductile iron castings. It can be used to evaluate riser design, gating system layout, and cooling conditions. It can also help reduce the number of trial pours and improve the yield of sound ductile iron castings. In my future work, I plan to extend the model to include variable thermophysical properties, gas porosity, and more detailed graphite expansion accumulation. I also plan to validate the model on a wider range of ductile iron castings with different wall thicknesses and alloy compositions.
In summary, I have developed a numerical simulation approach for shrinkage defect prediction in ductile iron castings. The approach combines microstructural modeling, temperature field calculation, isolated feeding zone search, volume change calculation, graphite expansion pressure, Darcy flow pressure loss, Dijkstra path search, and cross-zone feeding allocation. The temperature field was validated with measured cooling curves, and the shrinkage prediction was validated with actual pouring experiments. The prediction accuracy reached 94.09% in one bracket-type ductile iron casting. The results show that considering graphite expansion and cross-zone feeding is essential for accurate prediction of shrinkage defects in ductile iron castings. The model offers a useful tool for process optimization and quality control in ductile iron casting production.