Smoothed Particle Hydrodynamics Simulation of High Velocity Impact Dynamics of Molten Sand Particles

: Sand ingestion is highly detrimental for gas turbines because it leads to erosion and corrosion of engine components, accelerating material fatigue and contributing to global engine failure. In this paper the high velocity impact of a molten sand particle onto a solid wall is investigated by means of the Smoothed Particles Hydrodynamics method where the three phases are taken into account. Nominal conditions are a 25 µ m particle composed of molten sand (dynamic viscosity µ l = 11 Pa · s) impacting the wall at a velocity of 250 m/s. The influence of different parameters are explored such as the mechanical properties of the molten sand particle (density, viscosity, surface tension), the impact conditions (velocity magnitude, particle size and angle of impact) as well as the particle shape (sphere or cube with different geometrical features impacting the wall). It is found that the particles do not form a lamella during the impact but mostly conserve its initial shape. It is also confirmed that sharp features such as edges lead to a larger normal pressure at the impact location. Correlations to quantify (i) the spread factor, (ii) the maximum and mean impact force and impact pressure and (iii) the slip distance are derived for the first time based on the investigated parameters. The importance of these correlations is that they provide information needed to implement low-order models for studying impact and deposition of molten sand in engineering simulations.


Introduction
During mission operations of flying vehicles, gas turbines engines can ingest small solid particles during take-off or hovering, or when flying through a dust storm. The most common types of ingested particulates are volcanic ash, desert sand, dust and salt, in the form of calcia-magnesia-alumino-silicates (CMAS), where the specific composition of particulates often varies with geographical location. When ingested, these solid particles may impact engine components at high velocity, increasing the mechanical fatigue of the protective material and eventually lead to component failure. As particles pass through the engine combustor stage, where the operating temperatures exceed the melting point of sand (1300-1600 • C), they change phase. This leads to an increase in deposition onto engine component surfaces and enhances the potential for chemical reactions with the thermal or environmental barrier coating (T/EBC). In addition, when the ingested material sticks to the surface of components, it can lead to a degradation of the flow path, clogging of fuel nozzles or blockage of cooling holes [1].
Even though typical particles separators can remove particles larger than ≈80 µm, smaller particles pass through and are entrained into the turbine stage. Extensive use of particle separation systems also negatively impacts the overall system performance due to a reported high loss of inlet pressure and inlet mass flow [2]. To date, contaminated environments remain an enduring challenge for the aviation community, not only in desert regions, but also to enable safer flights near complex megacities and urban canyon environments across the world [3,4].
In order to develop superior protective systems that are resistant to CMAS attack, an integrated computational and experimental approach that leads to predictive models and insights of the interaction physics is essential. Over the last two decades multiple experimental research efforts have focused on developing novel protective materials to improve corrosion resistance to CMAS attack [5,6]. A more detailed understanding of how particles behave when they impact surfaces is needed to guide design of T/EBCs and, for example, to precisely define suitable targets for boundary surface energies and to tailor the surface wettability characteristics. Recently, simulations of CMAS particulates have been conducted to study the dynamics of droplet binary collisions and to characterize the turbulent kinetic energy process [7,8]. The results provide critical information needed for the development of CMAS coalescence and breakup models in coarse-grained sand-laden simulations. Bravo et al. [4] conducted analysis of CMAS particulates and their deposition for transonic flows through a linear cascade of blades and reported the sensitivity to the particle Stokes number. Recently, Singh & Tafti [9] developed a low-order deposition model for sand particles at high temperature based on the sticking probability. The sticking probability includes the effects of impact velocity (magnitude and angle) and temperature of the particle. The researchers took into account the change of state of the material due to the kinetic energy dissipation during the impact. They used their model within LES simulations of the impact of a sand-particle-laden jet on a flat coupon and compared it with the experiment [10], showing a satisfactory agreement. Goshal et al. [3] and Murugan et al. [11] conducted numerical investigations on the impact of spherical and sharp edged particles with Computational Fluid Dynamics (CFD) and a software for nonlinear dynamic impact analysis (LS-DYNA). Among their numerous findings, it was shown that sharp edged particles cause more damage than smooth particles and that the particle residence time was sufficiently high for the sand particle to melt. Jain et al. [12] reported the development of a new particle forcing model and conducted Smoothed Particle Hydrodynamics (SPH) simulations of a molten CMAS sessile droplet and computed the contact angle of the droplet on a Yttria-stabilized zirconia (YSZ) material.
The objective of the present work is to simulate the high-speed impact of a molten sand particle under realistic gas turbine operation conditions with a full multiphase approach based on first principles. To the author's knowledge, such a study has not been performed previously. The different mechanisms associated with the deposition of impurities on turbine blades depend on the Stokes number of the particles in the flow over the blade surfaces. For small particulates in flows with a Stokes number lower than one, these mechanisms are turbulent diffusion, eddy impaction, Brownian diffusion and thermophorese [1]. For larger particulates, typically larger than ≈10 µm in realistic gas turbine operation conditions [13], the Stokes number increases and, in turn, particle deposition is dominated by inertial impaction. The present study focuses on the early phase of inertial impaction. Subsequent phenomena, such as deposition, sticking and chemical reactions are not investigated here, and no correlations for spreading or rebound will be applied. The particle deformation and the evolution of the interface will be computed based on mass and momentum conservation. As this study is a first step, the focus will be on the mechanical outcomes of the impact. The energy equation will be neglected, such that the heating of the particle/TBC due to the dissipation of kinetic energy will not be taken into account. The selected numerical approach is the SPH method. It is a mesh-free Lagrangian method where the discretization is performed by elements which are called particles. These particles are moving at the fluid velocity and carry physical quantities such as pressure, mass, momentum and energy. As a Lagrangian particle method, it inherently takes advection and the corresponding acceleration into account. When simulating multiphase flows, each phase is represented by different particle types (e.g., gas, solid, liquid) featuring different physical properties.
All phases are computed with the same formalism. The SPH method was originally developed in the context of astrophysics [14,15]. It is nowadays successfully applied in many fields of fluid mechanics such as free surface flows [16], fluid-structure interactions [17], turbulent flows [18] and compressible flows [19]. In the context of liquid droplet impact on walls, the SPH method was used to investigate the influence of porosity and roughness on sessile droplet spreading [20,21]. The shape of the impinging droplet was investigated [22], and the liquid vaporization was taken into account to study the influence of the wall temperature [23]. The SPH method has also been applied to solid deformation and brittle failure of rocks during impacts [24].
The present paper is structured as follows. The numerical model is presented in Section 2, followed by a validation in Section 3 and the application to real conditions in Section 4.

Method
In the SPH method, any field and its derivative are expressed at a given particle location (subscript a) by an interpolation over its neighbor particles (subscript b) [25]: where V b is the volume of the neighbor particles. The term W is a weighting function, referred to as the kernel. It depends on the interparticle distance r b − r a and a characteristic length scale h (smoothing length). Its role is to promote the influence of closer neighbors (Figure 1, top). Only the neighbors encompassed in the so-called sphere of influence Ω a (Figure 1, bottom) whose radius is proportional to h are accounted for in Equation (1). The kernel used in this study is a quintic spline: where r = |r a − r b |/h is the distance between particle a and b normalized by the smoothing length h and is equal to the initial mean particle spacing ∆x. The radius of the sphere of influence Ω is R = 3h. The terms f (r a ), f (r b ) and W(r b − r a , h) are abbreviated by f a , f b and W ab in the rest of this work. We express the motion of the fluid phases (gas and liquid) by the weakly-compressible Navier-Stokes equations. Continuity is formulated using the particle volume and density [26]: where m a , the mass of particle a, is a constant. This expression of the density does conserve the mass exactly and eliminates any numerical diffusion of density at phase interfaces [27]. The conservation of linear momentum of particle a is expressed by: where u is the particle velocity and f a,p , f a,v , f a,st and f a,g are the forces due to pressure, viscosity, surface tension and gravity, respectively. They are expressed as: Equation (5a) exhibits a particular type of pressure gradient where the two pressures are summed. This formulation has the advantage of (i) conserving linear momentum locally and (ii) increasing the stability of the simulation. Its counterpart is a high sensitivity of the gradient accuracy on disordered particle lattices [28]. Note that there exists another type of gradient where the two pressures p a and p b are subtracted instead of being summed. This type of gradient leads to a better approximation on disordered particle lattices but it does not conserve momentum locally and cancel the stabilizing effect of the background pressure, presented later. The formulation of the viscous term (Equation (5b)) is taken from Cleary [29]. It conserves locally the linear and the angular momentum. The interparticle viscosity µ ab is proposed by Szewc et al. [30]: where ν = µ/ρ is the kinematic viscosity. This formulation, which combines an arithmetic average for the kinematic viscosity and a harmonic average for the density, has previously been shown to feature good stability characteristics, especially for predicting airblast atomization where the shear stress of the light phase is very high [31]. Surface tension (Equation (5c)) is adapted from the Continuum Surface Force (CSF) model [32] following Hu and Adams [27] and Adami et al. [33]. According to the weakly compressible assumption, the pressure p a at particle a is expressed by the equation of state as originally proposed by Cole [34]: where ρ 0 , γ and p back are the nominal particle density, the polytropic ratio and the background pressure, respectively. Note that the background pressure does not appear in the original formulation Cole [34] but was later added by Colagrossi and Landrini [35] for stability purpose. The artificial speed of sound (c) must satisfy the condition c 10 u max for the weakly compressible assumption, ensuring the density variation to be lower than 1% [36]. Furthermore, it was observed in practice that with the present multiphase formulation and equation of state, two additional conditions need to be imposed. First, it is recommended to have: in order to improve the numerical stability at the interface between medium 1 and medium 2. When Equation (8) is not satisfied, the computation of the pressure gradient at the phase interface leads to an exaggerated impact of the particle distribution from one medium on the other one. More details can be found in [31]. Second, it is necessary to limit the influence of the background pressure versus the effect of compressibility. This is achieved by keeping the ratios: as low as possible, where the subscript i refers to the phase. It was observed in practice that a value between 0.01 and 0.2 for Γ i is sufficient. For a too low value of Γ i , the simulation becomes numerically unstable whereas Γ i 1 leads to artificial particle acceleration due to the background pressure which is much larger than the acceleration due to the physics (pressure gradient, viscosity, etc).
Equation (4) is integrated in time by a predictor-corrector scheme. The time step of the numerical simulation is controlled by convection, viscosity, surface tension and the overall acceleration. It is computed at each iteration step: where Γ a is the magnitude of the particle acceleration and C 1 to C 4 are constants equal to 0.5, 0.25, 0.25 and 0.5, respectively.

Geometry
The geometry studied in the present work is illustrated in Figure 2. It consists of a circular plate (the substrate) surrounded by a layer of air. In order to reduce the computational costs, the air layer thickness H out is limited according to the initial droplet position, and an inlet tube is added for injecting the droplet. At the inlet of this tube a constant air velocity is imposed to represent the air entrained by the droplet. The bulk velocity is set to 10% of the droplet velocity. An outflow condition is imposed on the side of the domain, and a no-slip condition is imposed on the impact plate. At the top of the air layer, a wall with slip boundary condition is set as boundary condition. This choice, i.e., a wall and not an inlet or an outlet is motivated by the fact that the flow above the plate will be mostly radial, thus parallel to the surface of the boundary, which typically leads to numerical instabilities in the case of outlet/inlet boundary conditions. The dimension of the geometries for the validation and real condition cases are reported in Table 1.  Figure 2. Sketch of the geometry for the validation case.

Initial Conditions
The initial conditions of a SPH simulation can greatly affect the results [37]. This was also observed here as explained in the following. In case of the simplest initialization procedure, the particles are placed on a regular grid. Here, the particles were placed on a cylindrical grid, which describes concentric circles separated by ∆x in the radial and azimuthal directions. With such an initial condition, preliminary tests showed that the impacts of viscous liquid (µ l ≈ 0.1 Pa·s) led to splash in an hexagonal shape, which was not confirmed by any experiment. The discrepancy was attributed to the initial particle location. Therefore, a specific procedure was developed and applied to overcome this numerical artifact. In a first simulation, the domain was filled in with air particles only and flushed with a high velocity from the inlet tube. Second, the inlet velocity was set to zero, and the simulation was run for a few time steps to relax the pressure gradients due to the high velocity flushing. After this step, the particles inside the domain are not aligned on a regular lattice but randomly packed. By this initial particle arrangement, the unwanted effect of the pressure term according to the equation of state (Equation (7)) could by minimized. Finally, the particles located at the prescribed droplet initial locations were given the appropriate liquid properties and the initial droplet velocity.

Validation with a Liquid Droplet
In this section the numerical method is validated against a reference experiment [38] with viscous liquids to assess the fidelity of the numerical model for low Reynolds number impacts. Note that even with large viscosity liquids, the Reynolds number differs by one order of magnitude between the validation experiment (Re ≈ 10) and the real conditions (Re ≈ 1). However, these conditions always represent the same regime of laminar flow. Moreover, as it will be confirmed later, the very high viscosity of molten sand is expected to strongly limit the deformation of the sand particle. It is thus expected that the particle will not splash nor be disrupted into ligaments and that the spread will be small. Therefore, we focus in this part of the study on the early phases of the impact where the liquid spreads on the substrate.

Phenomenology and Scaling of the Early Phase
The early phases of the impact on a dry cold surface are described in the literature as follows. When the distance between the droplet and the surface tends to zero, the flow in between can be considered as a lubrication flow, and the pressure will rise locally, leading to the deformation of the droplet. Air can even be entrapped, generating tiny bubbles inside the liquid [39]. Then the droplet continues to deform as it spills on the substrate. In this phase, the evolution of the droplet deformation and the surface covered by the liquid is mainly driven by the initial diameter, the velocity and the shape of the droplet. Thus, it is called the kinematic phase where the surface covered by the liquid is roughly equal to the surface given by cutting a sphere by a plane at the location at the wall. If the kinetic energy is not fully dissipated during this phase, the liquid will start to be spilled over the substrate. In this case, the phenomenon is controlled not only by kinetic energy but also by viscosity and the wetting condition of the system (surface tension and surface roughness) which is typically quantified by the contact angle. It is very unlikely that a molten sand particle impacting a turbine blade will reach this phase, which will be corroborated by the numerical simulation.
The main quantity to be investigated in the present study is the liquid surface in contact with the substrate, which is discussed in terms of the spread factor. This corresponds to the nondimensional equivalent diameter of the contact surface: where D 0 is the diameter of the droplet prior to the impact, and S c is the liquid surface in contact with the substrate. In the context of turbine blade impaction by molten sand particle, the contact surface is important because it governs the chemical reactions between CMAS and TBC. Therefore, the spread factor will be investigated in the following. The transient of the early phase is expressed with respect to the nondimensional time t * = U 0 t/D 0 . During the kinematic phase, the spread factor is reported to evolve as D * ∝ (t * ) a with a close to 0.5.

Results
Five test cases were investigated. One was used for qualitative comparison and four for quantitative validation. These are summarized in Table 2. The cases contains approximately 4.7 millions particles and they were run on 125 cores for 150,000 time steps, amounting to 12 ms of physical time. For Case 1 to Case 4, the initiation of the solution was generated according to the procedure described in Section 2.3. The time series depicted in Figure 3 shows side by side images from the experiments in [38] (top) and the corresponding simulations (bottom). Case 0 was selected for a first qualitative comparison. The droplet surface of the numerical simulation was extracted using the α-shape algorithm [40]. An extensive description of our workflow for postprocessing was presented previously [31]. The numerical simulation is found to be in good agreement with the experiment in terms of overall shape and spilled surface. However, it does not correctly capture the rim forming at the front of the lamella, which is attributed to the too coarse spatial resolution of the simulations. In the scope of our study, as the lamella is very unlikely to form during the impact, no rim is to be expected. Therefore, the absence of this structure from the simulation is not detrimental for the validity of the results. One can observe some dispersed lumps of liquid at the periphery of the main region of fluid, especially at the two last time steps (t = 1.31 and 2.27 ms). These are due to the slight numerical instabilities arising on the phase interface during the early phase of the impact. Due to the rapid deformation of the particle lattice in the vicinity of the interface, some liquid particles are ejected outwards and eventually impact the substrate at larger radii. As it represents a negligible volume compared to the total liquid volume, this effect can be neglected in the evaluation of the spread factor.
The pressure was averaged over the symmetry axis of the droplet. It is shown in Figure 4 for early phases of the impact. For t * < 0.1, a checkerboard pattern is observed, which then relaxes and disappears at t * < 0.3. This pattern is typical for a weakly-compressible SPH simulation initialized with regularly distributed particles. Note that only Case 0 was initialized in this manner. The same field was checked on Case 1 (not depicted here) which was initialized with the procedure described in Section 2.3. In this case, the checkerboard pattern was not visible, which corroborates the fact that the checkerboard pattern is due to the initialization. However, vertical stripes were visible up to t * = 0.3, highlighting spurious pressure fluctuations which are typical with weakly-compressible SPH. As the pressure field shows oscillations in the early phase, the overall impact pressure will be then computed in terms of momentum time variation, as detailed later.   For Case 1 to Case 4, the transient of the spread factor was monitored and compared to the experiment. To monitor the transient, it is first necessary to determine the time of impact, which may be a large source of uncertainty in experiments [38]. To extract this incident from the simulations, we monitored the time step where at least one SPH particle is located closer than 0.75 dx to the wall. The spread factor D * was extracted by two means, illustrated in Figure 5. In the first method, the coordinates of the liquid volume in contact with the substrate are expressed in the local cylindrical coordinate system centered at the impact location, and the volume is cut into 360 sectors. In each sector, the particle of maximum radius R i is selected. The spread factor is computed as the surface of a disk whose radius is equal to the ensemble average of all R i . In the second method, the contact surface is defined as the volume of liquid V eq between z wall and z wall + ∆z divided by ∆z. In the present case, ∆z is set to the inter particle distance dx. The first method is valid only for a normal impact leading to a symmetrical spreading, whereas the second method can be used for any kind of impact resulting in a lamella of an arbitrary shape. The spreading computed with the two methods is compared to the experiment in Figure 6. It corresponds to the kinematic and the spreading phases. For all four cases, the spreading computed with method 1 is in good agreement with the experiment. The results given by the second method show a mean deviation of less than 6% for all cases except for Case 3 where it is 10%. We consider this deviation as acceptable, and therefore we will use the second method for the the study of realistic impact scenarios. The close agreement of the simulated and experimental results shown in Figure 6 indicate that both the kinematic and the spreading phases are correctly predicted. While the former was not a major challenge because it just required proper implementation of conservation of linear momentum and volume, the latter is a promising achievement. It demonstrates that the relative contributions of surface tension, viscosity and wall shear stress are correctly predicted.

Application to Realistic Operating Conditions
In this section, the numerical model is applied to real conditions of molten sand particles impacting gas turbine blades. The parameters of the reference case are given in Table 3. Sand particles can exhibit a smooth surface, edges with sharp angles and even corners. Therefore, spherical particles and cubic particles with different orientations to the substrate were investigated. We consider here molten sand particles with a diameter of 25 µm impacting on a solid and nondeformable wall with a velocity of 250 m/s. The mechanical properties are similar to the ones in [12]. The density is 2665 kg/m 3 . The sand particle temperature is estimated to be 1300 • C, leading to a viscosity of 11 Pa·s. The surface tension is set to 0.4 kg/s 2 . Air properties were taken at a temperature of 1300 • C and 30 bar, leading to a density of 10 kg/m 3 and a viscosity of 5.5 × 10 −5 Pa·s.The large viscosity and large velocity represent an atypical configuration characterized by a very small Reynolds number and a very large Weber number at the same time. This suggests a hierarchy of physics such as surface tension inertia viscosity. The surface tension is thus expected to play a negligible role, which will be confirmed subsequently.
The change of temperature ∆T of the particle due to the dissipation of kinetic energy can be estimated by considering that only the sand particle is heated, leading to ∆T = U 2 0 /(2c p ) where c p is the heat capacity of sand. In our conditions (c p = 800 J/K/kg), ∆T = 40 K, which can be be significant for the viscosity, as shown in Table 2 from [42]. However, in our parametric study, we investigate a variation up to ±50% of each parameter, including viscosity. Therefore, the maximum decrease of viscosity during the impact is taken into account in the parametric study. Concerning the practical conditions of stability (Equations (8) and (9)), the values selected for the speed of sound in air and in molten sand are 5000 and 2500 m/s, respectively. This leads to a ratio R = ρ g c 2 g /ρ ms c 2 ms of 0.015, Γ g = 0.8 and Γ ms = 0.012. Although R should be close to one, no numerical instability was observed in the simulation. This is attributed to the large viscosity of the molten sand, leading to a slow deformation of the particle lattice and hence a not too strong source term in the computation of the pressure gradient.
The dimensions of the domain are given in Table 1 and the simulations were initialized as explained in Section 2.3. The domain contains approximately 1.3 million particles and was run over 125 cores for a wall time of 2 h, achieving 30,000 time steps. This leads to a computational cost of 250 CPUhour to simulate 1 microsecond of physical time. The molten sand particle is composed of ≈65,000 particles. Because non-normal impacts were also investigated, the spreading may be nonaxisymmetrical and thus, it was estimated based on the volume in contact with the substrate (method 2).
The outcome of impact is illustrated in Figure 7 for a spherical particle and for a cubic particle impacting the surface at three different orientations. The particle are colored by the height (i.e., the z-coordinate) in yellow-green shades, and the particles in contact with the wall are marked in white. First, no spreading of the lamella is observed, and the particles conserve most of their initial shape. Second, when the particle is cubic, its final deformation depends on the orientation during the impact on the wall. Face impact leads to a very small deformation, whereas corner impact will lead to a strong deformation of the bottom part of the particle. The quantitative results to be presented in the following will demonstrate that the contact surface at the end of the impact is weakly dependent on the particle orientation.

At impact
After impact

Influence of Impact Conditions of a Spherical Particle during Normal Impact
First, the influence of mechanical properties as well as the magnitude of the velocity are investigated for the case of a normal impact of a spherical particle. The dimensional parameters were varied to −50%, −20%, 20% and 50% of the reference values as listed in Table 3, leading to a total of 21 different simulations. In the following, different physical quantities will be extracted and plotted versus time. The nondimensional spreading diameter D * is computed from the volume V c in contact with the substrate as given by Equation (11) where S c = V c /dx is the surface of the particle in contact with the substrate. Two different nondimensional times are used: where U N,0 is the normal component of the velocity. The nondimensional total force exerted on the substrate by the particle is computed from the nondimensional time derivative of the nondimensional particle momentum q * such as: where the subscript i refers to the normal (N) or the tangential (T) component of the velocity. The term q * i is expressed as the summation of the momentum of each numerical particle normalized by the initial momentum component: where the summation (subscript b) is performed on all the numerical particles that make up the sand grain. Since no splashing was observed, the mass m 0 of the sand grain can be considered constant during the impact. In addition, each numerical particle of molten sand has exactly the same mass, thus m b can be taken out of the summation and m 0 = N p m b , N p being the number of numerical particle composing the sand grain. It follows that: where the operator · refers to the ensemble average over the numerical particles of the sand grain. Substituting Equation (15) in Equation (13) and expressing the nondimensional time by its dimensional counterpart leads to: The appearance of two different scaling velocities in Equation (16) is justified by the fact that one stems from the nondimensional time t * where the total velocity U 0 is used and the other stems from the nondimensional velocity, where the consideration of the velocity component (normal or tangential) is important. The necessity of this choice will be justified when discussing oblique impacts in Section 4.4. The nondimensional contact pressure is estimated in the same manner by dividing the force by the contact surface: The computed influence of the sand particle diameter on the evolution of D * , F * and p * is shown in Figure 8 Left. First we describe the general shape of the curves. D * increases approximately with ∼ √ t * because the impact ends in the middle of the kinematic phase as the kinetic energy is completely dissipated by the large viscosity of the molten sand. The normalized force exerted on the substrate shows a maximum at t * ≈ 0.15, and the maximum normalized pressure is obtained at t * < 0.1. Globally, the particle diameter has a weak influence on the normalized quantities, except for a slight decrease of the final spread with D 0 and larger and earlier peak of p * with larger D 0 . The normal velocity (Figure 8, right) has a stronger effect on the normalized quantities. The influence of the viscosity and density is shown in Figure 9. The global shape of the curves is kept similar except for ρ = 1319 kg/m 3 (Figure 9, right) where D * decreases and F * N becomes negative at the same time, indicating the first phase of a rebound. However, the sand particle does not detach from the surface, as D * does not decrease to zero. Finally, the variation of the surface tension from −50% to +50% of the reference value led to exactly the same transients. This confirms the very weak influence of the surface tension, and thus this transient is not shown.

Correlation Characterizing the Normal Impact
The smooth and monotonic variation of the curves from Figures 8 and 9 suggest simple correlations between the output values of interest (D * , F * N and p * N ) and the physical parameters (D 0 , U 0 , µ l and ρ l ). However, we are not interested in finding correlations for the whole transient signal but we focus on characteristic values of the transient. Therefore, we define scalar quantities for each output value that we correlate to the input parameters. For D * , we define the final spread factor D * and T 95%,D * the time where the spread factor reaches 95% of its final value. For F * N and p * N , we define its peak value φ * max (where φ = F N or p N ), the time T max,φ when the peak value is obtained, the time T 5%,φ when the value decreases to 5% of the peak and φ * N,mean the time average between t * = 0 and T 5%,φ . For each quantity, we fit the coefficients (a, b, c, d, e) of the function: where the normalizing parameters and the results of the fitting are given in Table 4. To visualize the accuracy of these correlations, we plot the results of Equation (18) versus the values extracted from the simulations in Figure 10, where the black plain line represents y = x and the grey dashed lines represent a variation of 5%. For most quantities the correlations are within the ±5% limits of accuracy, except for F * N,mean , T 5%,F N and T max,p N . In the case of F * N,mean , it seems that an additive constant could lead to a better fitting while in the case of T max,p N it seems rather constant around 0.07. Apart from these three, the correlations (Equation (18) with Table 4) could be embedded in large scale simulations. Table 4. Correlations on peak values. Variables a to e refer to those in Equation (18). Reference values are expressed in SI units in the formulae.  Table 4 with the values extracted from the numerical simulation.

Influence of the Particle Shape
The impact of a cubic particle was investigated for different orientations towards the wall, thus leading to a study of the impact of different geometrical features summarized in Figure 7.
In Figure 11 Left we compare when the cubic particle impacts the wall with a face, an edge and a corner, to a sphere of the same volume. The first striking element in Figure 11 Left is that the final spread D * is almost independent of the type of geometrical feature impacting on the wall. Concerning the transients, it is found that the spread increases faster when the initial impact surface increases. The extreme cases are on one hand the cube corner where ideally the impact starts by a infinitesimal point touching the wall, and on the other hand the cube face where the initial spread is close to the final spread. Concerning the force exerted onto the wall, the shape of the curves is similar, and only the values of the peak and its temporal evolution depend on the geometrical features. This is due to the deformation of the particle during the impact. With a low amount of matter in the direction of the impact (e.g., the corner impact), the particle deformation is larger and occurs over a longer time, thus leading to lower strain and hence a lower force. To the contrary with a face impact, the amount of matter relative to the projected area is maximal (square impact surface vs. the same projected square area), so that the deformation is minimum, thus leading to a large impact force. The curve of the pressure p * N has a similar shape except for the corner impact where p * N peaks in the very early phase of the impact. This peak can be explained by the very small surface of the corner, which results in a high surface pressure of the predictions. Thus, in order to avoid a too large numerical uncertainty, the contact surface was taken into account in Equation (17) only when it contains more than three SPH particles. The presence of the peak of p * N for the corner corroborates the findings from [3,11] stating that sharp particles are more deleterious for the blade than smoother ones. Indeed for the same particle deceleration, the local pressure on the substrate is much larger for sharp particles which will lead to an augmented erosion. For the three other geometrical features, the peak of the pressure appears sooner than the peak of force, as for the results of Section 4.1. The transients for the different cube orientations considered are shown in Figure 11 Right. The orientations are identified by the rotation angles (in degree) of the cube around the axes defining the plane of the substrate passing by the cube center. With this nomenclature, (0, 0), (45, 0) and (45, 45) represent an impact with a face, an edge and a corner. The same trends as in Figure 11 Left are visible: the final D * is roughly similar with a mean value around 0.811, and the evolution of the transient is faster when the matter contained in the direction of the impact is larger. The shape of normal force F * N is also regular with different peak locations in (t * , F * N ). For the pressure, a peak at the early phase appears for small initial contact surfaces such as the cases (45, 45), (30,45) and (30,30).

Oblique Impacts for Spherical and Cubic Particles
The influence of the angle of impact θ i is investigated in this section. For the spherical particle, angles of 0, 15, 30, 45, 60 and 75 • were investigated while the cubic particle impacts were studied for angles of 0, 45, 60 and 75 • , with different orientations, as sketched in Figure 12. a b c d Figure 12. Illustration of the impact of the cube presenting its face (a), edge perpendicular (b) and parallel (c) to the tangential velocity and its corner (d). The arrows represent the tangential velocity.
For oblique impacts the tangential component of the velocity is nonzero. Therefore, the tangential component of the force and pressure exerted on the substrate can be studied. In the following we will drop the investigation of the force during impact and focus on the pressure. The tangential component of the pressure p * T is referred to as the shear stress τ * . It is recalled that t * (resp. t * N ) is normalized by the total velocity U 0 (resp. normal velocity U N,0 ) whereas normal (resp. tangential) forces and pressures are normalized by the normal (resp. tangential) component of the velocity.
In case of an oblique impact, the particle may slip along the substrate, until the kinetic energy corresponding to the tangential velocity component is dissipated by the friction at the interface. Thus, the final location of the particle is different from the impact location. In the following, the distance L between both locations is referred to as the slip and it is labeled in its nondimensional form as: The arbitrary representation of the nondimensional slip L * in Equation (19) is motivated by the results to be presented in the following. The transients of the oblique impact for the sphere are shown in Figure 13. It is found that the curves of D * and L * match well when they are plotted versus t * N instead of t * . The final value of D * is ≈ 0.79 for θ i ≤ 60 • while the final slip is ≈ 0.6 for θ i ≤ 45 • and ≈ 0.48 for θ i ≥ 60 • . Concerning the normal pressure, the scaling given by Equation (17) provides a satisfactory constant peak value of ≈ 11.6 at t * ≈ 0.08. The negative pressure for θ i = 75 • indicates the first stage of rebound, even though the particle stays always in contact with the substrate. The peak of the shear stress τ * decreases monotonically from six to three and occurs at various times between t * = 0.06 and 0.15. For angles between 15 • and 75 • , the following correlation for the peak of τ * is proposed: with θ in degree.  For the edge impacts (Figures 15 and 16), the curve of D * suggests that the cube is more likely to roll over when the edge is aligned perpendicular to the velocity. This can be explained geometrically by the fact that at the early phase of the impact, the edge in contact with the wall acts as a tipping line for the cube. The curves of the normal pressure and the shear rate collapse on each other for the perpendicular edge impact (Figure 15), confirming the appropriate scaling, whereas they are slightly stretched in time for the parallel edge impact (Figure 16). In both types of impact, the peak of p * N and τ * are between 12 and 15, and between 3.8 and 4.9, respectively. The curve of sand particle slip L * shows a roughly similar increase for all cases of both types of impact up to t * N = 0.5. Afterwards, the final values are confined to a regime between 0.4 and 0.7. The corner impact was investigated for different oblique impact angles. Results of the transients are shown in Figure 17, where a roll over is observed for θ i = 75 • . The final spread and slip are approximately independent of the impact angle with D * end ≈ 0.85 and L * end ≈ 0.60. The evolution of P * N,max and P * T,max is not monotonic with θ i , thus no correlations are proposed. It was seen that for large impact angle, the particle is suspected to roll on the substrate. Time series of cubic particles for an oblique impact of θ i = 75 • is presented in Figure 18 where the particle are colored by their unique identity tag. The time step between two snapshots is approximately 0.14 µs and the last snapshot corresponds to the end of the simulation when the particle has stopped rolling (≈1.5 µs after impact). Globally, as for a normal impact, particle deformation is larger when the amount of material touching the substrate is lower. Hence, a face impact (e.g., Figure 12a) weakly deforms the particle, whereas a corner impact (e.g., Figure 12d) creates a tail-like protuberance on the particle. The difference in the deformation has an influence on the number of particles deposited on the surface. When the particle is more deformed during the impact, it is more prone to deposit material on the substrate. Please note that since no microscopic interaction (e.g., chemical bounds, adsorption) between the sand and the substrate are modeled, the number of numerical particles deposited is not sufficient to quantify the mass deposition on the surface. In all cases, the particle has rolled over half a turn, which corresponds to three faces in contact with the substrate, i.e., approximately 50% of particle surface has touched the substrate at the end of the impact.

Conclusions
In this paper the high velocity impact of a molten sand particle (approximated as a very viscous liquid) of different shapes on a plane, solid substrate was investigated by means of the SPH method. The method was validated against experiments with impacts of water, glycerin and silicone oil droplets. By analyzing nondimensional quantities of real operating-condition impacts, the viscosity and the inertia of the sand particles were shown to be the governing physical quantities driving the impact scenario. Due to the high viscosity of molten sand, no spreading of lamella were observed; however, partial deformation of the molten particles was observed.
A series of simulations were performed varying particle mechanical properties as well as the impact velocity and the particle diameter by −50%, −20%, 20% and 50% of the reference value. From these results, correlations were derived which characterize the spread factor, the distance to impact and the maximum forces and pressures depending on many operating parameters and mechanical properties of the molten sand. Moreover, scaling laws were proposed in order to increase the applicability of the correlations.
Investigation of the effect of the particle shape under normal impact showed that only the transients are affected by the shape but not the final spread. Concerning the influence of the impact angle θ i , for θ i ≥ 75 • , all nonspherical particles were found to roll over on the substrate surface before reaching the final position at higher impact angles.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: