Trans. Nonferrous Met. Soc. China 25(2015) 2321-2330

Smoothed particle hydrodynamics modeling and simulation of foundry filling process

Wen-jiong CAO^1,2, Zhao-yao ZHOU², Fang-ming JIANG¹

1. CAS Key Laboratory of Renewable Energy, Institute of Energy Conversion, Chinese Academy of Sciences, Guangzhou 510640, China;

2. School of Mechanical and Automotive Engineering, South China University of Technology, Guangzhou 510640, China

Received 15 July 2014; accepted 30 September 2014

Abstract: A numerical model of foundry filling process was established based on the smoothed particle hydrodynamics (SPH) method. To mimic the constraints that the solid mold prescribes on the filling fluid, a composite treatment to the solid boundaries is elaborately designed. On solid boundary surfaces, boundary particles were set, which exert Lennard-Jones force on approaching fluid particles; inside the solid mold, ghost particles were arranged to complete the compact domain of near-boundary fluid particles. Water analog experiments were conducted in parallel with the model simulations. Very good agreement between experimental and simulation results demonstrates the success of model development.

Key words: smoothed particle hydrodynamics; foundry filling process; composite boundary treatment; water analog experiment

1 Introduction

As the first stage of foundry, the filling process largely influences the quality of the final casting products. Many defects such as gas porosity, flow mark, oxide inclusion, cold shut, and misrun may form during the filling process if they are not properly handled [1]. Modeling the filling process and predicting the relevant behaviors via numerical simulations can reveal the underlying fundamentals and design/optimize the process, hence are useful to shorten the trial period, to reduce the cost, and to enhance the quality of products. To date, most numerical models of the casting filling process are Eulerian, such as the finite difference method (FDM) [2], the finite volume method (FVM) [3], and the finite element method (FEM) [4]. These methods require meshes or grids to discrete the mold and the filling regions; the simulation results are very sensitive to the mesh quality and the artificial (numerical) diffusion may not be avoidable due to the extremely complex geometry of interest. Furthermore, the Eulerian methods cannot directly capture the free surface of fluid and the interface between materials [5], so the volume of fraction (VOF) function or level-set method merged with the Eulerian methods for modeling the filling process may further aggravate the artificial diffusion [6].

In recent years, the smoothed particle hydrodynamics (SPH) method has received increasing attention due to its mesh-free and Lagrangian properties [7]. Unlike traditional Eulerian methods, SPH method takes the continuum fluid as a set of Lagrangian particles. Each particle carries mass, momentum, enthalpy, and other physical quantities of fluids and evolves with time by following conservative laws. The free surface of fluid and the interface between materials can be naturally described by the distribution of particles. SPH has been proved to be more precise in modeling free surface flow with effective avoidance of the artificial diffusion [8].

In the past decades, many research works [9-11] focused on precision, stability, efficiency and other issues of SPH method and expended considerable efforts to develop more advanced SPH algorithms for free surface flow simulations. MONAGHAN [12] applied SPH to simulate incompressible flow with free surface, and a stiff equation of state linking the density and pressure fields was introduced in his work. The SPH method has also been successful in applications to the simulation of two-phase flow [13], seepage flow in porous media [14], and heat conduction [15]. CLEARY et al [16-18] established multi-dimensional SPH models for the high pressure die casting, and compared the SPH simulation results with the MAGMAsoft simulation results and the water analogue experimental results.

SPH simulations to the foundry filling process need to properly tackle the involved solid boundaries. Theoretically, the Lennard-Jones force [7] can be applied to mimic arbitrary solid boundary due to its flexibility. However, its accuracy is not high, as in the regions near the boundary, the compact domain of SPH particle is cut off. An other approach is the ghost particle method [7], which projects some SPH particles into the solid regions to complete the compact domain of SPH particles near the boundary. Although the ghost particle approach has been successfully used to model dam-break process [19], it is ascertained that this approach is only competent of tackling simple boundaries such as straight lines and flat planes. Recently, ADAMI et al [20] improved the ghost particle boundary treatment based on local force balance between wall and fluid particles. Their method shows advantages in dealing with complex geometries, but may not be able to completely avoid particle penetration when it is used to mimic solid boundaries involved in foundry filling process.

HE et al [21] reported an SPH model for mould filling process, in which the artificial viscosity and moving least square method were implemented to handle the numerical instability of fluid pressure. As a continuation, the present work aims at developing an effective SPH model of the foundry filling process. To this goal, a composite approach was proposed to tackle the involved boundaries of foundry mold and modify the link-list algorithm to improve the calculation efficiency. We also design and carry out water analog experiments to validate the effectiveness of the developed SPH model.

2 SPH modeling

2.1 Basic theory of SPH

The SPH method is based on a concept of smooth interpolation. Any general field quantity f(x) is approximated by an integration of f(x) multiplied by a weighting function W over a support domain. That is

                  (1)

where Ω represents the support domain at point x; h in the weighting or kernel function W denotes the smoothing length, which is used to define the influence region of W. Furthermore, the gradient of f(x) is formulated by replacing f(x) in Eq. (1) with  and implementing the integration by parts and divergence theorem.

             (2)

where  denotes the gradient of the kernel function W. As the kernel function W possesses the Dirac property and is also satisfied with the compaction and normalization conditions, the smooth interpolation is of second order accuracy. Many SPH kernel functions have been proposed [7,8]. The present work employed the cubic spline kernel function [7] owing to its very good  smooth property and narrow compact domain. It is formulated as

              (3)

where q=(x-x′)/h; α_d is a constant parameter which is dependent on the dimension of the problem. It is chosen as 1/h, 15/(7πh²), 3/(2πh³) for 1D, 2D and 3D simulations, respectively.

In SPH, the continuum is represented by a collection of discrete elements or pseudo-particles. These elements or particles are initially distributed with a specified density distribution and evolve in time according to the conservation equations (mass, momentum, energy, etc.). Flow and energy transport properties at a certain position are determined by the interpolation or smoothing of the nearby particle distribution using the weighting function or smoothing kernel W. We have

               (4)

           (5)

where m_j and ρ_j denote the mass and density of particle j, respectively; N denotes the total number of particles in the compact domain of particle i.

2.2 SPH formulae of foundry filling process

The foundry filling process can be described by a series of conservation equations, including the continuity equation, the N-S equation, and the energy conservation equation, which are summarized below in SPH format.

                           (6)

            (7)

  (8)

where  with  and  denoting speeds of particle; p is the pressure, μ is the dynamic viscosity,  is the shear stain rate; subscripts i and j are indices of particle, α and β are indices of coordinate axes.

In the momentum equation, Eq. (7), and the energy equation, Eq. (8), the symbol Π_ij denotes the artificial viscosity between particles i and j. The artificial viscosity term is important to the stability of SPH calculations, especially on the reduction of non-physical oscillations due to particle compactions at fluid/solid boundaries. The artificial viscosity term was introduced by MONAGHAN [12] and is expressed as follows:

    (9)

where f_ij is defined as

                              (10)

In Eqs. (9) and (10), the constants α_Π and β_Π are two parameters related to the strength of dissipation; c_ij and ρ_ij are averaged sound speed and density, respectively, over all the particles i and j; r_ij is the distance between particles i and j.

In the present work, the weakly compressible method [7] is employed to calculate the dynamic pressure p, in which the following equation of state for fluid is used:

                          (11)

where p₀ is the initial pressure under the reference density ρ₀, γ is chosen to be 7 to guarantee that the maximum density oscillation of fluid is less than 1% compared with the reference density.

The positions of SPH particles relate with the velocity field as

                                   (12)

The primary variables of the system, including the particle density, velocity, position, and internal energy, are obtained from the time integration of Eqs. (6)-(8) and Eq. (12). In the present study, the leapfrog algorithm [7] is adopted for the time integration due to its low requirement on memory storage and high computational efficiency. To ensure sufficient stability of SPH calculations, the time step is controlled by the Courant-Friedrichs-Lewy (CFL) condition [7]. Note that in the present work as a first step of model development, we only focus on the fluid flow and do not solve the energy equation (Eq. (8)).

2.3 Treatment of solid boundaries

Two types of boundary particles were set to mimic the solid boundaries, as schematically displayed in Fig. 1. Type I particles (cross-symbolized in Fig. 1) are placed exactly on the boundary surface and exert Lennard- Jones force on nearby fluid particles. Several layers of Type II particles, denoted with gray solid circles in Fig. 1, are fixed inside the solid mold close to the fluid domain to act as ghost particles. Setting the number of layers of ghost particles depends on the compact support domain of fluid particles. Type II particles have the same spatial separation, Δx, as the fluid particles, while that of Type I particles is Δx/2. The inner node approach is used to position the fluid particles and Type II ghost particles, i.e., the particles adjoining the boundary have Δx/2 distance away from the boundary. Note that in some cases such as for calculations of temperature and stress in the mold, Type II particles can act as true SPH particles too. This composite treatment to the solid boundary ensures the implementation of solid boundary conditions imposed by the foundry mold.

Fig. 1 Treatment of solid boundaries

To determine the properties of Type II particles, detecting points (triangle-symbolized in Fig. 1) associated with Type II particles are set in the fluid region. The detecting points are generated by following a projective method, exemplified for a two-dimensional problem as: 1) searching the two nearest Type I particles (i.e., particle D and E) to a Type II particle A; 2) projecting the particle A into the fluid region with respect to the boundary surface and getting the detecting point A′. During the projection operation, a criterion based on the inner-product of vectors relating the particles A, D and E is used. For the case described by Fig. 2(a), the inner-products d^DA·d^DE and d^EA·d^ED are both greater than zero, the detecting point of particle A is projected by following the line-symmetry with respect to line DE. If one of the inner-products, say d^DA·d^DE is less than zero (see Fig. 2(b)), the detecting point of particle A is projected by following the central-symmetry with respect to particle D. Similar treatment will be performed if d^EA·d^ED is less than zero. Obviously, both the inner-products will not be less than zero at the same time. Using the central-symmetry projection to generate detecting point generally avoids positioning detecting point into the solid region no matter how the boundary is deformed.

Fig. 2 Projective method for generating detecting points

The detecting points are position-fixed and do not participate in the kernel approximation of nearby fluid particles while their properties are updated at each time step in terms of smooth interpolation of nearby fluid particles. We calculate the primary variables such as pressure, density and velocity of a detecting point and then assign the values to the associated Type II particle.

        (13)

      (14)

             (15)

The index j denotes the fluid particles that are in the support domain of the detecting point A′. The calculated  can be decomposed into the normal component  and the tangential component . Assignment of  to the velocity of Type II particle A is dependent on the wall boundary condition assumed. For the free-slip wall boundary condition, the velocity of particle A is assigned to be , while for the no-slip condition, it is assigned as . It is worth pointing out that the normal vector  always points to the fluid domain, which also facilitates the implementation of Neumann and Robin boundary conditions.

Type II particles will participate in the kernel approximation of fluid particles near the boundary in each time step. However, for some high velocity flows, these Type II particles are not enough to fully prevent the fluid particles from penetrating the solid boundary. Using Type I boundary particles, which exert Lennard-Jones force on the nearby fluid particles, can guarantee the no-penetration boundary condition imposed by the solid boundary. The Lennard-Jones force is expressed as [7]

         (16)

where B₀ is a problem dependent parameter and has the unit of squared velocity; parameters a and b are usually taken as 12 and 4, respectively; the cutoff distance r₀ is an important parameter, which determines the influence distance of the repulsive force. In the present implementation, r₀ is set to be one quarter of the initial particle separation and B₀=20. Therefore, the fluid particles interact with Type II boundary particles normally, once a fluid particle approaches the boundary within a distance less than a quarter of the initial fluid particle separation, the Lennard-Jones repulsive force from Type I particles will act to prevent the fluid particle from penetrating the solid boundary.

LIU et al [22] proposed a coupled dynamic solid boundary treatment method, which is seemingly similar to the composite method detailed above. However, there exist major differences between the two method. 1) The repulsive force is formulated differently; 2) values of variables for ghost particles are obtained with in-situ Shepard or moving least square method in Ref. [22]; 3) the present projection method illustrated in Fig. 2 is unique, which is tested to be particularly suitable for modeling the foundry filling process of complex mould geometry.

2.4 Establishment of inlet fixed flow rate boundary

In the foundry filling process, fluid is injected at a fixed flow rate from the inlet into the mold. To establish the fixed flow rate boundary, we set four layers of virtual particles in the inlet region. The velocities of virtual particles are calculated from the imposed flow rate. The positions of the virtual particles are updated while their properties are kept constant during simulations. Once the inlet particles move forward at a displacement of Δx, we convert the forefront layer of inlet boundary particles as actual SPH fluid particles and one new layer of inlet particles is generated as the backmost layer of inlet boundary particles. This method can also be implemented in variable speed injection by setting the velocities of virtual particles as time-dependence.

2.5 Neighboring particle search

In SPH, interactions between fluid particles must be calculated at each time step. Commonly, the search of neighboring particles is based on a fictitious grid of Cartesian cells for two dimensional calculations [7]. The Cartesian cells are ah×ah square regions with a=2 for the specially used kernel function expressed by Eq. (3). The number of total searches is on the order of O(N), with N being the total number of all the fluid and boundary particles.

For foundry filling process, the fluid is gradually injected into the mold and some portions of the mold geometry may not have fluid particles, we upgrade the total SPH particles N at each time step. That is to say, during the search of neighboring particles, we exclude the Cartesian cells that are empty of fluid particles. In this way, the calculation efficiency especially for simulating early stages of the filling process is greatly improved.

3 Physical experiments

To validate the developed SPH model of foundry filling process, water analog experiments are performed. The experimental setup is schematically shown in Fig. 3. Colored water is pumped into a transparent mold made of perspex. A valve and a flowmeter are set in the pipeline and used to control and monitor the water flow rate. Figure 4 details the mold geometry. The fluid is injected from a gate located at the center of the mold bottom plane. A reversed T-shape obstacle is positioned inside the mold. In the present work, two flow rates with fluid injection velocity v=0.45 m/s and v=2.1 m/s, respectively, are taken for experiments. The dynamic fluid filling process is captured by using a high speed digital camera.

Fig. 3 Schematic diagram of experiment setup

Fig. 4 Geometry and dimensions of foundry mold (unit: mm)

4 Results and discussion

Physical properties of the filling fluid are tabulated in Table 1. The initial particle separation Δx is set to be 1 mm. Fluid particles are generated from the bottom gate (referring to Fig. 4) with inlet fixed flow rate boundary prescribed at the gate.

Table 1 Physical properties of fluid

4.1 Low speed filling process

Experimental results at four selected time instants are presented in Fig. 5. The images show that the foundry filling process is stably progressing. The SPH-calculated fluid particle positions and velocity distributions in the foundry mold at these four selected time instants are shown in Fig. 6, and the SPH-calculated pressure distributions are shown in Fig. 7. It is seen that the simulation results successfully capture the major features of free surface flow during low speed filling. Simulation results agree well with the experimental observations.

During early period of the low speed filling process, the injected flow gradually spreads around. Due to the action of gravity, the fluid forms a mushroom-like shape (see Fig. 6(a)). Fluid velocity at the top position of the free surface is zero, indicating an equal distribution of fluid to the horizontal directions. With the progress of the filling process, after the fluid reaches the side walls of the mold, the free surface becomes wavy (see Fig. 6(b)). With further progress of the filling process, once the free surface rises up to the position of the reversed T-shape obstacle, the free surface near the obstacle becomes concave (see Fig. 6(c)). After the fluid level surpasses the obstacle, the free surface is raised steadily and at 2.75 s during the filling process, the free surface approximately reaches half height of the mold (see Fig. 6(d)). Throughout the process, the pressure of fluid distributes smoothly without any pressure spurs in the flow field (see Fig. (7)), indicating very good stability of the calculation and the suitability of the composite treatment to solid boundaries.

4.2 High speed filling process

Experimental images at four time instants from the high speed filling process are shown in Fig. 8. Differently with low speed filling process, the fluid with high injection velocity possesses larger kinetic energy and exhibits more complex phenomena. The injected fluid directly strikes onto the downward surface of the reversed T-shape obstacle and splashes into two streams (see Fig. 8(a)). These two streams, called first grade streams, hit onto the two side walls of the mold, respectively. Each stream splits into two second grade streams, which flow upward and downward, respectively (see Fig. 8(b)). Subsequently, the two downward second grade streams travel back to the central region of the mold and merge with the main stream. During the process, some air is entrapped and two pores form in the bottom region of the mold. Due to the effects of gravity, the two upward second grade streams gradually slow down and fall back to the central region of the mold once the vertical velocity reduces to zero, leading to the formation of another two pores. Finally, there appear four pores inside the flow field (see Figs. 8(c) and (d)). The SPH-calculated velocity field and pressure field during high speed filling process are displayed in Figs. 9 and 10, respectively, both agreeable to experimental observations, which indicates the success of SPH modeling.

Fig. 5 Low speed filling experimental results

Fig. 6 SPH-calculated velocity distribution during low speed filling process

Fig. 7 SPH-calculated pressure distribution during low speed filling process

Fig. 8 High speed filling experimental results

Fig. 9 SPH-calculated velocity distribution during high speed filling process

Fig. 10 SPH-calculated pressure distribution during high speed filling process

5 Conclusions

1) SPH modeling of foundry filling process needs to properly handle the involved solid boundaries. In the present SPH model development for foundry filling process, a composite solid boundary treatment method is elaborately designed, which combines two solid boundary treatment measures. One is the L-S force type particles set exactly on solid boundary surfaces; the other is ghost particles located inside the solid mold region.

2) SPH simulations to two foundry filling processes of differing fluid filling rates are meaningful and reasonable results, demonstrating the successful implementation of SPH numerical methodology.

3) Water analog experiments were carried out, which gave results according well with SPH simulations, corroborating the validity and accuracy of the developed SPH model.

References

[1] SCHMID M, KLEIN F. Fluid flow in die-cavities: Experiments and numerical simulation [C]//Proceedings of the 18th International Die Casting Congress and Exposition. Indianapolis, 1995: 93-97.

[2] KERMANPUR A, MAHMOUDI S H, HAJIPOUR A. Numerical simulation of metal flow and solidification in the multi-cavity casting moulds of automotive components [J]. J Mater Process Tech, 2008, 206: 62-68.

[3] TOH T, TAKEUCHI E, YAMAMURA K. Multi-physics analysis for electromagnetic processing in continuous casting of steel by finite volume method approach [C]//Proceedings of the Seventh International Conference on CFD in the Minerals and Process Industries. Melbourne, 2009: 1-5.

[4] KIM K D, YANG D Y, JEONG J H. Adaptive refinement techniques based on tetrahedral and hexahedral grids for finite element analysis of mold filling in casting processes [J]. Comput Method Appl M, 2006, 195: 48-49.

[5] ZHAO Hai-dong, BAI Yan-fei, OUYANG Xiao-xian, DONG Pu-yun. Simulation of mold filling and prediction of gas entrapment on practical high pressure die castings [J]. Transactions of Nonferrous Metals Society of China, 2010, 20(11): 2064-2070.

[6] PANG S Y, CHEN L L, ZHANG M Y, YIN Y J, CHEN T, ZHOU J X, LIAO D M. Numerical simulation two phase flows of casting filling process using SOLA particle level set method [J]. Appl Math Model, 2010, 34: 4016-4122.

[7] LIU G R, LIU M B. Smoothed particle hydrodynamics—A meshfree particle method [M]. Singapore: World Scientific Publishing, 2003: 33-149.

[8] LIU M B, LIU G R. Smoothed particle hydrodynamics (SPH): An overview and recent developments [J]. Arch Comput Methods Eng, 2010, 17: 25-76.

[9] HU X Y, ADAMS N A. An incompressible multi-phase SPH method [J]. J Comput Phys, 2007, 227: 264-278.

[10] MORRIS J P. A study of the stability properties of smooth particle hydro-dynamics [J]. Publ Astron Soc Aust, 1996, 13: 97-102.

[11] MOLTENI D, COLAGROSSI A. A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH [J]. Comput Phys Commun, 2009, 180: 861-872.

[12] MONAGHAN J J. Simulating free surface flows with SPH [J]. J Comput Phys, 2003, 110: 399-406.

[13] COLAGROSSI A, LANDRINI M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics [J]. J Comput Phys, 2003, 191: 448-475.

[14] JIANG F M, OLIVEIRA M S A, SOUSA A C M. Mesoscale SPH modeling of fluid flow in isotropic porous media [J]. Comput Phys Commun, 2007, 176: 471-480.

[15] JUBELGAS M, SPRINGEL V, DOLAG K. Thermal conduction in cosmological SPH simulations [J]. Mon Not R Astron Soc, 2004, 351: 423-435.

[16] CLEARY P, HA J, PRAKASH M, NGUYEN T. 3D SPH flow predictions and validation for high pressure die casting of automotive components [J]. Appl Math Model, 2006, 30: 1406-1427.

[17] CLEARY P, HA J, PRAKASH M, NGUYEN T. Short shots and industrial case studies: Understanding fluid flow and solidification in high pressure die casting [J]. Appl Math Model, 2010, 34: 2018-2033.

[18] PRAKASH M, CLEARY P, GRANDFIELD J. Modelling of metal flow and oxidation during furnace emptying using smoothed particle hydrodynamics [J]. J Mater Process Tech, 2009, 209: 3396-3407.

[19] CHANG T J, KAO H M, CHANG K H, HSU M H. Numerical simulation of shallow-water dam break flows in open channels using smoothed particle hydro-dynamics [J]. J Hydrol, 2011, 408: 78-90.

[20] ADAMI S, HU X Y, ADAMS N A. A generalized wall boundary condition for smoothed particle hydrodynamics [J]. J Comput Phys, 2012, 231: 7057-7075.

[21] HE Yi, ZHOU Zhao-yao, CAO Wen-jiong, CHEN Wei-ping. Simulation of mould filling process using smoothed particle hydrodynamics [J]. Transactions of Nonferrous Metals Society of China, 2011, 21(12): 2684-2692.

[22] LIU M B, SHAO J R, CHANG J Z. On the treatment of solid boundary in smoothed particle hydrodynamics [J]. Sci China Tech Sci, 2012, 55: 244-254.

铸造充型过程光滑粒子流体动力学建模及数值模拟

曹文炅 ^1,2，周照耀 ²，蒋方明 ¹

1. 中国科学院 广州能源研究所，中科院可再生能源重点实验室，广州 510640；

2. 华南理工大学 机械与汽车工程学院，广州 510640

摘  要：建立了模具充型过程的光滑粒子流体动力学方法数值模型。为了表征模具固壁对充型流体的约束，建立了一种耦合固壁边界条件。在固壁表面设置了L-S斥力粒子；在固壁内部设置了若干层边界虚粒子，实现对流体粒子支持域的补充。通过冷态水模拟实验对计算模型进行了验证。计算结果与实验吻合良好，证明了本文模型的有效性。

关键词：光滑粒子流体动力学；铸造充型过程；耦合固壁边界；水模拟实验

(Edited by Yun-bin HE)

Foundation item: Project (2011006B) supported by the Open Project of National Engineering Research Center of Near-Shape Forming for Metallic Materials, China; Project (FJ) supported by the CAS “100 talents” Plan

Corresponding author: Fang-ming JIANG; Tel: +86-20-87057656; E-mail: fm_jiang2000@yahoo.com

DOI: 10.1016/S1003-6326(15)63847-X