On the derivation of SPH schemes for shocks through inhomogeneous media

Smoothed Particle Hydrodynamics (SPH) is typically used for the simulation of shock propagation through solid media, commonly observed during hypervelocity impacts. Although schemes for impacts into monolithic structures have been studied using SPH, problems occur when multimaterial structures are considered. This study begins from a variational framework and builds schemes for multiphase compressible problems, coming from different density estimates. Algorithmic details are discussed and results are compared upon three one-dimensional Riemann problems of known behavior.


INTRODUCTION
The Smoothed Particle Hydrodynamics (SPH) numerical method is widely used for the solution of a broad range of shock propagation problems [20]. SPH has been extensively applied to the simulation of hypervelocity impacts of solid material objects [12,13,16], which induce shock waves into the colliding objects and are characterized by impact velocities higher than the sound speed of the involved materials-typically 10 km/s. For hypervelocity impacts into inhomogeneous materials, it is commonly noted that SPH schemes produce large errors when shocks propagate through material interfaces [5,15,25,26]. Most hypervelocity impact studies with SPH opt for a homogenization of the inhomogeneous structure [7,6,26,30] and thus neglect reflections and transmissions which occur whenever a shock encounters a material interface [9,31]. Special treatments have been proposed [15,25,26], which are not within the SPH variational framework described in [20,24]. In order to find a treatment coherent to this SPH framework, one needs to recognize that shock propagation through inhomogeneous materials is effectively a multiphase problem; discontinuities appear in density and material parameters.
SPH solution of multiphase problems in the weakly compressible regime-i.e. where density varies maximally 1 % of its initial value-has received noteworthy attention [8,14,11,22], particularly due to the inherent abilities of SPH to perfectly describe advection perfectly and to keep track of interfaces between fluids. Although critisized [22], the number density estimate of [14] at some occasions performed better than standard formulations [8,22]. Therefore, two questions arise: whether shock propagation through inhomogeneous solids is possible within this general SPH framework; and subsequently which scheme is the most suitable.
Aiming at answering the questions above, the current article describes the development of three fully compressible SPH schemes starting from the general SPH framework of [20,23] and working this out for three different density estimates. Additionally, the correct treatment of steep density gradients and discontinuous material properties through the use of artificial dissipative terms is discussed. Finally, the behavior of each scheme is analyzed with three one-dimensional tests: an isothermal impact into a material of purely discontinuous density profile, the classical shock-tube benchmark test [21,28,29] and an isentropic equivalent of the first test.

DENSITY ESTIMATES
In SPH, continuous media are discretized into a finite number of particles. Each particle i ∈ [1, N] is assigned with the local properties of the medium it represents, and a particle mass m i := ρ 0 , i V 0 , i . This constant-through-time mass depends on the discretized medium's initial density distribution ρ 0, i and a chosen partitioning V 0, i of the initial physical domain Ω, satisfying Σ j V 0, i ≈ ∫ Ω dV. At any time instance, material density at a particle's position is readily provided by the traditional SPH estimate [21]: (1) or the number density estimate discussed in [14]: (2) while the next estimate-not directly exploitable-proves to be useful further on: The latter is based on time varying particle volumes 〈V i 〉 = m i /〈ρ i 〉 and the discretization of the kernel function's W ij : = W(|x i -X j |, h i ) scaling property: ∫ Ω WdV = 〈1〉. Other properties demand that the kernel is positive, sufficiently smooth and becomes zero at a finite distance from x i (i.e. compact support) [20]. For unbounded domains, the Gaussian function is an intuitive example, since it becomes practically zero at some distance from x i . The parameter which characterizes the radius over which the kernel decays is the smoothing length h i , and in the weakly compressible regime it is kept constant in space and time. However, in the fully compressible regime the variation of the kernel function with respect to the local particle concentration cannot be neglected. A way to introduce adaptivity of resolution in n-dimensional problems is to postulate: (4) with η = 1.2 the reported optimal value [20,23].
Equations (1) and (2) have to be solved simultaneously with (4) and since they discretize integrals, they are integral forms of mass conservation [20,24]. Although they are effective at resolving shocks [23], in case physical boundaries exist in the computational domain the assumption of ∫ Ω WdV = 〈1〉 is inaccurate and a flat density profile cannot be estimated accurately. Thus they are not popular in SPH simulations of hypervelocity impacts. A simple way to improve the density estimate on boundaries is to use differential forms. Temporal differentiation of (1) and (2) is quite straightforward, while considering 〈V i 〉 = m i /〈ρ i 〉, (3) becomes: Under the previous assumptions on W ij , the last RHS term becomes: (6) approximately equal to the first RHS term. This extra assumption is probably an explanation for the worse performance of the emerging SPH scheme at the tests studied.
With the sources of errors for each density estimate in mind, we drop the brackets denoting approximation. Thus the differential forms of estimates (1), (2) and (3) respectively write: (7) and only involve the temporal variation of the kernel function: (8) It is crucial to underline that the latter depends on particle distance (∇-terms denote explicit spatial differentiation) and smoothing length, as well.

VARIATIONALLY CONSISTENT SPH SCHEMES
The Lagrangian function is defined as the difference between kinetic and internal energy per unit mass e, which is assumed to be a function of two other thermodynamic quantities. Density and entropy per unit mass s are useful choices. The Lagrangian is: (9) and considering the medium's introduced volume partition V i = m i /p i its discrete form becomes: (10) Likewise a typical variational problem [2], the Euler-Lagrange equation leads to the equation of motion for each particle: and its solution provides the particle positions xt and particle velocities v i := ẋ i at any time t: (12) A relation for the internal energy is provided by the first law of thermodynamics de = Tds -Pdv, for temperature T, specific volume v and specific entropy s. The law's particle-wise application actually means that each particle is a small thermodynamic system in equilibrium. For ideal processes-reversible and adiabatic-there is no change in entropy and de = P/ρ 2 dp. The spatial variation of density appearing in (12) depends on the chosen density estimate. Similarly to the temporal differentiation in (7), the spatial variation of density solely includes the spatial variation of the kernel and correspondingly for estimates (1), (2), (3) is: These density estimates lead to the following SPH fully compressible schemes, which reduce to weakly compressible schemes by keeping h fixed and setting any related variations to zero: 2. n-scheme The evolution for each particle's specific internal energy complements the derived fully compressible SPH schemes: (22) It is a particle-wise application of the first law of thermodynamics and practically states that each particle is an adiabatic thermodynamic system capable of performing reversible processes. The equation of state P = P(ρ, e) is a modeling choice, which defines the state of the system and closes each of the schemes above.
Global mass conservation is guaranteed, since particle masses remain constant, while total momentum J = Σ i m i |v i | and total energy E tot = Σ i m i (v 2 i /2 + e i ) are preserved due to the Lagrangian structure of the equations, as proven in [20,24].

COMPARISON TO TRADITIONAL SPH SCHEMES
Variational consistency of equations is crucial in SPH and arbitrary combinations of mass and momentum discretizations introduce numerical noise in the results [3,8]. It has been demonstrated that variational consistency can be more important than devising higher order SPH approximation schemes [24].
Apart from the variational consistency, which is rarely satisfied by traditional SPH schemes for shocks through solids [12,18,16,26], traditional schemes differ substantially from the above derived schemes. Traditional schemes neglect terms coming from the W-h-ρ coupling 1 (Ω-terms in present notation or sometimes called ∇h-terms) and their equation of motion (considering only normal stresses) writes: (23) with . For the update of the smoothing length, the time evolution of (4) is used: . Comparing (23) to (17), note that in the traditional schemes, the contribution of each particle j to the force acting on particle i is averaged, instead of depending per se on the two local discretization scales h i and h j . The observations above have two immediate effects. Firstly, for the traditional schemes, their well-known problems of momentum conservation [16,12] are expected to worsen with increasing shock effects. In this case, hi and hj discretization scales differ substantially and Ω diverges from 1, making the result of (23) significantly different than (17). Secondly, for the m-scheme and particles of equal masses or the n-scheme in any particle setup, the present formulation actually allows for the dependence of the estimates on local particle concentration: u i = Σ j W ij . Consequently, the number of neighbors of each particle is expected to remain constant in time [24]. The latter may play a role in the treatment of numerical instabilities, since particle clumping relates to changes in the number of neighbors and is an exhibition of numerical instabilities in SPH [20,24]. Suggested treatments often involve ways to keep the number of a particle's neighbors constant [1] or a form of kernel adaptivity [27]. Therefore, the present formulation is anticipated to deliver beneficial effects on numerical instability issues.

ARTIFICIAL DISSIPATION
For problems involving shock propagation, special treatment is needed in order to dampen out spurious effects near shocked regions. To this end, artificial diffusion is added. The terms suggested by Monaghan [19] are heuristically shown to have a structure similar to the structure of approximate Riemann solvers. Price [24] discusses their difference to traditional SPH artificial viscosity and notes that artificial diffusion terms compensate for errors caused by the subtle assumption of a differentiable Lagrangian (10) function [23,24]. A generic form of artificial diffusion for each variable ϕ = {ρ, v, e*}, where is [20,24]: with a signal velocity v sig , ϕ and parameters a ϕ and β ϕ , all related to each conserved variable and particle sound speeds

Dissipative mass flux
In the following impact scenarios, the addition of artificial viscosity alone is not sufficient to remove nonphysical spikes that appear in the pressure distribution, when a shock encounters the interface of two materials. They are due to the assumption of differentiable density in (8). We, therefore, devised an artificial diffusive mass flux term, inspired by the structure of (24). However, since density discontinuities are natural characteristics of inhomogeneous materials, smoothing them out is ruled out. An expression in the same manner as (24) needs to account for pressure differences. From dimensional analysis point of view, changing density differences into pressure differences in the general form of diffusive terms (24) needs to be complemented with a division by a characteristic velocity vs;g;P. Thus, we suggest the following structure: From numerical experiments, the optimal value appears to be a p = 0.3, while β p = 0, since it has negligible effect to the results.

Viscosity
Similarly to [24], diffusion in the momentum equation is added in the form of (24): and is restricted to approaching particles (v ij x ij ≤ 0), as: with a v = 1:0 and β v = 2:0.

Thermal conductivity
For problems involving thermal effects, artificial thermal conductivity is added to the evolution of specific internal energy [20,24] as: (27) The above form refers to the physically appropriate conservation of the total energy e*, rather than internal energy a e . The parameter ae is found to be problem dependent, while β ρ = 0, due to its small effect to the results.
In SPH literature, a similar approach of adding diffusive terms in all evolution equations appears in [17,25]; it is termed conservative smoothing. The most note worthy differences are the following: 1) the present dissipative terms belong to a generic form and are designed such that they complement schemes coming from a variational SPH framework [20,23,24], 2) they incorporate the W-h-ρ coupling, 3) they are directly related to approximate Riemann solutions [19] and finally the new dissipative mass flux term is precisely designed to treat spurious spikes on the internal interface of inhomogeneous materials.

TIME INTEGRATION
Monaghan [20] notes that properties of the Lagrangian formulation are better preserved using symplectic integrators. Accordingly, in all following tests, time integration from timestep k to timestep k + 1 is achieved with a leap-frog scheme: A constant timestep Δt = 10 -4 τ, was used without the occurrence of instabilities, where τ is the problem's timescale. In one-dimensional tests, they typically show up when particles overtake each other near shocks.

TESTS AND DISCUSSION
Usually, particles of equal masses are used for shock problems [21,23,4]. Choosing Δx 0,  initial interparticle distances for low-density (ρ  ) regions, particles of masses m: = Δx 0 ,  ρ 0 , i are placed therein and particles of the same masses are placed at Δ 0 , h = m/ρ 0 , h initial interparticle distances in the high-density (ρ 0 , h ) regions. Alternatively, particles of unequal masses m  := ρ 0,  , Δx 0 and m h : = ρ 0, h Δx 0 may be placed at equal initial interparticle distances Δx 0 in the whole domain. The latter distances correspond to equal one-dimensional volumes V 0 = Δx 0 . Suffice to say that in two or three dimensions and arbitrary geometry, a centroidal Voronoi tesselation (for a review see [10]) places particles at their initial positions. In any case, purely discontinuous data are used. Therefore, each of the three derived schemes may be combined with two initial configurations of the particle system. Note that for a system of equal-mass particles, use of either the m-scheme or the n-scheme results to the same equations-for equal masses, m j in (1) can be taken out of the summation. Additionally, the v-scheme with equal-mass particles performed much worse than the m-scheme with the same particle configuration in our preliminary tests and therefore was omitted from the results presented. Four cases remain; the m-scheme is examined for both equal-mass and unequal-mass particles, while the n-scheme and the v-scheme are investigated for unequal mass-particles.
Focus is on how the derived schemes behave in capturing shocks through material interfaces; hence, simulations end long before shocks reach the domain's boundaries. The Gaussian function is used as kernel function, similarly to [21].
The following tests constitute Riemann problems, because they are governed by hyperbolic systems of conservation laws with piecewise constant initial data [29]. We are particularly interested in the pressure distribution through material interfaces, which is expected to be uniform from the Riemann structure of the problems [9,29].

ISOTHERMAL IMPACT INTO INHOMOGENEOUS STRUCTURE (TEST 1)
In the first test, the barotropic equation of state is introduced: The idealized test of an isothermal impact into an inhomogeneous structure consists of a discontinuity in velocity and two contact discontinuities described by: It models the impact of a projectile onto an impactor made of three layers. The first and last layers are of the same material as the projectile, while the medium layer is a material of four times lower density and twice higher sound speed. The impact velocity (projectile's velocity in the current setup) is equal to the sound speed of the projectile's material. At time t = 0, two shock waves of the same pressure are produced at the "impact site" x = 0.6, and start moving in opposite directions. This event is described as a symmetrical plate impact experiment in [9]. When the right-moving shock reaches the interface of the two materials, it splits into a wave that is reflected back to the original material and another that is transmitted through the interface into the adjacent material (unsymmetrical plate impact according to [9]). Two shocks or a shock and a rarefaction may appear, depending on the shock material parameters. Independently of the wave structure, two regions of piecewise constant pressure develop. One at the impact region and another one expanding to left and right of the interface (the socalled star region [29]).
For the equal-mass particle system, the initial interparticle distance in the low density region is Δx 0 ,  = 0.008 and in the high density region it is Δx 0 , h = 0.002, resulting in 550 particles in the problem domain. As a comparison, in [21] and [23], 400 and 450 particles of equal masses are used respectively, for the shock-tube problem in the unit-length domain. For the particle system of unequal masses, 550 particles (400 particles per unit length) are placed at equal initial interparticle distances.
In the upper set of plots in Fig. l, density, pressure and velocity are presented at particle positions with the m-scheme and particles of equal masses at t = 0.40. The solid line indicates the exact Riemann solution (solution details are found in [9]). In the lower set of plots, the same solution components are shown for the m-scheme and unequal-mass particles. A dissipative mass flux, with a ρ = 0.3, is used. We note that In the inset plots of the pressure distribution, the SPH solution in the region of the interface is shown, when no dissipative mass flux is used.
Both schemes capture the correct values of all magnitudes at the impact site and the interface region. The beneficial effect of the dissipative mass flux becomes evident in the pressure plots; the dissipative mass flux suppresses spurious spikes on the contact discontinuity. For the system of equal-mass particles (Fig. l upper plot), pressure distribution is flat, as expected from the Riemann problem. In the solution with the system of unequalmass particles (Fig. l lower plot), the pressure distribution through the interface is not perfectly flat, it has a small hump. In Fig. 2 the same solution components are presented for the v-scheme and n-scheme (respectively upper and lower set of plots) and particles at equal initial spacing, with dissipative mass flux added. Results are similar to the results of the m-scheme with the same initial configuration. Finally, the test is conducted in the same domain, with the same initial values, apart from the Mach number, that writes Ma = {1, 1, 0.5, 1}. The middle layer's sound speed is half the projectile's sound speed and double the impact speed. Results at t = 0.30, using particles of equal masses and the m-scheme, are depicted in Fig. 3 in the upper set of plots. They show that the system of equal-mass particles is unable to represent the flat pressure profile at the contact discontinuity even with the use of dissipative mass flux (Fig. 3 inset pressure plot). Furthermore, an overshoot in pressure and density appears on the tip of the reflected wave. On the other hand, the n-scheme with particles at equal initial distances manages better, with only a small hump on the contact discontinuity ( Fig. 3 lower plots). This effect may pose a significant drawback in the use of equal-mass particle systems for hypervelocity impact simulations.
Regarding total momentum, it is conserved up to machine precision for all schemes. In contrast, by using traditional schemes [12,16,18,25])-set Ω i = 1 in (20) or (16) and use (23) by evolving the smoothing length as dh i /dt =-h i /(nρ i )dρ i /dt-momentum conservation fluctuates at 0.2% of its original value. This accuracy problem seems insignificant (also noted in [12,16]), although it may create serious problems in more complicated impact scenarios or advanced physical models.

92
On the derivation of SPH schemes for shocks through inhomogeneous media Figure 2: Test 1 at t = 0:40. In the upper set of plots, the n-scheme is used with particles at equal intial interparticle distances and in the lower set of plots the v-scheme with the same particle conguration

SHOCK TUBE (TEST 2)
The derived schemes and the new artificial mass flux term are examined against the classical shock tube test [21,29,23]. Initial conditions in the left interval 0 ≤ x L < 0.5 are {ρ L , P L , ν L } = {1.0, 1.0, 0}, and {ρ R , P R , ν R } = {0.1, 0.1, 0} in the right interval 0.5 ≤ x R ≤ 1.0. The ideal gas law, P = [γl)ρe, closes the system, with γ = 5/3 the ratio of specific heat coefficients. A right moving shock and a contact discontinuity are present, along with a left moving rarefaction wave. The performance of SPH at this benchmark test is extensively studied by Price [23], who only focuses on schemes with particles of equal masses and highlights the unphysical spikes on the contact discontinuity, when the differential from of mass conservation is used. Accordingly, the m-scheme is employed with equal-mass particles placed initially at spacings Δx h = 0.001 on the right and Δx  = 0.008 on the left (450 particles in total). Results at t = 0.2 are plotted in Fig. 4 along with the exact solution of the Riemann problem (solid line). The addition of the dissipative mass flux (a ρ = 0.3) and artificial conductivity with a e = 0.5 and β = 0 dampen out the-reported by [23]-spurious spikes in the pressure and density distributions. Although a small hump persists in pressure distribution, the artificial heat term removes a spurious spike in internal energy that occurs near the shock.
It is argued [4] that SPH algorithms involving particles of unequal masses have problems in resolving strong shocks. Indeed, without adding the dissipative mass term in an unequalmass particle system, a particle near the shock is found to move faster than its frontal neighbors and overtakes them. The addition of dissipative mass flux manages to treat this unphysical effect in all schemes (m-scheme, n-scheme and v-scheme) when unequal-mass particles are used. Indicative results using the n-scheme with particles placed at equal initial interparticle distances (400 particles per unit length) are plotted as detail plots in Fig. 4. Although the representation of shocks is steeper than with the m-scheme and equal-mass particles, a larger hump appears on the interface for the density distribution.

ISENTROPIC IMPACT INTO INHOMOGENEOUS STRUCTURE (TEST 3)
The last test comprises of an ideal (isentropic) hypervelocity impact into an inhomogeneous structure using the nonlinear equation of state: Parameters and initial conditions are described by:  and are in the cm-gr-μs system, with all parameters (K 1 , K 2 , K 3 , B 0 and B 1 ) given in Mbar. Values of the projectile's material are reported for aluminium [13], while the values of the impactor's medium layer material are fictitious and relate to a material of density and sound speed lower than the projectile's material. Normal stresses of metals during hypervelocity impacts are typically modeled with this type of equations [13].
No exact solution exists in order to validate the numerical results, however the purpose of this test is to highlight the importance of the dissipative mass flux term and examine whether its interplay with the artificial heat conduction coefficient a e delivers a flat pressure distribution across shocked material interfaces. We are practically interested in the behavior of the n-scheme with particles of equal initial volumes and especially in its ability to suppress any spurious kinks in the distribution of computed variables.
In Fig. 5 results at t = 0.3 μs are shown, using a coarse discretization of 280 particles (200 per unit-length). Dissipative mass flux (a ρ = 0.3) and dissipative heat with a e = 4.0 are used. Note that the substantially increased value of a e = 4.0 is needed in order to suppress a large spurious blip appearing at the impact region in the initial stages of the test.
Similarly to the isothermal impact test, the use of dissipative mass flux smoothens out spurious kinks on the shocked interface and provides continuous and almost uniform pressure distribution through the interface. A flat pressure distribution through the interface is expected due to the Riemann structure of the problem. Furthermore, the dissipative heat term removes the energy which is localized at the impact site. This localization is due to the abrupt change of kinetic energy into internal energy during an impact of velocity as high as the material's speed of sound. The localized internal energy is smoothened out to a length of several smoothing lengths, does not affect the strength of the shock and works independently of the dissipative mass term. Apart from noting that the other two schemes (m-scheme and v-scheme with unequal-mass particles) produced a much larger blip on the contact discontinuity, it is worth focusing on the pressure distribution in the inset plot of Fig. 5, where the m-scheme and equal-mass particles are used with dissipative mass flux. As an effect of the shock, the non-linear equation of state (34) dictates changes to the speeds of sound of the materials and creates a catastrophic situation for the computation, similar to the effect observed in Fig. 3. The original claims for momentum and total energy conservation are true up to a precision of 10 -7 , as exhibited in Fig. 6 (upper and lower set of plots, respectively). Their global variations from initial values J o = Σ i m i v 0 , i and E tot, o = Σ i m i (v 2 i, 0 /2 + e 0; i ) are negligible. The lower accuracy of the energy conservation is due to its first order accurate time integration, in contrast to the second order accurate time integration of the momentum in 28.

CONCLUSIONS
Three fully compressible SPH schemes have been developed from a standard SPH variational framework and three different SPH density estimates; they all embed a wellknown density -smoothing length -kernel coupling. They incorporate differential forms of 96 On the derivation of SPH schemes for shocks through inhomogeneous media mass conservation, usually preferred for computations in domains that include material boundaries.
Overall, from the test involving an isothermal impact into an inhomogeneous medium and for the density and speed-of-sound ratios studied, we conclude that all examined SPH schemes are able to follow the exact Riemann solution of the problem. However, we find that an unequal-mass particle configuration with the n-scheme (based on number density) can handle a wider variety of inhomogeneous media configurations under shock. This scheme is able to fairly follow the solution of the classical shock-tube test, without the occurrence of any spurious effects.
Due to the strength of the examined shocks and the use of the differential form of mass conservation, a dissipative mass flux term needs to be added. The current study, provides such a term. Moreover, coefficients for the artificial dissipation terms for impact tests in the order of Ma = 1.0 are suggested such that no unphysical effects occur and results depict the physics of the discussed problems. These coefficients can later be used for simulations in two or three spatial dimensions.
Thus, particle systems of unequal masses are apt for one-dimensional shock propagation through inhomogeneous media. Especially for higher density ratios, this leads to reduced computational costs. Finally, these early findings remain to be validated for the simulation of hypervelocity impacts in two or three dimensions and with material models that incorporate shear stresses.