On the CFD Modelling of Slamming of the Metal Melt in High-Pressure Die Casting Involving Lost Cores

This paper uses computational fluid dynamics (CFD), in the form of the OpenFOAM software package, to investigate the forces on the salt core in high-pressure die casting (HPDC) when being exposed to the impact of the inflowing melt in the die filling stage, with particular respect to the moment of first impact—commonly known as slamming. The melt-air system is modelled via an Eulerian volume-of-fluid approach, treating the air as a compressible perfect gas. The turbulence is treated via a Reynolds-averaged Navier Stokes (RANS) approach. The RNG k-ε and the Menter SST k-ω models are both evaluated, with the use of the latter ultimately being adopted for batch computations. A study of the effect of the Courant number, with a view to establishing mesh independence, indicates that meshes which are finer, and time steps that are smaller, than those previously employed for HPDC simulations are required to capture the effect of slamming on the core properly, with respect to existing analytical models and empirical measurements. As a second step, it is then discussed what response should be expected when this force, with its spike-like morphology and small force-time integral, impacts the core. It is found that the displacement of the core due to the spike in the force is so small that, even though the force is high in value, the bending stress inside the core remains below the critical limit for fracture. It can therefore be concluded that, when assuming homogeneous crack-free material conditions, the spike in the force is not failure-critical.


Introduction
High-pressure die casting (HPDC) is an important process for the manufacture of high-volume and low-cost automotive components, such as automatic transmission housings, crank cases and gear box components [1][2][3]. Liquid metal, generally aluminium or magnesium, is injected through a complex system of gates and runners and into the die at speeds of between 50 and 100 ms −1 at the ingate, and under pressures as high as 100 MPa [4].
From an economic point of view, HPDC is typically constrained by a huge base investment in machinery and tooling, although this is alleviated by low incremental costs for each additional unit produced; in short, the process scales very well with increasing output. On the other hand, this complicates the task for the design engineer, who has to be sure about the viability of a process and manufactured parts before any budget is invested in tooling and machinery. One technological constraint to date is that there is no serial production, via HPDC, of parts with inlying hollow shapes or undercuts formed by lost cores, even though other casting techniques have employed lost cores for decades [1,4].
Nevertheless, several ideas for such products do already exist [5][6][7], and one material that has been put forward as a candidate for HPDC with lost cores is salt [6,8,9]. The basic idea of using salt cores is to block parts of the die volume by inserting cores as placeholders; in so doing, the melt will not penetrate into this space. The cores may then be removed after solidification and one creates undercuts or hollow sections with them, which may then later act as cooling or oil-flow channels [6,8,9].
One way to determine whether lost cores are indeed a viable option for a given geometry is to employ numerical simulation-in particular, computational fluid dynamics (CFD) [10][11][12][13]. However, whereas the papers just mentioned focused on the flow of a melt past a core in more general terms, the particular focus of this paper shall be on determining the load on the solid core during the melt's initial impact. This concept is often referred to as slamming, and is a phenomenon that is investigated more commonly in the context of naval architecture [14][15][16]; it occurs when a two-phase interface with one phase of high density -water -and another of lower density -air-hits an obstacle, for example, a ship. This can also be the case when investigating lost cores in high-pressure die casting, when the air surrounding the core is replaced by the much heavier melt. The particular scientific and technological issue of interest is the initial peak in the load that all previous simulations of high-pressure die casting and lost cores have predicted [10][11][12][13]; more precisely, we wish to investigate the core's response to this load peak in more detail, in order to determine whether the core will fail as a result of it.

Model Equations
As in our earlier CFD modelling of high-pressure die casting [12,17], we model the two-phase flow of molten metal and air by using the volume-of-fluid (VOF) method [18], wherein a transport equation for the VOF function, γ, of each phase is solved simultaneously with a single set of continuity and Navier-Stokes equations for the whole flow field; note also that γ, which is advected by the fluids, can thus be interpreted as the liquid fraction. Considering the molten melt and the air as Newtonian [19], compressible and immiscible fluids, the governing equations can be written as [20,21] where t is the time, U the mean fluid velocity, p the pressure, g the gravity vector, F s the volumetric representation of the surface tension force and T denotes the transpose. In particular, F s is modelled as a volumetric force by the Continuum Surface Force (CSF) method [22]. It is only active in the interfacial region and formulated as F s = σκ∇γ, where σ is the interfacial tension and κ = ∇ · (∇γ/|∇γ|) is the curvature of the interface; the significance of the term U r in Equation (3) is explained shortly in Section 2.2. The material properties ρ and µ are the density and the dynamic viscosity, respectively, and are given by where the subscripts g and l denote the gas and liquid phases, respectively. We take ρ l , µ g and µ l to be constant, but assume the air to be an ideal gas, that is, its density changes with pressure and temperature; hence, the equation of state for our model reads where M is the molecular weight for air, R is the universal gas constant and T is the temperature. This is a new unknown in the system and needs to be solved for via the heat equation [20,21], where K = 1 2 U · U is the kinetic energy, c v g and c v l denote the specific heat capacities at constant volume for the gas and liquid phases, respectively, α e f f is given by where k g and k l denote the thermal conductivities for the gas and liquid phases, respectively, and σ tur is the turbulent Prandtl number, whose value is set to 0.9 [23]. Note that α e f f resembles a phase-averaged thermal diffusivity that includes the contribution of turbulence, although it lacks a density term in the denominator. Furthermore, µ tur in equation (2) denotes the turbulent eddy viscosity, which will be calculated via the Menter SST k-ω model [24]; the equations for this model are summarized in the appendix. The major purpose of this paper is to calculate the forces on the core during slamming. Those are governed by the following formulae. For i = x, y, the total force acting on the core in direction i, F i,tot , is given by where F i,p and F i,s are pressure and viscous shear forces, respectively, and are given by where U = U x , U y as the velocity and its components and A c is the surface area of the core. The model two-dimensional (2D) geometry shown in Figure 1. This consists of two channel walls, an inlet, an outlet and a salt core that is located about the horizontal axis of symmetry of the channel. A summary of the mathematical expressions of the applied boundary conditions is presented in Table 1. In summary: the channel walls and the surface of the core are no-slip, zero heat-flux boundaries; the inlet velocity and temperature are prescribed; the outlet is at ambient pressure and heat conduction is negligible. In addition, there is always melt at the inlet, so that γ = 1 there; at all other boundaries, γ satisfies a zero flux condition.
As for the initial conditions of the die, the cavity is filled with air (γ = 0) that is at rest (U = 0), at ambient temperature (T = T amb ) and at atmospheric pressure (p = p amb ). The initial conditions for the turbulence model are set via the length scale of the largest eddies and the turbulence intensity. Those quantities were set to 2 mm and 5 %, respectively.
What is now still missing are the material properties of the phases. Those are presented in Table 2. Note that we have taken ρ l and µ l to be constant, even though they would in general depend on the temperature. The reason for taking them to be constant is that the temperature of the melt will not change appreciably from T melt during the filling, in view of the zero heat flux boundary conditions that are being applied; on the other hand, ρ g may change appreciably, as the gas will heat up from its initial ambient temperature and experience a significant increase in pressure, both of which are accounted for via equation (6). One additional comment concerns the specific heat capacities c p g and c p l : the index p means that their values have been measured and documented for constant pressure, which are more commonly tabulated rather than c v g or c v l . Equation (7), however, requires c v . This conversion can be conducted via the value of the isentropic expansion factor or heat capacity factor, which can for air be assumed to be constant and equal to c p g /c v g = 1.4 in the given temperature range [12,25,26]. For the aluminum melt, this additional step is not necessary, as it was treated as being incompressible; thus, c v l = c p l . Table 2. Model parameters. The parameters for gas are those for air; those for metal are for the alloy AlSi9Cu3 [12,19,27].

Simulation Methodology
To solve the above equations, the OpenFOAM software package [28][29][30][31] was used, due to the niche nature of the field of application, as well as to make use of its extendability. OpenFOAM has most of its applications designed for solving fluid mechanics problems [29,32] and applies the finite volume method [23,33,34] to solve them. With the later aim of scaling the model to complex industrial 3D-geometries in mind, OpenFOAM is also a very powerful tool, since the distribution of a case to a huge number of CPU-cores is not limited due to license restrictions, as is the case with commercial CFD codes. In particular, we used the compressibleInterFoam solver, which handles two compressible, non-isothermal immiscible fluids using a VOF interface-capturing approach. The solver also uses the multi-dimensional universal limiter with explicit solution (MULES) scheme for interface compression, which introduces a supplementary velocity field, U r , in the vicinity of the interface, as seen in equation (3); in doing so, the local flow steepens the gradient and the interface becomes sharper and more pristine. A typical form for U r is U r = min(U, max(U)), as given by Reference [20]. Further benefits that accrue from using OpenFOAM for this application that its implementation of the Menter SST k-ω model has previously been shown to be robust, and to give results that are in excellent agreement with experimental data [35,36]. Moreover, by use of a function object called forces that is already pre-implemented in OpenFOAM, the forces required in equations (10) and (11) can be readily calculated.
The corresponding mesh for the 2D-geometry shown in Figure 1 is illustrated in Figure 2. This is a structured mesh consisting of one layer of hexahedral cells. It is also parameterized and can therefore easily be refined for later mesh independence studies. For this purpose, the OpenFOAM utility refineMesh was used that, by default, halves the mesh spacing upon execution, that is, quadrupling the cell count for a 2D mesh as in Figure 2. This simplified model was used for efficiency reasons in order to establish that the newly installed tool yields sensible results. Lastly, as this is a rather complex flow problem in the context of CFD research, it is worthwhile to look at the solution algorithm a little more closely. In order to do that, Figure 3 is presented. Here, the standard PIMPLE algorithm, a combination of PISO (Pressure Implicit with Splitting of Operator) [37] and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) [38], was modified and, for this purpose, an additional step was added that solved the transport equation of the phase variable γ according to equation (3) at the beginning of the solution algorithm, that is, before the momentum predictor. Figure 3 shows the original PISO algorithm on the left and on the right the reader sees the implemented extension of it including the additional step for the phase field for handling the two-phase flow. The corrector loops can either be broken by a residual or a maximum number criterion. One additional comment for the sake of completeness is that the outermost loop for the increasing time step is omitted, as the chart only shows the processes that occur within each time step.

Results and Discussion
In this section, we present the results of the numerical simulations. The principal aim of these is to determine the magnitudes of the forces which act on the salt core, with a view to being able to infer whether an actual salt core would be able to withstand such forces in practice. As one of the diagnostics of the computations, we consider first the slamming factor, so as to link our results to analytical expressions for this quantity that are available in the field of marine technology.

The Concept of the Dimensionless Slamming Factor
Following the results of recently published articles on salt cores in high-pressure die casting [10,12,13,17], in almost every simulation where the core was treated as rigid, the signature peak mentioned in Section 1 appeared in a plot of the force exerted on the core as a function of time. The ambition here was therefore to investigate its nature in more detail. This will be done by introducing the dimensionless slamming factor, C s , which is defined as where F is the force, and is given by with F x,tot and F y,tot being computed according to Equations (9) to (11); C s also depends on the density, ρ, the undisturbed melt velocity U, the radius of the core, R, and its length, l. Analytical progress was made by von Karman [39] and Wagner [40] to determine this slamming factor in the context of marine technology, with the result that according to von Karman the slamming factor for a cylinder was π, while Wagner put it at 2π. More recent studies [14,41], however, indicate that the value according to von Karman represents a lower limit of the slamming factor, while the Wagner model acts as an upper cap. Both models also fail to provide an evolution of the value over time. Campbell and Weynberg [15,42] present a model that offers results based on empirical data to determine the plot over time, according to the formula The letters in Equation (14) represent the same quantities as in Equation (12). The introduced dimensionless slamming factor makes the calculated results independent of the ingate velocity or geometric dimensions. This is further illustrated in the following plot. Figure 4 shows that the calculated slamming factor is the same no matter whether it is being evaluated for an ingate velocity of 10 or 20 ms −1 . Those values for U in were used for U for calculating the slamming factor according to Equations (12)- (14).

Determining the Slamming Factor with Respect to Mesh Resolution
The model setup was calibrated with existing data on similar cases, such as the water entry of an axisymmetric projectile [43] with similar Reynolds numbers and also for the feature of the pile-up effect [41,44,45], wherein a jet is formed near the intersection with the initial water surface.
Applying the model to the previously introduced case, Figsures 5-7 show the temporal development of the phase distribution, the pressure and the velocity at impact on the core. Please note that, although the model is non-isothermal, it does not take the temperature difference between the core, the walls and the melt into account, in the sense that zero heat-flux boundary conditions are applied at all solid boundaries (see Table 1). The actual fluid flow and heat transfer in the near-wall area may therefore differ in reality due to heat transfer-induced effects.   After calibrating the setup, it was time to proceed to evaluating the forces on the core and to benchmark them with previously published data. The first thing to determine was the necessary mesh resolution in order to compute plausible results for the slamming factor. The results of this mesh resolution study are presented in Figure 8. A meticulous investigation of this phenomenon of slamming showed that, with generally used meshing standards in industry-oriented computations and mesh spacings in the range of 1 mm [10], the CFD simulation underpredicts the value of the slamming factor (see Figure 8). It becomes evident that the simulation requires a mesh as fine as a spacing of 0.3 mm in order for the computed result for C s to be above the lower limit of the slamming factor according to the von Karman model. This is significantly lower than the typical mesh spacing that was used in Reference [10] and requires a mesh for this rather small geometry to consist in total of 60,000 cells, albeit being a 2D mesh. Transferring these results onto a real-world casting geometry that includes salt cores and designing a 3D mesh according to these findings would, in turn, produce a mesh that it would be completely impractical to use with currently available computational capacity.  [39] and Wagner [40]. Here, U in = 20 ms −1 .
While Figure 8 showed that coarser meshes underpredict the peak in the force, the stationary value of the force tends to be overpredicted the coarser the mesh is. From a design engineer's point of view, this is a beneficial outcome as it includes a built-in safety factor into the design of a devised product. Note also that the stationary or steady-state solution in Figure 8 refers to the constant value for C s that is obtained once the melt front is well past the salt core; see, for example, Figure 4. Figure 9 shows this study's result in comparison with the findings of previously published articles [15,39,40,46,47]. The results are plotted for a mesh spacing of 0.025 mm. Technically, only a refinement of the cells in the near core area would be necessary. However, given the fact that then the structured nature and the numerical benefits that come along with it would have to be sacrificed, this opportunity was forgone here as the domain was sufficiently small and calculation efficiency was not a benchmark in this academic study. It should be noted that the analytical models start with t = 0 at the point of impact, whereas for our numerical simulation t = 0 denotes the time a which melt start to enter the cavity; for this reason, the simulation time values were therefore adjusted so that t = 0 in Figure 9 corresponds to the time when melt first hits the core, with the consequence that negative values of t appear in this figure.  Figure 9 shows also the values for several other models, among them the static von Karman [39] and Wagner [40] models, as well as more recent models that also represent the development of the force over time [46,47]. It can be concluded that this study's assessment of the slamming factor agrees very well with earlier authors' findings. One interesting observation of the result from the present volume-of-fluid model is that its increase is, compared to the mostly analytical models, more gradual. The other models feature a sharp step-like increase in value, resembling the morphology of a Heaviside function [48].
The reason for this gradual increase is to be found most likely inside the volume-offluid approach itself. As Equations (4) and (5) illustrated, the material properties are being averaged by their volume fraction, γ. Another reason for this feature is to be found in the numerical approach of this study, where the numerics are known to smear out the interface of volume-of-fluid simulations [49] and therefore artificial interface compression terms are to be introduced into the partial differential equation for the volume fraction [12,49], equation (3). This effect of a non-discrete interface causes additional problems in stability of the interFoam solver-family in the sense of spurious parasitic velocities at the interface, and several strategies to mitigate the effects towards a sharper interface have since been proposed [50][51][52][53]. We can conclude that, taking into account Figs. 8 and 9, one has, for the given boundary conditions, to maintain at least a mesh fineness corresponding to a spacing of 0.2 mm, in order to reach a value for C s above the limit defined by von Karman in Reference [39], if one wants to achieve a proper result for the slamming factor. Transferring this finding to more applicable designs in industry, it underscores the need for mesh refinement near a salt core boundary if the load is to be evaluated correctly.

Effect of turbulence
The aspect of turbulence treatment on the slamming factor in the simulations was also evaluated during the investigations. Figure 10 features this result. Two different turbulence models were benchmarked against each other for the given case and the different mesh spacings, as illustrated in Figure 10. The benchmarked turbulence models are Menter's SST k-ω model [24,54] and the RNG k-ε model [55]. The reason why those two were picked is that the Menter model was proven to be of excellent stability and accuracy proven in studies published by other authors [35] and our own results published in previous papers [11,12]. The RNG-model was selected as it is the standard of the commercial CFD casting software Flow-3D Cast. Both turbulence models belong to the RANS-model family (Reynolds-Averaged-Navier-Stokes equations). As is evident from Figure 10, the selection of the turbulence model is of minor importance. Even with different mesh spacings, the result of the slamming factor stayed more or less the same. This is plausible, as the main contribution to the slamming results from the pressure forces on the core-an effect documented in more detail in Reference [12]. For those forces, the momentum of the melt is more important compared to other factors in the fluid. This is also easy to comprehend with the introduced formulae in Equations (10) and (11). One sees here that the only quantity that is influenced by the turbulence model is µ tur , while all the others do not directly depend on it. This turbulent eddy viscosity, however, only appears in the formula for the viscous forces, which are of negligible importance as seen in Figure 4.

Effect of Courant number
One might now argue that the presented strong mesh dependency of the results is the consequence of the simulation time-step size being Courant number-controlled [56], and hence that an automatic reduction of the mesh spacing should lead to a shorter time step, ∆t, in consequence. In order to eliminate this possibility, studies with different Courant numbers were also conducted; in particular, each cell's local Courant number is calculated, and the solver then sets a value for ∆t so that these local Courant numbers are not greater than some maximum allowable value that is specified by the user. The result of this investigation is shown in Figure 11. In Figure 11a, the values for F stat are shown, whereas Figure 11b shows the values for F max . F stat calculates the slamming factor value via the average force computed while the denser phase surrounds the cores, while F max uses the peak force value. Although the latter in particular oscillates somewhat, there is no obviously visible trend that would support the hypothesis that only reducing ∆t would also lead to a more physical value for the slamming factor on a coarser mesh. The values for F stat and F max are the same for a Courant number range from 0.1 to 2, while at the same time the difference caused by the different mesh resolution was clearly visible. In conclusion, a mesh spacing of 0.6 mm with a Courant number of 2 leads to values much closer to reality for the peak in the force than a mesh resolution of 1 mm with a maximum Courant number of 0.1. But, even for a mesh resolution of 0.6 mm, the computed values for the peak in the force are significantly below the lower bound given by the von Karman model (π).

Response of the Core to the Spike-Like Force Impact
Although it is academically interesting to determine the value for the slamming factor as precisely as possible, the general underlying question is also whether such a small force-time interval is failure-critical for salt cores in high-pressure die casting in general. We will therefore in the following present a couple of considerations to assess the necessity for the engineer to keep the slamming factor below a certain level.
If one, for this purpose, leaves the field of 2D considerations and imagines how a beam of length 70 mm that is situated within the die facing the inflowing melt stream orthogonally would react to the inflowing melt and the slamming, one would make the following observations. The consideration is of course based on the assumption that the inflowing melt stream stays as thin as in the test case, a fact that has emerged in more application-oriented studies [12,17]. In general, any force impacting on the core at one point cannot travel through the solid body faster than the speed of sound, with v s being the speed of sound in salt, K being the bulk modulus for salt cores and ρ s being the density of salt cores. K is related to the documented quantities for salt cores [17] via References [57,58] with E as the Young's modulus and ν as the Poisson ratio. Assuming ρ s as 2056 kgm −3 , ν as 0.21 and E as 1.5×10 10 Pa [17], Equations (15) and (16) yield the speed of sound in salt as 2048 ms −1 . Based on Figure 4, one can now assume 0.02 s as a scale for the dimensionless filling time, yielding 3.5×10 −6 s as an absolute time scale during which the peak force acts on the core. Together with the calculated speed of sound, this means that a stress signal travelling a distance below 10 mm through the core, that is, referring back to the beam with 70 mm in length, would not even be able to reach the bearing.
If we further assume a slamming factor of 3 2 π, as reported in Figure 9, the absolute force according to Equation (12) is slightly below 1400 N. Now, applying Newton's second law, where F is the force, m is the mass and a is the acceleration, integrating (17) twice with respect to t, and assuming the initial velocity to be zero and a to be constant, to obtain where s denotes distance, one may get an idea of how far the centre of the core travels. Based on the calculations above for the distance travelled in two directions of the beam and the reported dimensions of the core in Figure 1, neglecting the round edges, one obtains a displaced mass of 0.02 kg leading, with the previously calculated value for the force, to an acceleration of 6.6×10 −6 ms −2 and in turn to a displacement of the core section of 4×10 −7 m in the direction of the mean flow, according to equation (18). Three-point bending tests with salt cores, conducted by the authors [17], show that the tested specimens manufactured via pressing and sintering [59] are, even at room temperature, capable of undergoing displacements of this scale without cracking, provided that they are free of internal cracks, that is, meaning that the continuum mechanics approach can be applied. In turn, this leads to the conclusion that, due to the short time scales of the slamming effects, they are not to be considered failure-critical for the salt core. Based on the presented knowledge, only a small segment of the core will be displaced by the spike in the force, causing the core in turn to vibrate. The time scale of the impact will be too small to reach the core bearings and produce a counterforce. Due to the displacement, the core will start to vibrate before being permanently deformed by the steady-state force. Since this steadystate force will be present for a longer period of time, this may in the end be failure-critical. To finally conclude the discussion, the reader should bear two additional things in mind. The results presented, although relying heavily on dimensionless numbers, are hugely geometry-, parameter-and, most notably, ingate/impact velocity-dependent. One must therefore be careful to transfer these findings to other geometries and setups; on the other hand, it is highly recommended to apply the presented CFD approach to the geometry at hand. Secondly, the assumption of the core being crack-free may also be flawed with respect to reality. It may therefore be possible that salt cores, as the material salt is of a ceramic nature, will most often contain small cracks. It thus provides scope for future research work to develop a model for assessing how long potential cracks are allowed to be in order to withstand the slamming impact. Another interesting direction will be to measure not only the displacement of a die casting, as was done in Reference [17], but also the forces on it.
As a practical conclusion, this also suggests that pre-heated cores that are warmer than room temperature may yield a higher fraction of castings with intact cores and thus yield channel geometries within the tolerance limit. Higher temperatures naturally make it more likely that the material will fail in a ductile way, rather than in a brittle way.

Conclusions
In the course of this paper, the presented research work showed that the slamming phenomenon is something that is in general underestimated by state-of-the-art CFD simulations in industry that have used typical mesh spacings of 1 mm. We have found that the mesh resolution has to be increased by a factor of 40 in order for the computed slamming factor to be in the range of existing analytical results of empirical measurements. Interestingly, the reduction of the time step size, ∆t, did not show any impact on the proper estimation of the slamming factor's value. The contribution of the turbulence in this study was found to be negligible and the mathematical justification for this was shown with the presented formulae; this echoes the findings and conclusions of more application-related studies [12]. It remains a topic for scientific debate as regards what response of the core the slamming peak in the force actually causes inside the material. This paper provides a reasoning based on simple continuum mechanics formulae that the core will, for the given parameters, be only slightly displaced and thus in turn only vibrate -thus rendering the slamming impact to be a non-failure-critical phenomenon. The reader is reminded that the provided results are only valid for the presented boundary conditions and transferring them to other setups has to be handled with care. It is never possible, without the detailed geometric constraints and boundary conditions, to decide whether a salt core solution for high-pressure die casting will be viable. Future research will be undertaken to measure the force inside a die during filling and also to combine the concept of slamming with the laws of fracture mechanics.

Data Availability Statement:
The data presented in this study are available on request from the first author.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.