An Accurate SPH Scheme for Dynamic Fragmentation modelling

We focus on the use of a meshless numerical method called Smooth Particle Hydrodynamics (SPH), to solve fragmentation issues as Hyper Velocity Impact (HVI) cases. Firstly applied to fluid flow simulations, this method can be extended to the solid dynamics framework. However it suffers from a lack of accuracy when evaluating state variables as the pressure field. And such inaccuracy generally generates non-physical processes (as numerical fragmentation). In the hydrodynamic context, SPH-ALE methods based on Riemann solvers significantly improve this evaluation, but increase the scheme complexity and low-Mach issues are difficult to prevent. We propose an alternative scheme called γ -SPH-ALE, firstly implemented to solve multi-regime barotropic flows, and secondly extended to solid dynamic cases. It relies on the combination of the SPH-ALE formalism and a finite volume stabilizing low-Mach scheme. Its characteristics are detailed and evaluated through a nonlinear stability analysis, highlighting CFL-like conditions on the scheme parameters. Finally, its implementation on several test cases reveals that the proposed scheme actually increases both stability and accuracy, in reduced computation time, with respect to classical solvers.


Introduction
The dynamic fragmentation is a process involved in various events as impacts on aeronautics and aerospace structures (hypervelocity impacts HVI, ice impacts), warhead fragmentation or shielding design. Combining large deformations as well as numerous interface productions, this mechanism is particularly challenging from a numerical point of vu. However, reliable predictive modelling for such events are crucial as efficient decision-support tools. Eulerian solvers are classically preferred to manage such computations but limited regarding interface tracking, accuracy or energy conservation. Lagrangian solvers, well suited for interface tracking, are then a relevant alternative. Among them is the meshless method called Smoothed Particle Hydrodynamics (SPH) and we focus in this paper on new robust SPH schemes.
The SPH is a Lagrangian particle method used to solve partial differential equations (PDE). The computational domain is discretized in a set of interpolation points interpreted as particles interacting between themselves through kernel functions, and carrying material properties. Originally developed by Lucy [1] for astrophysical problems and by Gingold & Monaghan [2][3][4] for hydrodynamic applications (incompressible free surface flows), the method was then adapted by Benz [5 ,6] to the solid mechanics by including material models.
In the hydrodynamics context, Monaghan [4] proposed an explicit SPH formulation based on weakly compressible hypotheses (WCSPH), to treat complex free surface flows. Such method is particularly attractive for its conceptual simplicity, robustness, ability to handle large arbitrary deformations, to not produce diffusive interfaces and for its monolithic nature. Still, SPH schemes suffer from well-known issues questioning their accuracy: lack of interpolation completeness, nonphysical oscillations, and tensile instability. When it comes to the solid dynamic context, such instabilities become particularly challenging by activating nonphysical processes as numerical fragmentation. Various strategies are employed to solve these issues as high order kernels [7][8][9], additional diffusive terms [2,10,11] or corrective pressure terms [12,13].
In order to deal with all these issues in a robust and consistent way, we focus on a different SPH framework introduced by Vila [14] based on Arbitrary Lagrangian Eulerian (ALE) considerations. The idea is to combine and benefit from both Eulerian and Lagrangian formulations through an arbitrary transport velocity field v 0 . Indeed, an Eulerian description is generally preferred to handle large deformations, and a Lagrangian one for interface tracking. Thus a smart choice of v 0 allows to increase the stability and robustness of the scheme (XSPH [3], Particle shifting [15][16][17]). Generally combined with Riemann Solvers [15], such formulation improves the scheme accuracy, but, at the cost of an increased solving complexity (use of a Riemann Solver) and low-Mach limitations [18]. Regarding this late aspect, in the Finite Volume (FV) framework some authors [19,20] proposed a stabilizing Low-Mach scheme. It relies on a corrective velocity term added to the continuity and momentum equations. A non-linear stability analysis provided them stability conditions (CFL-like) insuring robust and accurate computations whichever flow regime. Clayer et al. [21] adapted this Low-Mach scheme to the SPH framework and achieved similar convincing conclusions.
We propose an alternative scheme called γ-SPH-ALE relying on the combination of the SPH-ALE formalism and the FV low-Mach scheme. Its governing equations and the corresponding stability conditions are exposed in section 2. Firstly detailed for a hydrodynamic version, they are extended then to a solid formulation. Finally in section 3, the proposed γ-SPH-ALE scheme is evaluated through several validation cases (both hydrodynamic and solid versions).

γ γ γ γ-SPH-ALE
In this section the proposed γ-SPH-ALE scheme is detailed. The general SPH framework is firstly exposed and secondly applied to the ALE context. The governing equations of the proposed scheme and the corresponding stability conditions are thirdly exposed in the hydrodynamic version. They are then extended to a solid formulation.

SPH basics
The SPH method relies on approximation features allowing to discretize conservation laws on a particle set (x i ,ω i ) i∈P (respectively the particle position and weight, and P the particle domain). These approximations are achieved thanks to kernel functions W(||x i -x j ||,h)=W ij depending on the smoothing length h (radius of its compact support) and the distance between two particles. The kernel gradient can also be evaluated as ∇ ∇ ∇ ∇W ij =grad x W ij allowing to define then the smoothed particle approximation of a function f (1) and of its derivative (2) To ensure the conservative property of the scheme, alternative derivative operators are used in practice (3,4) and combined to symmetric kernels such that W ij =W ji and ∇ ∇ ∇ As introduced above, classical SPH suffers from a lack of interpolation completeness which can be compensated by using renormalized kernels. Vila [8] proposed to replace ∇ ∇ ∇ ∇W ij by A ij =B ij ∇ ∇ ∇ ∇W ij where B is the renormalization matrix (introduced by Johnson et al. [7] and Vila [8]) defined as Note that, as proved in [14], operator D h (3) combined to renormalization strongly approximates the gradient under the condition x/h=O(1) (where x is the particle spacing) which is a crucial point to ensure the scheme consistency (in standard SPH this convergence property depends on the ratio x/h and is not always enforced [22]).

ALE Framework
The idea is now to apply these approximation features to a conservation law in ALE formalism. We consider then the following law Where v 0 ∈ℝ 3 is a regular transport field, Φ∈ℝ 4 the vector of the conserved variables, 3 α the Eulerian flux vectors, S=(S l ) 4 l the source term and L v0 the transport operator associated to v 0 (8) In the particular case of Euler equations we have Providing the following equation set Where ∇ h * is the adjoint of an operator ∇ h approximating the derivative operator ∇.

Governing Equations
The idea is to choose ∇ h * in order to stay conservative and consistent. Several choices are admissible and we propose to work with the following scheme discretized in space With where the variables x, ω, ρ, v and p refer respectively to the position, volume, density, velocity and pressure of the particles. w correspond to the relative ALE velocity defined by w=v-v 0 . We also define p ij =p i +p j , Γ Γ Γ Γ ij corresponds to the stabilizing term coming from the FV low-Mach scheme and π π π π ij corresponds to the artificial viscosity. 〈.〉 defines the usual scalar product on ℝ 3 .
Disregarding Γ Γ Γ Γ ij and π π π π ij , the classical Lagrangian (w=0) SPH formulation with constant masses m i =ω i ρ i is recovered A vectorial version of Monaghan's stabilizing artificial viscosity [2] is proposed (22). It depends on a parameter α (23) to adjust in order to minimize the amount of artificial diffusion induced in the scheme. c 0 is the material sound speed and ρ ij =(ρ i +ρ j )/2 π π π π ij = α ij (v j -v i ) (22) α ij = αc 0 ρ ij ||A ij || (23) The coupling between the SPH version of the FV Low-Mach scheme and the SPH-ALE formalism is achieved by the term Γ Γ Γ Γ ij (24,25) through the mean relative ALE velocity w ij . We set With this ALE formulation particles are moved following the velocity field v 0 (12). Such field is said to be arbitrary in the sense that, as suggested by Vila [14], a smart choice of v 0 could increase both stability and robustness of the calculations without impacting the scheme stability conditions (29-32). Various methods are available in the literature as Monaghan's XSPH method [3] or Particle shifting [15][16][17]. They consist in reorganizing the particle set during calculations (quasi-Lagrangian mode) to prevent instabilities such that tensile instability [12] or the formation of Lagrangian particle structures [17].

Stability Conditions
In order to ensure the proposed scheme is conservative, robust, stable and consistent, stability conditions on the scheme parameters are exhibited thanks to a nonlinear stability analysis. The proof is not detailed in this paper but we refer to the work of Grenier et al. [19], Lavalle et al. [20] and Couderc et al. [23] for a similar process in a FV framework. The conservative property is straightforward considering that A ij = -A ji , w ij = w ji and H ij = -H ji . Regarding the stability property, performing an energy balance provides a CFL-like condition (29) as well as stability intervals for α (30) and γ (31). These conditions ensure the limited behavior of the scheme total energy E (32) where T is the end time of computation.
Note that conditions (29-31) are obtained in practice after some calculations. The theoretical results are not detailed here for sake of clarity. Besides, they are actually too restrictive due to a loss of optimality during their evaluation. Also they ensure the stability of a first order time integration scheme whereas a second order scheme is used in practice (allowing relaxing the time step restriction). Thus, we choose to work with a CFL condition (29) classically used in the SPH literature [3]. It worth noticing that these stability conditions are defined regardless the expression of v 0 revealing its arbitrary property previously introduced. However they depend on a geometrical constant meaning that the stability is ensured while the particle distribution stays regular enough. Having volume bounds (according to (13)), we can actually show that CFL (29) also ensures the existence of density bounds which gives the expected robustness of the scheme. Considering now the energy and density bounds, we can finally achieve the weak consistency of the scheme thanks to a Lax-Wendroff like theorem. We refer to [14] where Vila details the employed technique.

Solid Formulation
The proposed γ-SPH-ALE scheme is then directly extended to handle solid cases. o do so the momentum equation (21) is enhanced by the deviatoric stress tensor σ' contribution. Its rate of change is then classically evaluated thanks to Jaumann's formulation of the Hooke's law. See [24,25] for similar detailed processes.

Validation
Having the stability conditions of the proposed γ-SPH-ALE scheme, we now validate it on two test cases: isentropic shock tube and hypervelocity impact.

Isentropic Shock Tube
The isentropic shock tube is a classical academic test case with a known exact solution. A large amount of configurations and data are available in the literature allowing then to efficiently evaluate the γ-SPH-ALE scheme. We focus here on a simple barotropic shock The fluid is taken as having a referential density ρ 0 =1000 kg/m 3 and a sound speed c 0 =1466 m/s. The tube is taken with a length L x =0.1 m and a width L y =0.01 m. The discontinuity is initially defined at x=0.05 m. All computations are performed with cubic B-splines kernel functions and a fixed smoothing length. The exact solution used as a reference is given by an exact 1D Riemann solver with 2000 cells. SPH computations are performed in 2D with periodic boundary conditions to mimic a 1D case. We present the results given by an ALE SPH HLLC Riemann solver with a MUSCL scheme, and by the γ-SPH-ALE scheme. The initial spacing is set to x=2.10 -4 m, giving 50.000 particles.
The γ-SPH-ALE scheme is used with α=0.02. As depicted on Fig. 1 the discontinuity generates a left rarefaction wave and a right shock wave. The solution is then characterized by three pressure levels (2.1492 8 Pa, 1.0493 8 Pa, 0 Pa). Figure shows the comparison between the exact solution, the ALE SPH Riemann solver and the γ-SPH-ALE scheme with γ=1.0. Both schemes are able to reproduce the pressure jumps. The ALE SPH Riemann solver manages to capture the jumps more efficiently than the γ-SPH-ALE scheme, generating narrower pressure variations. However the ALE SPH Riemann solver also produces significant overshoots the proposed scheme completely damps. Regarding the choice of γ, Fig. 1 results are obtained with the value corresponding to the threshold below which overshoots appear. Fig. 2 shows the results for different γ. We can see that increasing its value stabilizes the scheme by damping the oscillations in a very efficient way, but it also diffuses the shocks (larger jumps). The γ-SPH-ALE scheme is then able to faithfully reproduce the shocks and the use of γ allows to easily control its stability and accuracy without the Riemann solver complexity. Note that the γ-SPH-ALE scheme produces identical results in Eulerian description (v 0 =0).

Hypervelocity Impact test (HVI)
Hypervelocity Impacts are events occurring at velocities around 1 to 10 km/s. They are particularly challenging for space vehicles design when generated by orbiting debris. Regarding the experimental limitations (to reach the impact velocity), providing reliable predictive modelling is then crucial to evaluate the vehicle vulnerability. We focus here on the case of an Aluminium sphere (diameter 9.53 mm) impacting Aluminium plates with different thickness at 6.7 km/s. The experiment is described in [26]. Two cases are reproduced: a-thickness = 0.8 mm, b-thickness = 4.039 mm. The γ-SPH-ALE scheme is implemented in full 3D with an Elastic perfectly Plastic material model as well as a Mie-Grüneisen equation of state. The spalling process is managed by a pressure threshold above which the deviatoric stress is relaxed (only pressure remains). Such simple criterion is used as the idea is to evaluate the γ-SPH-ALE scheme ability to reproduce the cloud of debris.
The constitutive parameters are the following: ρ 0 = 2703 kg/m 3 , Young modulus 76 GPa, Poisson ratio 0.3, Yield stress 260 MPa, Linear Hugoniot Slope coefficient 1.338, Grüneisen gamma 2.0. The numerical results are achieved with 1.5.10 6 SPH particles. Fig. 3 and Fig. 4 show respectively the experimental and numerical results for both configurations described above. The experimental results are achieved thanks to X-ray flashes. We can see that in both configurations the numerical model is able to faithfully reproduce the cloud of debris characteristic features: an ejecta veil on the plate front side, an expending debris cloud "bubble" on the plate rear side and a disc-like structure inside the expending bubble (appearing clearly on case thickness 0.8 mm).  A similar configuration is now studied and corresponds to an oblique HVI. The sphere diameter is 3 mm with a velocity 4050 m/s and the plate thickness is 2 mm. Both elements still are in Aluminium (same constitutive parameters). However the plate is now 32°t ilted with respect to the sphere impact direction. Fig. 5 shows the experimental results at different times achieved by Thiot Ingenierie (introduced in [27]), and Fig. 6 shows the γ-SPH-ALE numerical results full 3D with 6.5.10 6 SPH particles. We can see that the proposed numerical model is again able to faithfully reproduce the cloud of debris shape.

Conclusion
We have proposed a new mesh-less scheme called γ-SPH-ALE based on the combination of both SPH-ALE formulation [14] and FV low-Mach scheme [19,20]. PDE are written in ALE formalism (7) allowing for a simple treatment of the arbitrary velocity field v 0 , ruling the particle motion (12). A stabilizing term extracted from the FV low-Mach scheme is also introduced by means of a velocity corrective term (depending on a parameter γ). Proportional to the pressure gradient it improves the state variables evaluation. An artificial viscosity (depending on a parameter α) is combined with this γ velocity to ensure the global stability of the scheme.
We evaluated the performances of the γ-SPH-ALE scheme (for both hydrodynamic and solid contexts) on two test cases: isentropic shock tube and hypervelocity impact. For the shock tube case the γ-SPH-ALE scheme managed to faithfully reproduce the pressure field by efficiently damping the spurious oscillations (in both Lagrangian and Eulerian descriptions). Regarding the HVI cases, the γ-SPH-ALE scheme managed to reproduce the cloud of debris for both normal and oblique configurations.
Future works will be the enhancement of material models and fracture modelling to handle more complex materials and configurations.