# 3D numerical simulation and experimental validation of resin-bonded sand gravity casting: Filling, cooling, and solidification with SPH and ProCAST approaches

Mohammad Zarbini Seydani, Abdelkader Krimi, Sofiane Khelladi, Marie

Bedel, Mohamed El Mansori

► To cite this version:

Mohammad Zarbini Seydani, Abdelkader Krimi, Sofiane Khelladi, Marie Bedel, Mohamed El Mansori. 3D numerical simulation and experimental validation of resin-bonded sand gravity casting: Filling, cooling, and solidification with SPH and ProCAST approaches. Thermal Science and Engineering Progress, 2024, 47, pp.102329. 10.1016/j.tsep.2023.102329. hal-05324621

HAL Id: hal-05324621

https://hal.science/hal-05324621v1

Submitted on 21 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Archive Ouverte - Open Repository

Science Arts & Métiers (SAM)

is an open access repository that collects the work of Arts et Métiers Institute of
Technology researchers and makes it freely available over the web where possible.

This is an author-deposited version published in: https://sam.ensam.eu
Handle ID: http://hdl.handle.net/10985/26215

To cite this version :

Mohammad ZARBINI SEYDANI, Abdelkader KRIMI, Sofiane KHELLADI, Marie BEDEL, Mohamed EL MANSORI - 3D numerical simulation and experimental validation of resin-bonded sand gravity casting: Filling, cooling, and solidification with SPH and ProCAST approaches - Thermal Science and Engineering Progress - Vol. 47, p.102329 - 2023
# 3D numerical simulation and experimental validation of resin-bonded sand gravity casting: Filling, cooling, and solidification with SPH and ProCAST approaches

Mohammad Zarbini Seydani<sup>a, *</sup>, Abdelkader Krimi<sup>b</sup>, Sofiane Khelladi<sup>a</sup>, Marie Bedel<sup>c</sup>,
Mohamed El Mansori<sup>c</sup>

## ARTICLE INFO

**Keywords:**

Gravity casting

SPH

Filling time and morphology

Solidification time

3D Printing mold

## ABSTRACT

This article provides a comprehensive investigation into the resin-bonded sand gravity casting process, with a focus on the filling, cooling, and solidification steps. The research combines numerical simulations and experimental validation in a three-dimensional (3D) configuration, utilizing a realistic filling system. The study employs three approaches—experimental tests, Smoothed Particle Hydrodynamics (SPH) simulations, and ProCAST simulations—to analyze the filling, cooling, and solidification steps. Two different molds were used in the experiments. The first mold has a transparent glass component in front of the plate, enabling observation and recording of the filling process, while the second mold, used solely for thermal analysis, did not incorporate any glass component. The SPH approach yields more accurate results for the filling time, liquid level height, and morphology when compared to ProCAST. The discrepancy in final filling time for the desired casting part between the experiment and SPH is 5.13 %, and the difference between the experiment and ProCAST is 15.38 %. Additionally, the cooling and solidification steps are investigated through an analysis of cooling curves. The numerical methods demonstrate slightly higher cooling rates and deviations in solidification times compared to the experimental data, mainly due to the thermal Neumann boundary condition. Furthermore, the average discrepancy in solidification time for five points of the intended casting component between the experiment and SPH is 8.60 %, whereas when compared with ProCAST, the discrepancy increases to 9.13 %.

## Introduction

Gravity casting is a metal pouring process where molten metal is transferred from a crucible into a mold solely through the force of gravity [1]. This technique comprises two primary steps: filling the mold with the liquid metal and solidification. The common objective of various industries is to improve product quality and minimize costs, and one effective approach to accomplish this is by leveraging the capabilities of numerical simulation [2–4]. Researchers have dedicated considerable efforts over the past few decades to investigating the correlation between the filling step and the quality of castings through a combination of numerical simulations and experimental validations [2,5]. Effectively forecasting and analyzing the behavior of molten metal during the filling, cooling, and solidification steps necessitates the

integration of computer modeling with practical investigations [6]. Extensive research has been conducted in the fields of thermal analysis, predictions of oxidation and free-surface flow to enhance our comprehension of the casting process [6,7].

These research studies were conducted using two different approaches: mesh-based and meshless methods [8]. The mesh-based method has certain challenges, including a high number of grid cells, lengthy computation time, and a lack of fluid history. On the other hand, the meshless technique is better suited for predicting free surfaces and oxidation defects compared to mesh-based methods. Smoothed Particle Hydrodynamics (SPH) is a Lagrangian method that is particularly suitable for simulating the filling process and solidification in gravity casting [9,10]. This method offers unique advantages, such as the ability to model the behavior of free surfaces and material interfaces, as well as

<footnote><sup>a</sup> Arts et Metiers Institute of Technology CNAM, LIFSE, HESAM University F-75013 Paris, France</footnote>

<footnote><sup>b</sup> Department of Civil, Geological and Mining Engineering, Polytechnique Montréal, Montréal, QC H3T 1J4, Canada</footnote>

<footnote><sup>c</sup> MSMP Laboratory, EA-7350, Arts et Metiers ParisTech, 2 Cours des Arts et Metiers, 13617 Aix en Provence, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: Mohammad.zarbini_seydani@ensam.eu (M.Z. Seydani).</footnote>
simulate complex physics involving solidification and multiphase flow [11–13]. While there have been numerous studies on casting modeling using SPH, only a few of them have successfully incorporated filling, cooling, solidification, and have validated their results using experimental data. During the past 20 years, much more information has become available on simulating casting processes using the SPH method especially high-pressure die casting [14].

Cleary et al. [15] and Ha and Cleary [16] conducted studies on the application of smoothed particle hydrodynamics (SPH) in two-dimensional high-pressure die casting (HPDC) and compared the results with water analog experiments. Ha and Cleary [16] performed isothermal SPH and volume-of-fluid (VOF) fluid simulations for HPDC, and the numerical results were compared with water flow tests for three different geometries. The SPH method demonstrated better prediction of separation from the corners, the timing of separation void disappearance, and the smoothness and shape of primary jet free surfaces [16]. For three-dimensional isothermal numerical simulation, Cleary et al. [17] developed a modeling approach for HPDC. The simulations revealed flow separation and significant fragmentation of free surfaces, leading to the formation of transient voids. To compare the SPH methodology with commercial software, Cleary et al. [18] contrasted the findings of MAGMAsoft simulation and SPH methodology using a water analog experiment for HPDC. The SPH model exhibited remarkably good agreement with fill pattern, speed prediction, and void trapping timings, unlike the MAGMAsoft simulation. Additionally, by incorporating an enthalpy formulation for solidifying metals, heat transfer and solidification predictions using SPH aligned well with experimental observations. He et al. [19] investigated two- and three-dimensional modeling of the filling stage in HPDC using smoothed particle hydrodynamics with artificial viscosity and the moving least squares method to address pressure oscillations. When comparing the results of SPH and the finite difference method (FDM) with experimental data, it was found that the SPH results were more consistent with the FDM. Cleary et al. [20] demonstrated the utility of SPH as a computational tool to predict die-filling behavior in complex thin-walled high-pressure die casting and compared the findings with water analog studies at similar Reynolds numbers. The rate of flow movement, the shapes and positions of voids during early filling, and the location and rate of the backfilling process matched perfectly. Hu et al. [11] employed smoothed particle hydrodynamics to develop a three-dimensional simulation of HPDC for a realistic die. The experimental results confirmed the excellent stability and reliability of the SPH model in evaluating the flow rate of molten metal during the cavity filling process. In their study, Tokunaga et al. [21] described a particle method simulation for the generation and flow of cold flakes in the high-pressure die casting process.

It is worth mentioning that smoothed particle hydrodynamics (SPH) has also been applied to other casting methods, including gravity and low-pressure casting. Ha et al. [9] conducted a study using SPH to model the filling process in two-dimensional gravity die-casting. They considered isothermal conditions and the thermophysical properties of water in two different shapes. The results showed good agreement with experiments in predicting the overall structure of the filling step. The SPH simulations successfully captured the behavior of free surface waves and provided detailed information about the flow. In the work of Cleary [8], an SPH approach was developed to predict shrinkage defects in two-dimensional simulations of the low-pressure casting process for an engine block. The model considered oxide production, feeding, solidification dynamics, and the distribution of residual pressure in the solidified metal. The ability to predict the amount and duration of liquid metal exposed to air was found to be crucial in estimating oxide concentration. Cao et al. [22] proposed a novel methodology for the foundry-filling process, utilizing composite solid boundary treatment and creating numerical models based on the SPH method. A traditional method was employed for the solidification model, incorporating temperature-dependent viscosity. However, the results obtained using this approach were not satisfactory. Lysenko et al. [23] utilized the SPH

approach to simulate the filling and crystallization stages of low-pressure casting. They found that the SPH method accurately simulated the filling process and could be applied in foundry manufacturing. The authors noted that under high pressure, supercooling reduced the solidification time of aluminum alloys. While several studies have focused on modeling the filling, solidification, and prediction of oxidation defects in casting processes, there has been a relatively limited amount of research that simultaneously investigates these processes in gravity and low-pressure casting. Zarbini et al. [24] utilized the SPH approach in numerical simulations and validated its accuracy through experimental tests to replicate the filling, cooling, and solidification processes of resin bonded sand gravity casting in a 2D scenario. By comparing the results obtained from both the experimental findings and simulations conducted with SPH and ProCAST, they effectively demonstrated the superior performance of SPH in accurately representing crucial aspects such as filling time, filling morphology, and volume flow rate.

This article presents a comprehensive study on resin bonded sand gravity casting, incorporating 3D numerical simulation and experimental validation. The goal of the research is to assess the accuracy and reliability of the Smoothed Particle Hydrodynamics (SPH) approach in comparison to ProCAST software. By conducting real gravity casting tests using transparent glass molds, both with and without glass, the researchers confirm the numerical simulation's ability to accurately replicate the filling and solidification processes.

## Gravity casting experiment

The experiments were carried out to validate the 3D SPH code and compare it with the ProCAST software while utilizing the actual filling system used for casting a plate. To conduct these tests, a 3D printed mold was employed, ensuring both high thermal stability and dimensional accuracy. The mold was filled with molten metal using gravity casting, which facilitated controlled and consistent flow. This combination of factors allowed for the precise replication of complex geometries and ensured reliable casting outcomes. The primary objective of these tests was to affirm the effectiveness of the 3D SPH code in simulating the gravity casting process with a resin bonded sand mold.

## Casting and molding design

A test case, as shown in Fig. 1, is designed to validate the accuracy and performance of the 3D SPH code in simulating gravity casting with a resin bonded sand mold.

The geometry shown includes all the essential elements required for a traditional casting experiment. These components include a pouring basin, a circular sprue, a runner, and an ingate. The plate is filled from the bottom through the ingate located at the center of the plate. The plate itself measures 200 mm in length, 100 mm in height, and 8 mm in thickness. It is important to note that the plate's thickness is intentionally much smaller than its length and height to ensure a uniform filling in the vertical direction and reduce the simulation's runtime. For the filling process, an ingate with a height of 5.6 mm and a length of 22.3 mm is utilized. The ingate is positioned 10 mm away from the end of the runner to limit the amount of sand that can be carried along with the molten metal. The runner, which connects the sprue to the ingate, has a square shape with dimensions of 11.1 mm and a length of 210 mm. The circular sprue employed has a diameter of 12.6 mm and a height of 100 mm. It is designed to have the same height as the plate to ensure complete filling. Additionally, a cylindrical pouring basin with a diameter of 80 mm and a height of 80 mm is considered. The top of the plate is intentionally left open to allow trapped air to escape without affecting the filling process.

The molds were produced using an ExOne S-Print Furan printer, which has a job box size measuring 800 × 500 × 400 mm³ [25]. The resin-bonded sand molds consist of a casting cavity, a base, a stopper,
<figure><img src="image_3.png" /><figcaption>Fig. 1. Schematic of geometry intended for gravity casting (mm).</figcaption></figure>

and a transparent glass. It is important to note that 3D printed resin-bonded sand molds need to exhibit sufficient casting resistance and allow for the release of gases. Casting simulation was employed using commercial software to determine the appropriate thickness of the sand mold surrounding the plate. Two different molds were printed for the experiments. The first mold included a transparent glass component in front of the plate, which allowed for observation and recording of the filling process of AlSi13 into a 3D printed sand mold bonded with resin. On the other hand, the second mold used for thermal analysis solely consisted of a 3D printed resin-bonded sand mold, without the incorporation of glass. The experimental mold with glass contains three components (positive side, negative side, and the glass), while the mold

without glass has just two parts (positive side and negative side). The positive side of the mold is larger in both cases and makes up the majority of the filling system, whereas the negative side is located in the lower part of the mold. Fig. 2 represents designing mold with glass and without glass with Catia and 3D printing resin - bonded sand mold with glass.

## In situ instrumentation

A resistance furnace of NPG.45 type, having a power output of 12 kW, voltage of 380 V, and current of 19A, is employed for the purpose of melting around 5 kg of AlSi13 (also referred to as EN AC-4400 by the

<figure><img src="image_4.png" /><figcaption>Fig. 2. Geometry of mold a) Designing mold with glass by CATIA b) Designing mold without glass by CATIA c) 3D printing resin - bonded sand mold with glass.</figcaption></figure>
Aluminum Association's US designations) at a temperature of 700 °C. The temperature of the molten metal when it enters the mold is not the same as that of the furnace, but it consistently stays higher than the temperature of the liquid metal. The metal has a superheat of around 75 °C, while the pouring temperature is approximately 655 °C. It's important to note that AlSi13 is slightly hyper eutectic due to its silicon content of 13 wt%, exceeding the eutectic point of the Al-Si alloy, which is 12.6 wt% [26]. AlSi13 is widely used in foundries because of its favorable properties, including good fluidity, castability, corrosion resistance, and crack resistance [27]. According to the ProCAST commercial software, the liquidus and solidus temperatures of AlSi13 are 572 °C and 570 °C, respectively, with a temperature difference of only two degrees. Consequently, the solidification process occurs at a relatively constant temperature, simplifying the analysis of solidification models.

An acquisition system and strategically positioned thermocouples are employed to measure the filling and solidification times at various stages, allowing for a comprehensive examination of these processes. To gain insights into the flow dynamics during the filling stage, an acquisition system equipped with electrical contacts is utilized to track the fluid movement in specific sections. Steel rods, each with a diameter of 1 mm, are precisely positioned at six distinct locations within the mold, as illustrated in Fig. 3. These rods serve as monitoring points to track the time it takes for the molten metal to fill the mold, starting from the top of the sprue and progressing towards the end of the plate. As the molten metal comes into contact with the rods, completing the electrical circuit, the resistance decreases, resulting in a corresponding drop in voltage. The output voltage data is collected at a fixed frequency of 100 Hz.

Furthermore, in order to examine the temperature changes within the molten metal, K-type thermocouples are strategically placed at six specific locations, as indicated in Fig. 3. Given that the maximum temperature recorded during the experiment is 655 °C, the thermocouples have a precision of ±4.91 °C (0.0075 × |T'|).

For the experiment using glass the acquisition system was put in nodes (N1, N2, N3, N4, N6, N7) to capture filling time (Fig. 4-4). The acquisition system and thermocouples were placed respectively in the nodes (N1, N2, N3, N4, N6, and N7) and (N1, N5, N6, N7, N8, and N9) for the experiment without glass in order to record the filling time and temperature (Fig. 4-4). N1 is positioned at a distance of 17.6 mm from the top of the sprue. N2 is 5.55 mm from the sprue bottom in the central portion of its thickness. To clarify further, a reference steel rod is positioned at the top of the sprue (N1), and the circuit is closed when the metal connects N1 to N2. This time represents "t = 0 s".

N3 is situated at the bottom of the runner, at the ingate entrance. N4 is in the middle of the ingate. By taking into account 20 mm, 50 mm, and 90 mm heights from the bottom of the plate, respectively, N5, N6, and N7 are installed at the middle of the thickness and length of the plate. By considering 50 mm and 25 mm from the plate's border, two additional points, defined N8 and N9, are embedded at the middle thickness of the plate and 50 mm height from the bottom of it. It should be noted that the filling process of casting is visible thanks to the use of 8 mm thick heat-resistant transparent silicone glass. The high-speed CR1000x3 camera can record the casting process filling. The 1280 x 512 pixel camera's field of view enables the recording of images while the casting is being filled. The camera is configured to take pictures at a 100 Hz frequency.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Location of thermocouples and tracking of molten metal (mm).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Voltage versus filling time for different points of experimental test.</figcaption></figure>

## Experimental results

### Filling time curves

Fig. 4 illustrates the voltage decline over time for two different experiments using glass (N1, N2, N3, N4, N6, and N7) and without glass (N1 to N7). Fig. 5 displays time the metal reaches the positions of the points N2 to N7 (bottom of sprue (N2), in the runner (N3), in the ingate (N4), and 20 mm, 50 mm and 90 m height of plate (N5, N6, and N7) respectively). Fig. 6 shows the height of the liquid at the center of the 200 mm-long plate.

### Point by point filling time analysis

The filling times at various points in the experiments were analyzed. The filling times at N3 and N4 were similar for both experiments, as shown in Figs. 4 and 5, with a relative error of less than 5 %. Nonetheless, there was a lack of data for the glass experiment at point N5, whereas the experiment without glass had recorded the filling time. Notably, at point N6, a notable difference of approximately 13 % was observed between the two experiments, with the glass experiment filling

at a faster rate. This difference can be attributed to the glass having a higher slip coefficient than the sand mold, which enables quicker filling. Based on the analysis of filling time at N7, it can be estimated that the experiments with and without glass would have filling times of 1.935 s and 1.965 s, respectively.

**Analysis of liquid level height:** Initially, between 15 mm and 40 mm, the plate fills up in a quick time interval of 0.2 s, indicating a high vertical filling velocity of approximately 0.125 m/s. The experiment results show that the liquid's height oscillates until it reaches 40 mm. When liquid reaches a height of 40 mm, it remains at this level height for 0.3 s (0.5 s to 0.8 s). This stability is due to the balance between liquid pressure and gravity, causing the metal to transition from vertical to horizontal flow. Between 0.7 s and 1.0 s, oscillations are observed between a height of 40 mm and 55 mm. From a height of 55 mm, the liquid exhibits a linear increase between 1.1 s and 1.93 s, with an average velocity of 0.042 m/s and minimal oscillations.

### Temperature time curves

In Fig. 7, you can observe temperature variations at six positions

<figure><img src="image_7.png" /><figcaption>Fig. 5. Points position versus filling time for different points of experimental test for whole part.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6. Liquid level height versus filling time of experimental test for plate.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Experimental cooling curves of AlSi13 alloy for points (N1, N5, N6, N7, N8 and N9).</figcaption></figure>

denoted as N1, N5, N6, N7, N8, and N9. These cooling curves exhibit three distinct zones that correspond to phase transitions in the alloy. The experimental cooling rates for Zones 1 and 3, as well as the solidification

<table><caption>**Table 1**
Experimental cooling rate and solidification time for (N5, N6, N7, N8, and N9).</caption><thead><tr><th>Position of thermocouples</th><th>Cooling rate (Liquid state) (T.s<sup>-1</sup>)</th><th>Solidification time (s)</th><th>Cooling rate (solid state) (T.s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>N5</td><td>16.27</td><td>83.36</td><td>0.91</td></tr><tr><td>N6</td><td>18.29</td><td>119.5</td><td>1.27</td></tr><tr><td>N7</td><td>17.38</td><td>109.17</td><td>1.07</td></tr><tr><td>N8</td><td>13.91</td><td>117.62</td><td>1.22</td></tr><tr><td>N9</td><td>13.60</td><td>94.63</td><td>0.95</td></tr></tbody></table>

time for Zone 2, are detailed in Table 1. The slopes in each of these zones were calculated using the highest and lowest temperature values along with their respective time points. Notably, the zero time point in Fig. 7 represents the initial closure of the electrical circuit at the base of the sprue, identified as N2. N1, located at the top of the sprue, is used to determine the casting's starting temperature. Two temperature analyses were conducted on the plate. The first analysis examines temperature changes vertically using N5, N6, and N7. These positions are situated at distances of 20 mm, 50 mm, and 90 mm from the bottom of the plate, respectively, at the middle point of the 4 mm thickness and 100 mm length. The second analysis focuses on temperature changes along the plate's length and considers N6, N8, and N9. These positions are located at the midpoint of the 50 mm height and the thickness, with distances of 100 mm, 50 mm, and 25 mm in the length direction, respectively.
Temperature analysis in height direction (N5, N6, and N7)

N6 exhibits a more significant temperature decrease in the liquid cooling area compared to N5 and N7, with drops of 12.41 % and 5.43 %, respectively. N5 and N6 initially have similar temperatures and enter the solidification zone simultaneously. However, N5 fills before N6, and N7 experiences a higher rate of temperature decline. Surprisingly, even though N6 is filled later, it undergoes a faster temperature change than N7 due to its lower initial temperature. In the solidification zone, N6 experiences the longest solidification time compared to N5 and N7, largely unaffected by height boundary conditions. N7, on the other hand, takes longer to solidify than N5 because it's filled to a height of over 100 mm, reducing the air's influence on N7. During the cooling of the solid state, N6 encounters a more rapid temperature decline than N5 and N7. N7 experiences a faster temperature drop than N5 due to N6's extended solidification, resulting in a quicker temperature drop as it enters the solid cooling zone. When N6 enters zone 3, the surrounding material has already solidified, leading to an accelerated temperature drop at N6.

Temperature analysis in length direction (N6, N8, and N9)

It is observed that the thermocouples at these points display the maximum temperature value at the same time. The maximum temperature decreases as one moves toward the mold's boundary, with N6 having a maximum temperature of 625 °C and N9 having a maximum temperature of 595 °C. Differences in cooling rates are observed in the first zone due to variations in maximum temperature values. N6 takes longer to solidify than N8 and N9, indicating a sharper decrease in solidification time approaching the mold wall. In the third zone, N6 undergoes a greater rate of temperature decline compared to N8 and N9. The cooling drop values show a small difference between N6 and N8, likely due to their similar solidification times. It is noteworthy that N6 experiences a higher temperature drop than N9, despite being closer to the mold wall, by approximately 25.2 %.

Numerical methodologies

This article discusses a research study that focuses on modeling the filling and solidification processes in 3D sand gravity casting. The study utilizes two different numerical approaches: the meshless method, specifically Smoothed Particle Hydrodynamics (SPH), and the mesh-based method, using commercial software, for simulation. The accuracy and reliability of these models are subsequently confirmed through experimental validation.

SPH method

Discretization of SPH approach

SPH is a technique that employs particles to simulate various phenomena, including fluid flow, heat transfer, and intricate physics like multiphase flow and solidification [13,14]. Unlike conventional Eulerian methods [14,18], SPH operates on a Lagrangian continuum principle, allowing each particle to automatically track complex flows. To achieve this, SPH employs interpolation kernels [14,18] to smoothen particle density, pressure, and other individual properties based on neighboring particles. The SPH approach expresses the interpolated value of a field function A at any given location r [14,18].

$$A(r_i) = \sum_{j=1}^{N} \frac{m_j}{\rho_j} A(r_j) W_{ij} \quad (1)$$

This equation signifies the summation of all j particles within the support domain of the i particle, which extends up to a radius of 2h from r. $m_j$, $\rho_j$, and $W_{ij} = W(r_i - r, h)$ correspond to the mass, density and the kernel function, respectively [14]. The constant value of the smoothing length, denoted as “h,” is determined based on the initial inter-particle

distance $\delta x_0 (h = 2\delta x_0)$ [12,14]. The initial volume of the particles, represented by $V_0 = \delta x_0^d$, is dependent on the number of space dimensions, denoted as d. The Kernel function adheres to certain properties, including smoothness, compact support, normalization, and regularity, as stated in references [11,12,14].

$$W(r, h)_{3D} = \alpha_{3D} \begin{cases} \left(3 - \frac{r}{h}\right)^5 - 6^* \left(2 - \frac{r}{h}\right)^5 + 15 \left(1 - \frac{r}{h}\right)^5 & 0 \le \frac{r}{h} < 1 \\ \left(3 - \frac{r}{h}\right)^5 - 6^* \left(2 - \frac{r}{h}\right)^5 & 1 \le \frac{r}{h} < 2 \\ \left(3 - \frac{r}{h}\right)^5 & 2 \le \frac{r}{h} < 3 \\ 0 & \frac{r}{h} > 3 \end{cases} \quad (2)$$

Where $\alpha_{3D}$ is equal to $\frac{3}{359\pi h^3}$. The fluid properties are determined by employing smoothing gradients. The gradient of the function A is then given by [14]:

$$\nabla A(r_i) = \sum_{j=1}^{N} \frac{m_j}{\rho_j} A(r_j) \nabla W_{ij} \quad (3)$$

Several enhancements have been implemented to enhance the SPH method, aiming to simulate free surface flows with reduced numerical noise and tensile instability [28,29]. In this study, the $\delta$-SPH technique [30,31] is employed to enhance accuracy and minimize numerical noise when simulating casting process in the Navier-Stokes and energy equations. The continuity, momentum, energy, and equation of state were modeled using weakly compressible and Newtonian assumptions. Further details can be found in [30,31].

Discrete form of governing equations

To accurately model the physics involved in the sand gravity casting process during filling and solidification, it is essential to address and solve the continuity, momentum, and energy equations [32,33]. The following equations represent the discrete form of continuity, Navier-Stokes, energy, and displacement equations for the $\delta$-SPH approach [24]:

$$\left\{ \begin{aligned} \frac{D\rho_i}{Dt} &= -\rho_i \sum_j (V_j - V_i) \cdot \nabla_i W_{ij} \cdot V_j' + h c_0 \sum_j \delta D_{ij} \cdot \nabla_i W_{ij} \cdot V_j' \\ \frac{DV_i}{Dt} &= -\frac{1}{\rho_i} \sum_j (\rho_i + p_j) \nabla_i W_{ij} \cdot V_j' + h c_0 \frac{\rho_0}{\rho_i} \sum_j \alpha \pi_{ij} \nabla_i W_{ij} \cdot V_j' + g_i + \frac{A(\phi_s) V_i}{\rho_i} \\ \tilde{C}_i \frac{DT_i}{dt} &= \nabla \cdot (\langle k_i \rangle \nabla T_i) \\ \tilde{C} &= \langle C \rangle + \rho_i L \frac{d\phi}{dT} \\ \frac{D \vec{r}}{Dt} &= \vec{V} \end{aligned} \right. \quad (4)$$

In Cartesian coordinates, the substantial derivative is denoted by $\frac{D(r)}{Dt} \cdot \nabla_i$ represents the gradient operation, commonly known as the coordinates of the particle. Furthermore, This symbol $\langle \cdot \rangle$ defines as follows:

$$\langle \cdot \rangle = (1 - \phi)(\cdot)_s + \phi(\cdot)_i \quad (5)$$

$\rho$, $p$, $V$, $r$, $T$, $k$, $V_j'$, $g$, $L$ and $C$ represent density, pressure, velocity vector, position vector, temperature, thermal conductivity, volume, gravitational acceleration vector, latent heat and volumetric heat capacity, respectively.

In the continuity equation, the diffusive term is represented by the symbol “$D_{ij}$.”
$$D_{ij} = 2[(\rho_j - \rho_i) - \frac{1}{2}(\langle\nabla\rho\rangle_j^L + \langle\nabla\rho\rangle_i^L) \cdot (r_j - r_i)] \frac{(r_j - r_i)}{|r_j - r_i|^2} \quad (6)$$

The symbol $\langle\nabla\rho\rangle_i^L$ indicates the renormalized density gradient that has been defined as [29]:

$$\langle\nabla\rho\rangle_i^L = \sum_j (\rho_j - \rho_i) L_i \nabla_i W_{ij} V_j^{-1} \quad (7)$$

$$L_i = \left[ \sum_j (r_j - r_i) \otimes \nabla_i W_{ij} V_j^{-1} \right]^{-1} \quad (8)$$

$\delta$ is an important parameter that plays a significant role in determining the extent of density diffusion [30]. The permissible range for this dimensionless parameter is restricted, and it is typically set to 0.1 in $\delta$-SPH.

In the momentum equation, $\pi_{ij}$ serves as an artificial viscosity to enhance stability [34], while $\alpha$ is a constant, dimensionless parameter ranging from 0.01 to 0.05 [30].

$$\pi_{ij} = \frac{(V_j - V_i) \cdot (r_j - r_i)}{(r_j - r_i)^2} \quad (9)$$

The variable $A(\phi_s)$ is a linear function employed to introduce solidification into the momentum equation, whereas $\phi_s$ represents the solid fraction [35].

$$A(\phi) = -C_A(1 - \phi) \quad (10)$$

In this study, $C_A$ represents a constant parameter with a value of $10^6$. This linear function has the advantage of being computationally inexpensive. The effect of the porosity function on the momentum equation is elucidated in three consecutive steps as follows [24].

* In the entire liquid region (when $\phi = 1$), the function $A$ becomes zero and therefore does not affect the momentum equation.

* The porosity function gradually changes between the liquid and solid regions, particularly in the mushy/regularization region (when $0 < \phi < 1$).

* In the complete solid region (i.e., when $\phi = 0$), the function $A$ assumes an extremely large value. Consequently, it overwhelms all other terms in the momentum equation when considering velocity changes over time. As a result, the predicted velocity is set to zero.

Energy equation is written based on the temperature and volumetric heat capacity which is totally explained in previous research. In energy equation, smooth fraction function address both isothermal and non-isothermal phase transition problems [36,37]. The smooth fraction function can be represented by a sine function with respect to temperature [24].

$$\phi(T) = \begin{cases} 0 & T < T_r - \epsilon \quad (\text{solid}) \\ \frac{1}{2} \left\{ 1 + \sin\left(\pi\left(\frac{T - T_r}{2\epsilon}\right)\right) \right\} & T_r - \epsilon \le T \le T_r + \epsilon \quad (\text{mushy}) \\ 1 & T > T_r + \epsilon \quad (\text{liquid}) \end{cases} \quad (11)$$

In Eq. (11), $\epsilon$ denotes a smoothing parameter that is used with a Heaviside step function within a region of $2\epsilon$ when considering isothermal phase transition. The reference temperature, denoted as $T_r$, varies depending on the type of process, whether it is isothermal or non-isothermal. In the case of an isothermal phase change process with a pure material, $T_r$ is equal to the melting temperature ($T_r = T_m$). On the other hand, in non-isothermal phase change phenomena, $T_r$ is determined as the average of the liquidus and solidus temperatures ($T_r = \frac{T_r + T_f}{2}$). For more details, refer to [24,30,31].

The WCSPH technique [38], known as weakly compressible smoothed particle hydrodynamics, is utilized in this research. The pressure calculation in this study is based on the equation of state presented as follows:

$$P = P_r \left( \frac{\rho}{\rho_0} - 1 \right) + P_b \quad (12)$$

Where $\rho_0$ represents the reference density, $P_r$ indicates the reference pressure and $P_b$ shows the background pressure. The reference pressure is determined by the reference density and the artificial sound speed $c_0$ in the linear constitutive equation of state given by Eq. (12).

$$P_r = \rho_0 c_0^2 \quad (13)$$

After considering the stability criteria, employing the actual speed of sound as a point of reference leads to an exceedingly small time step. Consequently, the determination of $c_0$ is accomplished in the following manner.

$$c_0^2 \ge \max \left\{ \frac{U_0^2}{\delta\rho}, \frac{\|g\| L_0}{\delta\rho}, \frac{\mu U_0}{\rho_0 L_0 \delta\rho} \right\} \quad (14)$$

The velocity reference, length reference, and dynamic viscosity are denoted as $U_0$, $L_0$, and $\mu$, respectively. The dimensionless density variation, represented by $\delta\rho$, is set at 1 % according to reference [39]. In order to simulate free surface flow, the background pressure is typically set to zero ($P_b = 0$). However, in confined flow situations, a positive value is assigned to this parameter to prevent tensile instability. It is worth noting that the maximum density variation, ($\delta\rho$), is approximately 1 %.

### Boundary condition

#### Fluid boundary condition:

In this study, we employed an adaptation of the generalized wall boundary condition method proposed by [39,40], which accounts for multiphase fluid flow. To put this method into practice, three layers of dummy particles were positioned along the wall interface in a perpendicular manner, as shown in Fig. 8. This method allows for the establishment of either free-slip or no-slip wall boundary conditions. When calculating the viscous forces exerted by the fluid, the free-slip boundary condition is employed by ignoring the viscous interaction between the fluid particle and the neighboring dummy particles.

As suggested by [39,40], the pressure in the dummy-wall particle is

<figure><img src="image_10.png" /><figcaption>Fig. 8. Visual Representation Depicting Different Parameters Involved in the Generalized Wall Boundary Condition [39].</figcaption></figure>
determined based on the pressure in the adjacent fluid particles j.

$$p_w = \frac{\sum_{i}^{n_f} p_j W_{ij} + g \sum_{j}^{n_f} \rho_j r_{wj} W_{wj}}{\sum_{j}^{n_f} W_{ij}} \quad (15)$$

### Thermal boundary condition

The previous work has introduced the non-homogeneous thermal neumann boundary condition $\nabla T \cdot n = q$ [24]. In order to express the thermal Neumann boundary conditions on the material particles (including liquid, mushy/regularization, and solid regions), the following equation has been developed, utilizing the appropriate wall particles.

$$T_{wall_i} = \frac{q \sum_{j}^{n_f} r_{ij} \frac{\partial W}{\partial r_j} V_j + d \sum_{j}^{n_f} T_j V_j \nabla W \cdot n_i}{d \sum_{j}^{n_f} V \nabla W \cdot n_i} \quad (16)$$

Here $q$ shows a heat flux ($q = 0$ for the adiabatic conditions) and $n_i$ represents the normal vector to the wall.

### ProCAST simulation

ProCAST is a commercially available software that utilizes the finite element method to simulate the complete casting process, including both the filling and solidification stages. In order to explore the characteristics of gravity casting in a mold made of resin-bonded sand, a thorough three-dimensional thermal-hydraulic simulation was carried out using the ProCAST software [41,42]. The simulation focused on the AlSi13 alloy and treated the liquid alloy as an incompressible Newtonian fluid [43,44]. By solving the 3D continuity and Navier-Stokes equations, along with incorporating heat conservation principles, the filling and solidification processes of gravity casting were modeled [45,41,42].

$$\text{Continuity equation : } \nabla \cdot \vec{V} = 0 \quad (13)$$

$$\text{NavierStokes (Momentumequation) : } \rho \frac{D\vec{V}}{Dt} = \rho \vec{g} - \nabla p + \nabla \cdot (\mu \nabla \vec{V}) \quad (17)$$

$$\text{Energy equation : } \rho c_p \frac{DT}{Dt} = \nabla \cdot (k \nabla T) + \rho L \frac{\partial \phi_s}{\partial t} \quad (18)$$

$$\text{Solidification Path : } \phi_s = \left( \frac{T - T_l}{T_l - T_s} \right) \quad (19)$$

In the ProCAST simulation, a linear solidification path is assumed. The time steps for filling and solidification have minimum values of 0.01 s and maximum values of 0.1 s and 1 s, respectively. The thermal

<caption>**Table 2**</caption>

Thermal properties of the resin-bonded sand and the AlSi13 alloy in ProCAST database.

<table><thead><tr><th>Temperature (°C)</th><th>Density (kg·m<sup>-3</sup>)</th><th>Specific heat (J·kg<sup>-1</sup>·K<sup>-1</sup>)</th><th>Conductivity (W·m<sup>-1</sup>·K<sup>-1</sup>)</th><th>Viscosity (kg·m<sup>-1</sup>·s<sup>-1</sup>)</th></tr></thead><tbody><tr><td colspan="5"><strong>Resin-bonded sand properties</strong></td></tr><tr><td>20</td><td></td><td>670</td><td>0.71</td><td></td></tr><tr><td>100</td><td></td><td>800</td><td>0.68</td><td></td></tr><tr><td>200</td><td></td><td>920</td><td>0.64</td><td></td></tr><tr><td>300</td><td></td><td>883</td><td>0.60</td><td></td></tr><tr><td>400</td><td>1590</td><td>1006</td><td>0.56</td><td></td></tr><tr><td>500</td><td></td><td>1006</td><td>0.53</td><td></td></tr><tr><td>600</td><td></td><td>1006</td><td>0.50</td><td></td></tr><tr><td>700</td><td></td><td>1006</td><td>0.60</td><td></td></tr><tr><td>900</td><td></td><td>1006</td><td>0.73</td><td></td></tr><tr><td colspan="5"><strong>AlSi13 alloy properties</strong></td></tr><tr><td>300</td><td>2600</td><td>1020</td><td>151.6</td><td>100</td></tr><tr><td>450</td><td>2578</td><td>1130</td><td>147</td><td>100</td></tr><tr><td>550</td><td>2563</td><td>1200</td><td>144</td><td>100</td></tr><tr><td>570</td><td>2560</td><td>1210</td><td>143.5</td><td>100</td></tr><tr><td>572</td><td>2478</td><td>1210</td><td>80</td><td>0.001</td></tr><tr><td>750</td><td>2431</td><td>1210</td><td>80</td><td>0.001</td></tr></tbody></table>

properties of the resin-bonded sand and AlSi13 alloy can be found in Table 2 of the ProCAST database. The latent heat value of AlSi13 is 522 kJ/kg, while the liquidus and solidus temperatures are recorded as 572 °C and 570 °C, respectively. The density of the resin-bonded sand is assumed to be constant. When the molten metal enters the mushy area below 570 degrees and at a velocity of less than 100 kg.m/s, the viscosity of AlSi13 increases, resulting in the prioritization of the viscosity term and halting fluid movement to simulate the solidification process.

### Initial and boundary conditions

#### Initial condition

**Hydraulic:** It's crucial to highlight that, in fluid dynamics, the fluid begins in a state of rest, as was observed in the conducted experiment. This occurs when the appropriate initial and boundary conditions are introduced. The flow of the fluid is then naturally initiated as a consequence of gravity.

**Thermal:** The thermocouple (N1) is situated at the topmost section of the sprue. The initial condition for numerical simulation is set based on the approximate maximum temperature recorded at this location, which is 655 °C. In contrast, the initial temperature of the mold wall is 20 °C, matching the temperature of the surrounding environment.

#### Boundary condition

**Hydraulic:** The WSHEAR (wall shear) technique is employed in the ProCAST simulation to consider the presence of a velocity boundary layer along the mold wall. The parameter WALLF is assigned a value of 0.99, indicating a higher degree of wall slips. Additionally, the SPH method defines Perfect slip as the boundary condition between the wall and the fluid.

**Thermal:** To model the cooling and solidification process, both WCSPH and ProCAST simulations utilize the Neumann boundary condition (Conduction). This condition helps estimate the heat extraction at the interface between the metal and the mold. The heat flux leaving the mold surface is assumed to be constant across three different levels, as depicted in Table 3. These levels correspond to the cooling of molten metal above the liquidus temperature, the solidification process between the liquidus and solidus temperatures, and the cooling of solids at temperatures below the solidus temperature. The values for these levels are determined through one-dimensional calculations based on experimental data. The following method demonstrates how to calculate the value of "q" using one-dimensional equations. $T_0$, $T_L$, $T_S$, and $T_{end}$ shows initial, liquidus, solidus, and temperature at 150 s. Furthermore, the average time intervals between specific temperatures in each location (N5 to N9) were employed in order to calculate these values.

*Due to the temperature gradient between the mold and the molten material in the first period, when the mold is still cold and the metal is liquid, $q$ is equal $195000W/m^2$, the heat extraction is high. Then, in the second and third stages, the value decreases to $44000W/m^2$ and $11300W/m^2$, respectively. Fig. 9 shows $q$ values of Neumann boundary condition versus temperature and time for SPH and ProCAST method, respectively.

<caption>**Table 3**</caption>

One-dimensional equation to obtain Neumann boundary conditions in three different zones.

<table><tr><td>Zone (1)</td><td>Q = -K <sup>∂T</sup> / ∂x = <sup>ρeC<sub>p</sub>(T<sub>0</sub> - T<sub>L</sub>)</sup> / <sub>t<sub>exp</sub></sub> ≈ 195000W/m<sup>2</sup></td></tr><tr><td>Zone (2)</td><td>Q = -K <sup>∂T</sup> / ∂x = <sup>ρeL + ρeC<sub>p</sub>(T<sub>L</sub> - T<sub>S</sub>)</sup> / <sub>t<sub>(exp)</sub></sub> ≈ 44000W/m<sup>2</sup></td></tr><tr><td>Zone (3)</td><td>Q = -K <sup>∂T</sup> / ∂x = <sup>ρeC<sub>p</sub>(T<sub>S</sub> - T<sub>end</sub>)</sup> / <sub>t<sub>(exp)</sub></sub> ≈ 11300W/m<sup>2</sup></td></tr></table>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. q values of Neumann boundary condition versus temperature and time for SPH and ProCAST method.</figcaption></figure>

### Mesh and convergence

CATIA V5 software is used to create the CAD models of the part, the gating system, and the sand mold before being imported into ProCAST.

For the part and mold, the model is meshed using linear tetrahedral elements with a mesh size of 2 mm. The mesh size is established via a number of mesh refinement trials in order to confirm that ProCAST and SPH simulations are mesh-independent. The trend of mesh dependence

<figure><img src="image_13.png" /><figcaption>Fig. 10. Mesh independency of filling step for ProCAST Simulation.</figcaption></figure>
for filling is shown in Figs. 10 and 11. According to ProCAST simulation results, a larger filling velocity value is expected when mesh is refined, and this value converges as the cell size goes down. Fig. 10 shows that when the cells are 2 mm in size or smaller, variation is less than 5 % in total filling time with reducing cell size. As a result, mesh is generated in a 2 mm size. This leads to the selection of 596,578 cells, which correspond to 2 mm in size. With the SPH method, the trend is reversed. When there are more particles, the velocity is decreased. The results of the SPH technique are considered to become independent of the number of particles once the number of particles reaches 184,839 corresponding to 1.8 mm between each particle. Fig. 10 shows that when the distance between each particle is between 1.7 mm and 1.8 mm, variation is less than 3 % in total filling time.

Fig. 12 represents the final mesh used for ProCAST simulation. Fig. 13 also illustrates the SPH method’s mesh generation. The geometry must be created using CATIA by designing the desired part’s outer shell in order to generate mesh for the SPH method. Then, four particles should be positioned in the outer direction perpendicular to the shell, with an embedded thickness proportional to the space between each particle. The appropriate geometry file is given to a program that can generate a triangular mesh, such as ProCAST, QuikCAST, Abaqus, etc. The file is sent to MATLAB software because the generated triangular mesh may not have exactly 2 mm between each point. The distance between the points is organized in MATLAB software, and an embedded geometry shell is created with a thickness equal to four times the distance between each particle. In the context of the SPH approach, it is worth mentioning that the time step has a value of 2.5e-4 s.

# Results and discussion

## Study of filling step

Three different approaches are employed to analyze the filling results in this study. Firstly, the filling time of six points in the mold cavity is considered, and experimental data, SPH simulation, and ProCAST simulation are presented in Fig. 14 to compare the results with and without glass. Secondly, the filling time is correlated with the liquid level height in the plate, focusing on the middle section of the plate, as shown in Fig. 15. Lastly, the filling morphology is examined using a 2D analysis, with a camera placed in front of the plate. Fig. 16 illustrates the filling morphology using three techniques (experiment with glass, SPH, and ProCAST), as well as a comparison of the volume flow rate over time for each method in the 2D flow.

## Filling time analysis

Table 4 displays the time it takes to fill points N3, N4, N5, N6, and N7 in seconds, based on the average values from two experiments conducted with both glass and without glass, using the SPH method and ProCAST simulation.

The results obtained from the acquisition system shown in Fig. 14 demonstrate that it takes approximately 0.15 s to go from N2 (bottom of the sprue) to N3 (beginning of the ingate) based on experimental data. The SPH approach carefully estimates the filling time as 0.15 s, while the simulation using ProCAST slightly overestimates it as 0.16 s. The experimental filling time between N3 (beginning of the ingate) and N4 (middle of the ingate) is only 0.012 s, with corresponding values of 0.045 s for SPH and 0.06 s for ProCAST. According to the experimental results, it takes 0.088 s for the filling time between N4 and N5, which is 20 mm above the bottom of the plate. Both SPH (0.045 s) and ProCAST (0.04 s) show a negligible underestimation of the filling times between N4 and N5. The filling times for N6 (50 mm height from the bottom of the plate) and N7 (90 mm height from the top of the plate) reveal that the numerical approaches SPH and ProCAST overestimate and underestimate the filling time compared to the experimental results, respectively. The outcome indicates that the discrepancies between the experiment and SPH for points N6 and N7 are approximately 2.94 % and 5.13 %, respectively, while the differences between the experiment and ProCAST for these two points are 11.11 % and 15.38 %, respectively. It is worth noting that the flow initially becomes turbulent as the liquid enters the plate through the small ingate, characterized by high velocity and changing cross-section. The flow remains turbulent until the free surface becomes flat. After 0.6 s, when the liquid becomes flat, the Reynolds number drops below 1000, indicating a transition to the laminar phase of flow.

## Liquid level height

Fig. 15 represents the level of liquid height in the center of the plate versus filling time for the experiment and numerical (SPH and ProCAST) results.

In Fig. 15, the initial position of each curve varies for the three methods because of different filling times. The starting point is determined by examining the initial image on the plate, through both numerical analysis and experimentation. Initially, the liquid fills the plate and rises to a height controlled by pressure, subsequently spreading longitudinally after passing through the ingate located at the center of the plate. ProCAST simulation accurately predicts the filling time, evident from the comparison of liquid level up to 0.5 s, including appropriately predicted oscillations. However, in ProCAST simulation, the middle height of the liquid oscillates instead of remaining constant.

<figure><img src="image_14.png" /><figcaption>Fig. 11. Mesh independency of filling step for SPH Simulation.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 12. Schematic of mesh of ProCAST simulation.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. SPH geometry schematic with 1.8 mm distance between particles.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Filling time of several points in the mold cavity for Experiment, SPH method and ProCAST simulation.</figcaption></figure>
Liquid level height - filling time

<figure><img src="image_18.png" /><figcaption>Fig. 15. Liquid level height in the center of plate versus filling time for Experiment, ProCAST simulation and SPH method.</figcaption></figure>

In contrast, the SPH method maintains a nearly constant liquid level of around 27 mm during the 0.35 s to 0.65 s period of longitudinal filling. Afterwards, liquid level variations become oscillatory in both the experiment and ProCAST simulation, but follow a nearly linear trend in the SPH method. Notably, the SPH and experiment results are almost identical after 0.75 s. Comparing the SPH simulation to the experiment, the error in predicting liquid level height relative to filling time is below 10 % between 0.75 s and 2.05 s.

Filling morphology

Fig. 16 shows the filling morphology of the sand gravity casting process using the SPH method, experiment, and ProCAST simulation. The figure displays the filling process at 0.1-second intervals from 0 s to 1.7 s. The zero reference time in this figure is when the liquid enters the 200 mm wide plate through the ingate. The calculations focus on the filling height, up to 90 mm from the plate's bottom. Image analysis enables the determination of the liquid volume entering the plate by calculating the area of liquid and assuming uniform filling across the thickness. Fig. 17 presents the volume flow rate versus time.

At the initial stage (0.0 s), the plate fills up from the bottom through the right side ingate. The ProCAST simulation showed a higher volume of liquid entering the plate at 0.1 s compared to SPH and the experiment. However, the distribution of the liquid differed, with the experiment showing vertical filling and SPH predominantly filling along the plate's length. At 0.2 s, the volume flow rate was higher in both SPH and Pro-CAST compared to the experimental data. The filling morphology also differed, with the numerical methods reaching the wall in the longitu-dinal direction, while the experiment did not. ProCAST distributed the liquid more evenly along the wall than SPH. By 0.3 s, the experimental results showed the fluid reaching the wall on the right side of the plate, while the numerical methods differed in their extent of wall filling. The discrepancy between experimental and SPH results in volume flow rate is significant. Volume flow rate for Experiment, SPH, and ProCAST are $0.81 \times 10^{-4}m^3/s$, $1.06 \times 10^{-4}m^3/s$, and $1.30 \times 10^{-4}m^3/s$, respectively. At 0.5 s, the SPH method successfully predicted the free surface of the liquid with two dips. The volume flow rate differences between the experiment and numerical methods were reduced. At 0.6 s, the SPH approach showed a nearly flat fluid surface, while ProCAST and the experiment displayed oscillations. The volume flow rate discrepancies between SPH and the experiment reversed, and ProCAST still over-estimated the flow rate. By 0.7 s, all methods showed a nearly horizontal liquid surface, and the experimental results revealed a rise beyond the

plate's bottom. The volume flow rates decreased in the experiment and SPH, while increasing in ProCAST. From 0.8 s to 1.7 s, the liquid morphology changed only in the vertical direction, and the volume flow rates gradually declined. The SPH approach properly modeled the effects of metal height changes. It should be noted that the minimum volume flow rate of ProCAST simulation (at time 1.7 s) is higher than the maximum volume flow rate of SPH and Experiment (at time 0.8 s). Hence, it is obvious that ProCAST simulation grossly overestimates volume flow rate during casting process filling. The results demonstrate that the two numerical approaches differ in filling time, morphology of filling, and volume flow rate during the filling stage of gravity casting process. The variations between the results of these two numerical simulations can be attributed to a number of factors like physical models and numerical methods.

In physical model, ProCAST is mesh-based method (Finite Element Method), while SPH is mesh-less method. Mesh-based techniques are typically Eulerian: they simulate fluid flows using a fixed grid. SPH, on the other hand, employs a Lagrangian approach. This simply shows that particles move in accordance with the flow rather than being fixed in space. SPH is a great method for examining casting process filling because there is no mesh distortion, which enables efficient simulation of the free surface. In general, the SPH allows for a better understanding of the fine details of the flow, such as the filling time and morphology.

When it comes to the numerical method, weakly compressible flow is taken into account in the SPH method, whereas incompressible flow is taken into account in the ProCAST method. The equation of state is used by the SPH method, whereas ProCAST does not. This means that each method has a different way of calculating pressure. Because of the various pressure and density equations, the characteristics of the flow change and are diverse in two ways.

Using the two distinct numerical approaches, the pressure gradient and diffusion term of the Navier-Stokes equation are derived and discretized in different ways. Due to this, the filling results of gravity cast parts for both numerical methods could differ. Examining mathematical operators such as gradient, divergence and Laplace and discretization methods used in FEM and SPH show that the discretization of Navier stocks and continuity equation is different.

Study of cooling and solidification step

Fig. 18 displays the cooling curves for points N5, N6, N7, N8, and N9, while Table 5 provides the associated cooling rates and solidification
<figure><img src="image_19.png" /><figcaption>Fig. 16. Morphology of filling process versus time – a) SPH (Left) b) Experiment (middle) c) ProCAST (Right).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 16. (continued).</figcaption></figure>

times for three distinct zones. It's important to highlight that the qualitative comparison between the experimental and numerical results is influenced by the absence of experimental data regarding the thermal properties of both the sand and the alloy.

N5: Comparing the cooling rates in Zone One of Fig. 18 and Table 5 at N5, located in the middle of the plate with a height of 20 mm, reveals that the SPH and ProCAST methods exhibit higher experimental cooling rates by 14.19 % and 17.70 % respectively. The examination of

experimental and numerical data for solidification time (Zone 2) indicates that the SPH approach overestimates it, whereas ProCAST underestimates it compared to the experimental values. Similarly, numerical simulations tend to overestimate the cooling rate in Zone 3, as observed in Zone 1.

Fig. 18 demonstrates the influence of thermal boundary conditions as a function of time on the ProCAST simulation, resulting in a step-like behavior in temperature evolution within the third zone. At 100 s, the
<figure><img src="image_21.png" /><figcaption>Fig. 16. (continued).</figcaption></figure>

metal remains in the second zone at the center of the plate, while it reaches the third zone near the mold boundary. This moment coincides with a sudden change in the value of *q* in the boundary conditions, leading to a gradual decrease in the rate of heat exit at the mold borders. The metal in the center acts as a heat source, transferring heat to the metal near the boundaries. Consequently, heat transfer from the center to the boundaries attempts to balance out, resulting in an almost zero rate of temperature drop.

N6: The SPH method provides predictions that are approximately 2 % closer to the experimental results than ProCAST. However, in zone 1 at N6, which is positioned at the center of the plate in terms of width, thickness, and height, the experimental method still exhibits a higher cooling rate. The solidification time values are comparable for all three methods at N6. ProCAST underestimates by roughly 1.5 %, while SPH overestimates by about 2 %. Notably, N6, located at the plate's center, has the longest solidification time. In the third zone, similar to the first
<figure><img src="image_22.png" /><figcaption>Fig. 16. (continued).</figcaption></figure>

<caption>**Table 4**</caption>

Filling time (s) of N3 to N7 in the mold cavity for Experiment, SPH method and ProCAST simulation.

<table><thead><tr><th>Methods/Points</th><th>N3</th><th>N4</th><th>N5</th><th>N6</th><th>N7</th></tr></thead><tbody><tr><td>Experiment</td><td>0.15</td><td>0.16</td><td>0.25</td><td>0.99</td><td>1.95</td></tr><tr><td>SPH</td><td>0.15</td><td>0.2</td><td>0.29</td><td>1.02</td><td>2.05</td></tr><tr><td>ProCAST</td><td>0.16</td><td>0.22</td><td>0.26</td><td>0.88</td><td>1.65</td></tr></tbody></table>

zone, SPH forecasts a slower cooling rate, whereas ProCAST shows a significantly higher value, which is abnormal.

N7: Like other points, Nè shows a higher cooling rate in the experimental results compared to the numerical ones. Positioned at 90 mm from the plate's bottom and at the midpoint of its width and thickness, N7's initial temperature differs from the SPH and ProCAST methods. The SPH method predicts solidification time more reliably with a 7 % deviation from the experimental results, while ProCAST's solidification time is 22.5 % lower. Both methods indicate a lower cooling rate in the third zone, but SPH outperforms ProCAST in prediction accuracy. It's worth noting that, similar to N5, a change in boundary condition at 100 s causes a temperature drop, followed by a rise from 100 s to 140 s. This

highlights the limitation of using Neumann's temperature boundary conditions over time.

N8: The experimental cooling rate at N8, located at a height of 50 mm in the middle of the thickness and 75 mm from the transverse wall, exceeds SPH and ProCAST by 6.6 % and 14.8 % respectively. In the second zone, the predicted solidification time values closely resemble N6 due to its proximity to the plate's center, with both numerical methods slightly overestimating the value. Notably, in the ProCAST simulation, the solidification time at N8 is unexpectedly higher than N6, which is abnormal. In the third zone, the SPH method predicts a temperature drop approximately 15.1 % lower than the experimental results, while ProCAST forecasts a cooling rate that is 12.9 % higher than the experimental findings.

N9: The SPH method shows an overestimation of 1.5 % at N9, which is located 50 mm away from the mold in the width direction and in the center of height and thickness. On the other hand, the ProCAST method underestimates by 12.5 %. Both approaches significantly overestimate the solidification time. In zone 3, both numerical methods indicate a higher cooling value. The SPH method predicts 7.4 % more than the experimental results, while the ProCAST method predicts 16.8 % more. The notable difference between the numerical and experimental

<figure><img src="image_23.png" /><figcaption>Fig. 17. Volume flow rate versus filling time into the plate for EXP, SPH and ProCAST method.</figcaption></figure>
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 18. Temperature evolution of point N5, N6, N7, N8, and N9.</figcaption></figure>

**Table 5**

<table><caption>Comparison of temperature drop rates (C°/s) in liquid (zone 1) and solid (zone 3) and solidification time (s) (zone 2) for N5, N6, N7, N8 and N9 using experimental data, SPH and ProCAST.</caption><thead><tr><th>Methods/Points</th><th colspan="3">N5</th><th colspan="3">N6</th><th colspan="3">N7</th><th colspan="3">N8</th><th colspan="3">N9</th></tr><tr><th>Region</th><th>1</th><th>2</th><th>3</th><th>1</th><th>2</th><th>3</th><th>1</th><th>2</th><th>3</th><th>1</th><th>2</th><th>3</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td>Experiment</td><td>16.27</td><td>97.16</td><td>1.04</td><td>18.29</td><td>119.5</td><td>1.27</td><td>17.38</td><td>109.17</td><td>1.07</td><td>15.32</td><td>117.62</td><td>1.22</td><td>13.60</td><td>94.63</td><td>0.95</td></tr><tr><td>SPH</td><td>13.96</td><td>105.4</td><td>0.95</td><td>14.41</td><td>121.95</td><td>1.03</td><td>14.72</td><td>116.80</td><td>0.97</td><td>14.30</td><td>118.81</td><td>1.06</td><td>13.80</td><td>117.81</td><td>1.02</td></tr><tr><td>ProCAST</td><td>13.39</td><td>86.8</td><td>1.02</td><td>14.12</td><td>117.7</td><td>1.33</td><td>14.34</td><td>84.49</td><td>0.94</td><td>13.04</td><td>119.96</td><td>1.40</td><td>12.58</td><td>114.58</td><td>1.11</td></tr></tbody></table>

outcomes for the solidification time in the second zone might be the reason for this overestimation in zone 3.

In the initial zone, it was observed that experimental tests exhibited a higher rate of temperature reduction compared to numerical simulations in most cases. This dissimilarity can be attributed to the presence of air in the experimental tests, which was not accounted for in the numerical simulations. In the second zone, solidification time in experiments consistently proved to be shorter than the predictions of the SPH method at all points. This discrepancy arises from the calculation of the “q” value in the Neumann boundary condition for each zone in the numerical simulations, which utilizes an analytical equation based on specific assumptions. Specifically, these assumptions involve complete contact between the liquid and the mold, as well as infinite thermal conductivity

of the molten alloy. Consequently, these assumptions lead to an underestimation of the “q” value in the boundary conditions, resulting in longer solidification times in the numerical simulations.

Applying Neumann boundary conditions as a function of time leads to different slopes in the cooling zone for liquid and solid regions with varying thicknesses. There are also issues when applying Neumann’s boundary conditions in terms of temperature. The “q” value experiences three-step decreases from zone 1 to zone 3, as illustrated in Fig. 9 and Table 3. When the temperature of the metal near the boundaries falls below the solidus temperature, the metal transitions from zone 2 to zone 3. Consequently, the thermal boundary conditions in the walls undergo a sudden change, causing a sharp drop in the value of “q.” At this point, the metal in the center of the plate remains in zone 2, leading to an
inaccurate prediction of solidification time in the middle of the plate.

By employing the thermal Neumann boundary condition as a function of temperature, the temperature evolution in the SPH method becomes more uniform across all points of the plate compared to the experimental results. As shown in Table 5, the solidification time values for all points in the SPH method are closer to each other than they are to the experimental results. In terms of cooling rate, this conclusion holds true for both the first and third zones. A comparison of results obtained by applying Neumann boundary conditions as a function of time and temperature indicates that applying boundary conditions based on temperature yield more accurate results.

## Conclusion

This article focused on the 3D numerical simulation and experimental validation of filling, cooling, and solidification for resin bonded sand gravity casting. The objective was to assess the accuracy and reliability of the Smoothed Particle Hydrodynamics (SPH) approach in modeling the filling, cooling, and solidification steps of sand gravity casting in a closed system, compared to commercial software like ProCAST.

The study conducted a real 3D gravity sand casting test using two molds—one with transparent glass and the other without—to facilitate the observation of the filling stage. By comparing the experimental observations with the numerical simulations using image analysis, it was demonstrated that the SPH approach could accurately replicate the filling process. The shape of filling, filling time, and volume flow rate were analyzed and found to be in good agreement between the numerical and experimental data. Additionally, the study examined the application of Neumann’s boundary conditions to temperature in the absence of glass. The SPH results showed that this approach improved the prediction of temperature evolution across all three zones, indicating a more homogeneous modeling of the process.

Through these real experimental gravity sand casting tests, the 3D SPH code developed in this thesis was validated. The results confirmed the reliability of the SPH approach in simulating the filling, cooling, and solidification processes of rapid gravity casting. Moreover, a comparison with ProCAST highlighted the superior performance of the SPH method in accurately modeling filling time, filling morphology, and volume flow rate.

In conclusion, the utilization of glass in real gravity casting experiments enabled the direct observation and validation of the filling process, supporting the results obtained through numerical simulations using the 3D SPH code and commercial software such as ProCAST. This research contributes to the evaluation and comparison of different modeling approaches, affirming the SPH method as an effective tool for simulating the complex dynamics of filling and solidification in resin bonded sand gravity casting.

## CRediT authorship contribution statement

**Mohammad Zarbini Seydani:** Investigation, Data curation, Methodology, Writing – original draft. **Abdelkader Krimi:** Methodology, Conceptualization. **Sofiane Khelladi:** Methodology, Conceptualization. **Marie Bedel:** Investigation, Methodology, Conceptualization. **Mohamed El Mansori:** Supervision, Conceptualization, Methodology.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

The data that has been used is confidential.

## Acknowledgements

The authors acknowledge the contribution of colleagues. Thanks are due to J. Bourgeois and J. Négre of the Arts et Metiers ParisTech for their technical support. Funding: Institut Carnot ARTS.

## References

[1] R.A. Harding, M. Wickins, H. Wang, G. Djambazov, K.A. Pericleous, Development of a turbulence-free casting technique for titanium aluminides, Intermetallics 19 (2011) 805-813.

[2] M.V. Abhilash Viswanath, S. Manu, U.T.S.P. Savithric, Numerical simulation and experimental validation of free surface flows during low pressure casting process, J. Mater. Process. Technol. 244 (2017) 320-330.

[3] P.D. Ingle, B.E. Narkhede, A literature survey of methods to study and analyze the gating system design for its effect on casting quality, Mater. Today.: Proc. 5 (2018) 5421-5429.

[4] H.U. Rasheed, S.I. Zeeshan, B. Ali, Q. Shah, R. Ali, Implementation of shooting technique for Buongiorno nanofluid model driven by a continuous permeable surface, Heat Transfer 52 (4) (2023).

[5] B. Ali, S. Jubair, Rheological properties of Darcy-Forchheimer hybrid nanofluid flow with thermal emission and heat source over a curved slippery surface, Pramana 127 (2023).

[6] K. Ellingsen, T. Coudert, M. M'Hamdi, SPH based modelling of oxide and oxide film formation in gravity die castings, IOP Conf. Ser.: Mater. Sci. Eng. 84 (2015).

[7] M. Prakasha, P. Cleary, J. Grandfield, Modelling of metal flow and oxidation during furnace emptying using smoothed particle hydrodynamics, J. Mater. Process. Technol. 209 (2009) 3396-3407.

[8] P.W. Cleary, Extension of SPH to predict feeding, freezing and defect creation in low pressure die casting, App. Math. Model. 34 (2010) 3189-3201.

[9] Joseph Ha, Paul Cleary, Vladimir Alguine and Thang Nguyen, "Simulation of die filling in gravity die casting using sph and magmasoft," Second International Conference on CFD in the Minerals and Process Industries, p. 423-428, 1999.

[10] T. Ye, D. Pan, C. Huang, M. Liu, Smoothed particle hydrodynamics (SPH) for complex fluid flows: recent developments in methodology and applications, Phys. Fluids 31 (1) (2018).

[11] M.Y. Hu, J.J. Cai, N. Li, H.L. Yu, Y. Zhang, B. Sun, W.L. Sun, Flow modeling in high-pressure die-casting processes using sph model, Int. J. Met. 12 (2018) 97-105.

[12] X. Niu, J. Zhao, B. Wang, Application of smooth particle hydrodynamics (SPH) method in gravity casting shrinkage cavity prediction, Comput. Particle Mech. 6 (2019) 803-810.

[13] P. Cleary, J. Ha, M. Sawley, Modelling industrial fluid flow application using SPH, in: H.B. Muhlhaus (Ed.), Bifurcation and Localization Theory in Geomechanics, Swets and Zeitlinger, Lisse, 2001.

[14] M.B. Liu, G.R. Liu, Smoothed Particle Hydrodynamics (SPH): an Overview and recent developments, Arch. Comput. Methods Eng. 17 (2010) 25-76.

[15] S.R. Idelsohn, E. Onate, To mesh or not to mesh? That is the question, Comput. Methods Appl. Mech. Eng. 195 (2006) 4681-4696.

[16] J. Ha, P.W. Cleary, Comparison of SPH simulations of high pressure die casting with the experiments and VOF simulations of Schmid and Klein, Int. J. Cast Met. Res. (2000).

[17] P.W. Cleary, J. Ha, Three-dimensional smoothed particle hydrodynamics simulation of high pressure die casting of light metal components, J. Light Metlas (2002) 169-183.

[18] P.W. Cleary, J. Ha, M. Prakash, T. Nguyen, 3D SPH flow predictions and validation for high pressure die casting of automotive components, App. Math. Model. 30 (2006) 1406-1427.

[19] H. Yi, Z. Zhao-Yao, C. Wen-Jiong, C. Wei-Ping, Simulation of mould filling process using smoothed particle hydrodynamics, Trans. Nonferrous Met. Soc. Chin. 21 (12) (2011) 2684-2692.

[20] P.W. Cleary, G. Savage, J. Ha, Mahesh. Prakash., Flow analysis and validation of numerical modelling for a thin walled high pressure die casting using SPH, Comput. Particle Mech. 1 (2014) 229-243.

[21] H. Tokunaga, Y. Motoyama, T. Okane, Particle method simulation for formation and flow of cold flakes in high-pressure die casting, Int. J. Met. 13 (2019) 897-904.

[22] W.-J. Cao, Z.-Y. Zhou, F.-M. Jiang, Smoothed particle hydrodynamics modeling and simulation of foundry filling process, Trans. Nonferrous Met. Soc. China 25 (2015) 2321-2330.

[23] T. Lysenko, Y. Morozov, K. Kreitser, E. Kozishkurt, Using the sph method for modeling the crystallization process of aluminum alloys., Eng. Sci. 3 (55) (2020) pp.

[24] M. Zarbini Seydani, A. Krimi, M. Bedel, S. Khelladi, M. El Mansori, A 2D Filling and Solidification Benchmark Test: Validating Smoothed Particle Hydrodynamics (SPH) Simulations for Sand Gravity Casting, Int. J. Adv. Manuf. Technol. (2023).

[25] N. Coniglio, T. Sivarupan, M. El Mansori, Investigation of process parameter effect on anisotropic properties of 3D printed sand molds, Int. J. Adv. Manuf. Technol. 94 (2018) 2175-2185.

[26] M. Zamani, Al-Si Cast Alloys -Microstructure and Mechanical Properties at Ambient and Elevated Temperature, Jönköping University, School of Engineering, Dissertation, 2017.

[27] M.B. Saada, M. El Mansori, Assessment of the effect of 3D printed sand mold thickness on solidification process of AlSi13 casting alloy, Int. J. Adv. Manuf. Technol. 114 (2021) 1753-1766.
[28] P.N. Sun, A. Colagrossi, S. Marrone, A.M. Zhanga, The δplus-SPH model: Simple procedures for a further improvement of the SPH scheme, Comput. Methods Appl. Mech. Engrg. 315 (2017) 25–49.

[29] M. Antuono, A. Colagrossi, S. Marrone, Numerical diffusive terms in weakly-compressible SPH schemes, Comput. Phys. Comm. 183 (2012) 2570–2580.

[30] D.D. Meringolo, S. Marrone, A. Colagrossi, Y. Liua, A dynamic δ-SPH model: How to get rid of diffusive parameter tuning, Comput. Fluids 179 (2019) 334–355.

[31] D.D. Meringolo, Y. Liu, X.-Y. Wang, A. Colagrossi, Energy balance during generation, propagation and absorption of gravity waves through the δ-LES-SPH model, Coast. Eng. 140 (2018) 355–370.

[32] B. Ali, N.K. Mishra, K.R.S. Jubair, Z. Mahmood, S.M. Eldin, Mixed convective flow of hybrid nanofluid over a heated stretching disk with zero-mass flux using the modified Buongiorno model, Alex. Eng. J. 72 (2023) 83–96.

[33] H.A. Othman, B. Ali, S. Jubair, M.Y. Almusawa, S.M. Aldin, Numerical simulation of the nanofluid flow consists of gyrotactic microorganism and subject to activation energy across an inclined stretching, Sci. Rep. 13 (2023).

[34] J.J. Monaghan, R.A. Gingold, Shock Simulation by the Particle Method SPH, J. Comput. Phys. 52 (1983) 374–389.

[35] A. Brent, V. Voller, K. Reid, Enthalpy-porosity technique for modeling convection-diffusion phase change: application to the melting of a pure metal, Numer. Heat Transfer, Part A Appl. 13 (3) (1988) 297–318.

[36] G.H. Meyer, Multidimensional stefan problems, SIAM J. Numer. Anal. 10 (3) (1973) 522–538.

[37] P.W. Egolf, H. Manz, Theory and modeling of phase change materials with and without mushy regions, Int. J. Heat Mass Transf. 37 (18) (1994) 2917–2924.

[38] J.J. Monaghan, simulating free surface flows with sph, J. Comput. Phys. 110 (2) (1994) 399–406.

[39] A. Krimi, M. Rezoug, S. Khelladi, X. Nogueira, M. Deligant, L. Ramírez, Smoothed Particle Hydrodynamics: A consistent model for interfacial multiphase fluid flow simulations, J. Comput. Phys. 358 (2018) 53–87.

[40] S. Adami, X. Hu, N. Adams, A generalized wall boundary condition for smoothed particle hydrodynamics, J. Comput. Phys. 231 (21) (2012) 7057–7075.

[41] P. Lan, J. Zhang, Study on the mechanical behaviors of grey iron mould by simulation and, Mater. Des. 53 (2014) 822–829.

[42] L.u. Su-Ling, F.-R. Xiao, S.-J. Zhang, Y.-W. Mao, B.o. Liao, Simulation study on the centrifugal casting wet-type cylinder liner based on ProCAST, Appl. Therm. Eng. 73 (2014) 512–521.

[43] Hua. Wu, Haiming. Shi, Haifeng. Liu, Zhenjia.Xia, "Numerical Simulation of Flow Field and Temperature Field on Aluminium Alloy Engine Cylinder in Casting Process," Mater. Sci. Forum, Vols. 704-705, pp. 50-57, 2012.

[44] L.P. Ning, Li,, Low pressure casting technology and forming process analysis of metal mold based on ProCAST FEA procedure, Mech. Adv. Mater. Struct. (2020).

[45] B. Ali, M. Mashael, AlBaidani 2 Sidra Jubair 3,4 Abdul Hamid Ganie5, "Computational framework of hydrodynamic stagnation point flow of nanomaterials with natural convection configured by a heated stretching sheet," ZAMM Z. Angew. Math. Mech. (2023).