A Comprehensive Framework for Predicting Shrinkage in Casting: Integrating SPH and FEM for High-Pressure Die Casting Simulation

The pursuit of high-integrity, near-net-shape components in manufacturing industries, particularly automotive and aerospace, has established high-pressure die casting (HPDC) as a critical process. Its ability to produce complex, thin-walled parts with excellent dimensional accuracy at high production rates is unparalleled. However, the intrinsic nature of the process—characterized by ultra-high velocity filling and rapid solidification under intense pressure—inevitably leads to the formation of internal defects. Among these, shrinkage in casting, manifesting as dispersed microporosity, is a predominant concern that severely compromises the mechanical performance and pressure tightness of the final component. Therefore, the ability to accurately predict the location and severity of shrinkage in casting during the product design phase is of paramount importance for optimizing process parameters, improving yield, and ensuring component reliability.

Numerical simulation has emerged as an indispensable tool for virtual prototyping in casting. Traditional simulation techniques often rely on grid-based methods like the Finite Volume Method (FVM) or Finite Element Method (FEM) for solving the governing equations of fluid flow, heat transfer, and stress. While powerful, these methods can struggle with the extreme distortions of the free surface, splashing, and jetting phenomena inherent to the HPDC filling stage. This is where meshless methods, particularly Smoothed Particle Hydrodynamics (SPH), offer a significant advantage. Conversely, for accurately calculating the thermal stresses and distortions during solidification under pressure, the FEM remains highly efficient and robust. This article presents a detailed, first-person perspective on the development and application of a coupled SPH-FEM numerical framework specifically designed for the comprehensive prediction of shrinkage in casting in high-pressure die casting processes.

The Challenge of Shrinkage in Pressure Die Casting

Shrinkage in casting defects are a direct consequence of the volumetric contraction of metal as it transitions from liquid to solid state. In pressure die casting, although the applied pressure acts to suppress pore formation, the rapid cooling often leads to isolated liquid pools that are cut off from the feeding source. When the localized solidification shrinkage cannot be compensated by incoming liquid metal, microporosity or shrinkage in casting forms. These defects typically reside in the last-to-freeze regions, such as thermal centers of thick sections, junctions, and areas adjacent to cooler dies where solidification is rapid and feeding paths are constricted. Predicting these areas requires an accurate simulation of the entire process chain: the turbulent flow during filling, the transient heat transfer, and the coupled thermomechanical effects during solidification.

The SPH Methodology for Filling and Thermal Analysis

Our approach employs the SPH method to simulate the initial stages. SPH is a Lagrangian, meshless particle method where the fluid domain is discretized into a set of interacting particles. Each particle carries field properties (mass, velocity, density, energy, temperature). The core principle is the kernel approximation, where the value of a function at a particle is approximated by a weighted average of the values at neighboring particles within a smoothing length \( h \). We adopt the cubic B-spline kernel function \( W \):

$$ W(R, h) = \alpha_d \begin{cases}
\frac{2}{3} – R^2 + \frac{1}{2}R^3 & \text{for } 0 \leq R < 1 \\
\frac{1}{6}(2-R)^3 & \text{for } 1 \leq R < 2 \\
0 & \text{for } R \geq 2
\end{cases} $$

where \( R = |\mathbf{r}_i – \mathbf{r}_j|/h \) is the normalized distance, and \( \alpha_d \) is a normalization constant dependent on dimensionality (\( 1/h, 15/(7\pi h^2), 3/(2\pi h^3) \) for 1D, 2D, and 3D, respectively). The discretized form of the Navier-Stokes equations for a weakly compressible fluid in SPH formalism are solved:

$$ \frac{d\rho_i}{dt} = \sum_j m_j \mathbf{v}_{ij} \cdot \nabla_i W_{ij} $$

$$ \frac{d\mathbf{v}_i}{dt} = -\sum_j m_j \left( \frac{P_i}{\rho_i^2} + \frac{P_j}{\rho_j^2} \right) \nabla_i W_{ij} + \sum_j m_j \frac{\mu_i + \mu_j}{\rho_i \rho_j} \mathbf{v}_{ij} \frac{1}{r_{ij}} \frac{\partial W_{ij}}{\partial r_{ij}} + \mathbf{F}^{(s)} $$

$$ \frac{de_i}{dt} = \frac{1}{2} \sum_j m_j \left( \frac{P_i}{\rho_i^2} + \frac{P_j}{\rho_j^2} \right) \mathbf{v}_{ij} \cdot \nabla_i W_{ij} $$

Here, \( \rho, \mathbf{v}, P, \mu, e, m \) represent density, velocity, pressure, dynamic viscosity, internal energy, and mass, respectively. \( \mathbf{F}^{(s)} \) encompasses body forces like gravity and the plunger pressure. The boundary condition is enforced using a dynamic virtual particle model, which creates layers of fixed particles to represent the die walls, interacting with fluid particles via the same equations.

The temperature field is simultaneously calculated by solving the energy equation. Crucially for predicting shrinkage in casting, we record the thermal history of each SPH particle. Upon complete filling, we have a highly accurate representation of the initial temperature distribution for the solidification analysis, capturing any cold shut regions or localized hot spots from the fluid flow.

Predicting Shrinkage in Casting with the Niyama Criterion within SPH

To directly identify regions prone to shrinkage in casting from the SPH thermal data, we integrate the well-established Niyama criterion into our particle-based framework. The Niyama criterion \( Ny \) is a local thermal parameter that has proven effective in predicting microporosity in alloys. It is defined as the ratio of the temperature gradient \( G \) to the square root of the cooling rate \( \sqrt{R} \) at the end of solidification:

$$ G = \sqrt{ \left( \frac{\partial T}{\partial x} \right)^2 + \left( \frac{\partial T}{\partial y} \right)^2 + \left( \frac{\partial T}{\partial z} \right)^2 } $$

$$ R = \frac{T_L – T_S}{t_L – t_S} $$

$$ Ny = \frac{G}{\sqrt{R}} $$

where \( T_L \) and \( T_S \) are the liquidus and solidus temperatures, and \( t_L \) and \( t_S \) are the times when the local temperature passes these points. Regions where \( Ny \) falls below a critical threshold \( Ny^* \) are predicted to contain shrinkage in casting. In our SPH implementation, for each particle representing the casting metal, we track its temperature over time during the solidification phase. The spatial derivatives for \( G \) and the time derivatives for \( R \) are calculated using the kernel approximation inherent to SPH. This provides a direct, particle-based map of the Niyama value across the entire casting, allowing us to flag particles below \( Ny^* \) as potential shrinkage in casting sites.

The FEM Methodology for Coupled Thermomechanical Solidification Analysis

While SPH excels at free-surface flow, FEM is more computationally efficient for solving the complex thermal stress and deformation during solidification under pressure. After the filling stage simulated by SPH, we map the particle-based temperature field onto a FEM mesh for the solidification stress analysis. The governing equation for the dynamic thermomechanical problem in FEM is:

$$ \mathbf{M}\ddot{\mathbf{u}}(t) + \mathbf{C}\dot{\mathbf{u}}(t) + \mathbf{K}\mathbf{u}(t) = \mathbf{Q}(t) $$

where \( \mathbf{M}, \mathbf{C}, \mathbf{K} \) are the global mass, damping, and stiffness matrices, \( \mathbf{Q} \) is the load vector, and \( \mathbf{u}, \dot{\mathbf{u}}, \ddot{\mathbf{u}} \) are displacement, velocity, and acceleration vectors. We utilize a lumped mass matrix for computational efficiency in explicit dynamics.

The key to accurate stress calculation is the constitutive model. We employ a thermal visco-elastoplastic model. The total strain increment \( \{\Delta \varepsilon\} \) is decomposed:

$$ \{\Delta \varepsilon\} = \{\Delta \varepsilon^{el}\} + \{\Delta \varepsilon^{th}\} + \{\Delta \varepsilon^{in}\} $$

The thermal strain \( \{\Delta \varepsilon^{th}\} \) is calculated from density changes. For the inelastic strain \( \{\Delta \varepsilon^{in}\} \) in the mushy and solid state, a rate-dependent model based on the hyperbolic sine law is used to capture the creep and plastic deformation of the semi-solid and solid material:

$$ \dot{\varepsilon} = A [\sinh(\alpha \sigma^*)]^n \exp\left(-\frac{Q}{R_g T}\right) $$

where \( \dot{\varepsilon} \) is the strain rate, \( \sigma^* \) is the flow stress, \( Q \) is the activation energy, \( R_g \) is the gas constant, and \( A, n, \alpha \) are material constants. This model allows us to compute the evolving stress field as the casting cools and contracts against the rigid die walls and under the applied plunger pressure.

A critical coupling effect considered is the pressure-dependence of the alloy’s solidification temperature, governed by the Clausius-Clapeyron relation:

$$ dT = \frac{T \Delta V}{\Delta H} dp $$

where \( \Delta H \) is the latent heat and \( \Delta V \) is the volume change. The FEM-calculated pressure field from the plunger action and solidification shrinkage therefore directly influences the local solidification temperature, refining the thermal field used for the Niyama-based shrinkage in casting prediction.

The SPH-FEM Coupling Strategy and Data Transfer

The seamless integration of the two solvers is vital. The coupling is sequential but bidirectional in terms of field data influence.

  1. Stage 1 (SPH): The filling process and initial solidification cooling are simulated entirely in SPH. Outputs include the final particle positions and their associated temperature history up to a transfer time.
  2. Stage 2 (Data Transfer): The temperature field from SPH particles is mapped onto the nodes of the FEM mesh. This is achieved using a radial-basis-function-based interpolation. For each FEM node, a search radius is defined, and all SPH particles within this radius contribute to the node’s initial temperature, weighted by their distance.
  3. Stage 3 (FEM): The solidification stress analysis is performed using FEM, with the mapped temperature as the initial condition and thermal load. The pressure from the plunger and thermal stresses are calculated. The modified solidification temperature due to pressure (from Eq. (16)) is fed back to refine the SPH-calculated thermal history for the final shrinkage in casting evaluation.
  4. Stage 4 (Final Shrinkage Prediction): The refined thermal history from the coupled simulation is used within the SPH post-processor to compute the final Niyama criterion map and identify shrinkage in casting regions.

This coupling leverages the strengths of both methods: SPH’s robustness for complex filling and FEM’s efficiency and accuracy for stress analysis.

Simulation Case Study: AM60B Magnesium Alloy Bracket

We applied our developed framework to simulate the HPDC process of a structural bracket made of AM60B magnesium alloy. The 3D model of the part and its gating system is complex, featuring thin walls and thicker sections. The key process parameters and material properties used in the simulation are summarized below.

Table 1: Key Process Parameters for HPDC Simulation
Parameter Value Unit
Plunger Pressure 11 MPa
Plunger Velocity 4 m/s
Pouring Temperature 680 °C
Die Initial Temperature 220 °C

Table 2: Thermophysical Properties of AM60B Alloy and H13 Die Steel
Material Property Value Unit
AM60B Liquidus Temperature 596 °C
Solidus Temperature 468 °C
Density (Liquid) 1,780 kg/m³
Specific Heat 1,050 J/(kg·K)
Latent Heat 370 kJ/kg
H13 Steel Density 7,850 kg/m³
Thermal Conductivity 28.6 W/(m·K)
Specific Heat 460 J/(kg·K)

Filling Pattern and Thermal Field Results

The SPH simulation captured the filling sequence in detail. The metal split into multiple streams upon entering the cavity, with the central stream impacting the far wall and spreading. The filling was completed in approximately 0.031 seconds. The simulated flow pattern showed excellent agreement with results from commercial FVM-based software (like ProCAST), particularly in predicting the flow fronts and potential cold shut locations. More importantly, the final temperature distribution at the end of filling, which serves as the starting point for the solidification stress analysis, was obtained with high fidelity, revealing localized hot spots in the thicker sections.

Stress Field and Refined Solidification Analysis

The FEM analysis of the solidification stage, initialized with the SPH temperature field, computed the evolution of stress. The results indicated that the areas experiencing the lowest pressure during solidification were not necessarily the geometric hot spots, but regions where the feeding path was long and restricted. The pressure from the plunger was effectively transmitted through the still-liquid channels, but diminished in isolated liquid pools. The Clausius-Clapeyron effect (Eq. (16)) was calculated, showing a slight but non-negligible increase in the local solidus temperature in high-pressure zones, slightly altering the local solidification time.

Prediction of Shrinkage in Casting

The final step was the application of the Niyama criterion to the refined thermal history. The critical Niyama value \( Ny^* \) for the AM60B alloy was determined through calibration. The resulting shrinkage in casting prediction map from our SPH-FEM coupled simulation was then compared against the prediction from a commercial software package using a traditional fully-FEM approach.

Table 3: Comparison of Shrinkage Prediction Results
Prediction Method Primary Shrinkage Locations Identified Remarks
SPH-FEM Coupled Framework Thermal centers of thick-walled bosses and hubs. Edges adjacent to thin, rapidly cooling die sections where feeding is restricted. Predictions are particle-based, offering a granular view. Incorporates stress/pressure effect on solidification temperature.
Commercial FVM/FEM Software Thermal centers of thick sections. Some isolated spots near thin walls. Predictions are element-based. Generally shows good agreement on major hot spots.

The comparison revealed a very high degree of concordance between the two methods regarding the major areas prone to shrinkage in casting. Both methods predicted the defect in the thermal centers of the thicker sections. Our coupled model provided additional insight by highlighting specific regions near thin walls. Cross-referencing with the FEM-calculated low-stress zones confirmed that these were areas where the feeding pressure was insufficient to overcome the friction and geometrical resistance of the long, tortuous, and partially solidified feeding paths, leading to inadequate compensation for solidification shrinkage and ultimately causing shrinkage in casting.

Discussion and Broader Implications

The successful development and application of this SPH-FEM coupled framework demonstrate its robustness and accuracy for predicting shrinkage in casting in high-pressure die casting. The key advantages of this integrated approach are:

  1. Superior Filling Physics: SPH naturally handles the free-surface fragmentation and splashing in HPDC, leading to a more accurate initial temperature distribution than traditional methods, which is critical for all subsequent calculations.
  2. Accurate Thermomechanics: FEM efficiently solves the complex, non-linear stress-strain problem during solidification under pressure, accounting for the visco-elastoplastic behavior of the semi-solid alloy.
  3. Enhanced Prediction Fidelity: The two-way coupling, where stress affects solidification temperature and vice-versa, refines the thermal history used for the Niyama criterion, leading to a more physically sound prediction of shrinkage in casting.
  4. Granular Defect Mapping: Representing shrinkage in casting at the particle level can be advantageous for subsequent analyses, such as quantifying local porosity for fatigue life prediction.

The framework is not without challenges. The data transfer between the disparate SPH particle and FEM mesh representations must be accurate and conservative. The computational cost of the SPH filling simulation, while advantageous for accuracy, is higher than that of a purely grid-based method for the same stage. Future work will focus on optimizing the coupling algorithm, exploring adaptive particle refinement in SPH, and integrating microstructure prediction models to link the shrinkage in casting prediction directly to mechanical property degradation.

In conclusion, the integration of Smoothed Particle Hydrodynamics for filling and initial solidification with the Finite Element Method for detailed thermomechanical stress analysis presents a powerful and comprehensive numerical tool for the foundry industry. It provides a deeper physical understanding of the formation mechanisms of shrinkage in casting defects in pressure die casting, moving beyond empirical guesswork to a science-based design and optimization paradigm. This capability is essential for manufacturing lighter, stronger, and more reliable cast components, particularly from lightweight alloys like magnesium, where defect tolerance is often low.

Scroll to Top