Approximate Riemann solvers for the Godunov SPH (GSPH)

Accepted Manuscript Approximate Riemann solvers for the Godunov SPH (GSPH) Kunal Puri, Prabhu Ramachandran

PII: DOI: Reference:

S0021-9991(14)00242-3 10.1016/j.jcp.2014.03.055 YJCPH 5177

To appear in:

Journal of Computational Physics

Received date: 19 August 2013 Revised date: 24 March 2014 Accepted date: 28 March 2014

Please cite this article in press as: K. Puri, P. Ramachandran, Approximate Riemann solvers for the Godunov SPH (GSPH), Journal of Computational Physics (2014), http://dx.doi.org/10.1016/j.jcp.2014.03.055

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Approximate Riemann Solvers for the Godunov SPH (GSPH) Kunal Puria,∗, Prabhu Ramachandrana a

Department of Aerospace Engineering, Indian Institute of Technology Bombay, Powai, Mumbai 400076

Abstract The Godunov Smoothed Particle Hydrodynamics (GSPH) method is coupled with noniterative, approximate Riemann solvers for solutions to the compressible Euler equations. The use of approximate solvers avoids the expensive solution of the non-linear Riemann problem for every interacting particle pair, as required by GSPH. In addition, we establish an equivalence between the dissipative terms of GSPH and the signal based SPH artificial viscosity, under the restriction of a class of approximate Riemann solvers. This equivalence is used to explain the anomalous “wall heating” experienced by GSPH and we provide some suggestions to overcome it. Numerical tests in one and two dimensions are used to validate the proposed Riemann solvers. A general SPH pairing instability is observed for two-dimensional problems when using unequal mass particles. In general, the Ducowicz Roe’s and HLLC approximate Riemann solvers are found to be suitable replacements for the iterative Riemann solver in the original GSPH scheme. Keywords: GSPH, Approximate Riemann Solvers, Euler Equations

1. Introduction Smoothed Particle Hydrodynamics (SPH ) is a mesh-free particle method for the numerical solution of the equations of continuum physics. It has evolved considerably since the pioneering work of Gingold and Monaghan [1], who used it as a numerical tool for astrophysical problems. SPH has been found suitable for problems in compressible hydrodynamics [2, 3, 4, 5, 6], astrophysics [7, 8, 9] and magneto-hydrodynamics [10, 11, 12, 13, 14]. Despite recent [15, 16] and historical [17, 18, 19] interest in artificial viscosity, SPH is often criticized for relying on an explicit artificial viscosity to stabilize the scheme. The viscosity is typically added as a viscous pressure to the momentum equation to dissipate kinetic energy at shocks. Depending upon the perspective, the form of the viscosity can be thought of as either a von Neumann type viscous pressure [2], a source of entropy [9], or wave jumps across suitably defined eigenvectors as an analogy with Riemann solvers [20, 14]. A novel approach was taken by Inutsuka [21] who developed an artificial-viscosity-free SPH scheme, by a reformulation of the SPH convolution integrals. The scheme introduces the ∗

Corresponding author Email address: [email protected] (Kunal Puri)

Preprint submitted to Elsevier

April 3, 2014

minimum and sufficient amount of dissipation in the equations by using an iterative solution to a non-linear Riemann problem (Godunov’s method) at an imaginary interface between an interacting particle pair. The scheme was dubbed the Godunov SPH (GSPH) and is attractive because it eliminates parametrization and hence, user intervention, typically associated with the SPH artificial viscosity. Moreover, the dissipation introduced by GSPH can be thought of as a “generalized” SPH artificial viscosity. Indeed, the popular “signal” based artificial viscosity proposed by Monaghan [20] can be interpreted as a GSPH viscosity using the Roe’s approximate Riemann solver (see Sec. 4). This relation to artificial viscosity aids in understanding the nature of solutions generated with GSPH. Conversely, the investigation into approximate Riemann solvers can be used to study the artificial viscosity employed in SPH. The GSPH method has been successfully applied to problems in compressible hydrodynamics [21, 22, 23] and magneto-hydrodynamics with the Method of Characteristics (MOC) [13]. We have found GSPH to be robust and as accurate as the Monaghan, Price and Morris (MPM) scheme [20, 14] for the compressible Euler equations [24]. In GSPH, a one-dimensional Riemann problem is solved between interacting particles, by considering them to define equivalent left and right states abutting an imaginary interface normal to the line joining them. The Riemann solver is used to determine the intermediate pressure, p∗ and velocity, u∗ , which is used in the GSPH equations of motion. A time consuming operation for the GSPH solution procedure is the exact solution of the nonlinear Riemann problem for every interacting particle pair. The number of neighbours for a particle scale with the smoothing length, which is usually adaptively adjusted to the local   d1 particle volume as h ≈ η mρ , where d is the number of dimensions. For an SPH kernel, the number of neighbours N can be heuristically estimated as 2ζη, π(ζη)2 and 43 π(ζη)3 in one, two and three dimensions respectively, where ζ is the kernel scaling factor (ζ = 3 for the Gaussian) [14]. Comparing with finite volume computations, GSPH requires the solution of an order of magnitude more Riemann problems per computational element. Furthermore, the solutions are weighted with the kernel gradient which means that we must expend the same computational effort for neighbours sufficiently far away although they have a minimal contribution. Thus, significant gains in efficiency can be obtained if we forgo the exact solution of the Riemann problem in favour of an approximation, subject to the same stability and accuracy requirements. Approximate Riemann solvers, pioneered by Roe [25] have been widely used by the CFD community as an efficient alternative to exact solutions. Moreover, the exact solution to the non-linear equations is often unknown or intractable which necessitates approximate solutions. This is our primary motivation for an investigation into using approximate Riemann solvers for GSPH solutions to the compressible Euler equations. CFD is abound with approximate solvers for the Euler equations and these are historically well documented [26, 27, 28]. For our GSPH implementation, we are interested in the approximate solution to the Riemann problem in Lagrangian coordinates. To this end, we use the Lax-Friedrichs [29], HLLE [30], Roe [25], Ducowicz [31] and HLLC [32, 33] approximate Riemann solvers when applied to canonical test problems in 1 and 2 dimensions. While this paper was under review, Sirotkin and Yoh [34] have published results for SPH with approximate Riemann solvers for simulating explosions. Their derivation is based 2

on discretizing the divergence of the flux vector using SPH approximation techniques and subsequently using the Lax-Friedrichs and HLL approximate Riemann solvers to evaluate the fluxes between particle pairs. This differs from the approach adopted herein where we use the GSPH framework to introduce and study the dissipation introduced through approximate Riemann solvers. We would like to point out that Vila [35] had developed a particle weighting scheme (SPHALE) which is commonly used for SPH simulations for incompressible flows. In these schemes, Riemann problems are solved on moving control volumes similar to the Finite Volume Particle Method [36]. Typical implementations of SPH-ALE use the HLLC approximate Riemann solver [27] with second order, MUSCL-type reconstructions [37, 38, 39]. The principle advantage sought from the SPH-ALE scheme is a robust treatment of boundary conditions for incompressible flow calculations. In contrast, our test problems are idealized in the sense that while the boundary conditions are simpler, the compressible nature of the flow imposes a more stringent requirement on the discretization. We therefore do not consider the SPH-ALE scheme in our discussions. This paper is organized as follows. In Sec. 2, we describe the GSPH procedure, motivating the need for the solution to the Riemann problem. Sec. 3 briefly describes the different Riemann solvers we have used. We proceed to a discussion on the results using these schemes in Sec. 5. Finally, we conclude with a summary and suggestions for using approximate solvers for different problems in Sec. 6. In the interest of reproducibility, we make available our code (SPH2D)1 that was used in the preparation of this manuscript. 2. SPH and Godunov SPH (GSPH) We are interested in solving the Euler equations of gas dynamics, expressing the conservation of mass, momenta and energy of an ideal inviscid gas: dρ = −ρ (∇ · v) (1a) dt dv 1 = − (∇p) (1b) dt ρ de p = − (∇ · v) (1c) dt ρ The equations are closed with an equation of state (EOS). In this work, we use the ideal gas equation: p = (γ − 1)ρe (2) ∂ In the above equations, dtd = ∂t + v · ∇ is the material derivative, e is the specific thermal energy and γ is the adiabatic index. In SPH, points or “particles” are used as the computational elements for the numerical discretization of Eqs. 1. The discretization is written as a weighted summation of pairwise interactions between particles, which can be thought to represent the “flux” contributions as expressed by the governing equations. SPH as an interpolation technique is used to represent field functions through convolutions 1

SPH2D available at https://bitbucket.org/kunalp/sph2d

3




N(a)


h

a

Figure 1: The set of neighbors (Na ) for a given particle (a). The radius of interaction is determined by the kernel scale parameter ζ and smoothing length h.

with a smooth approximation to the Dirac delta distribution function, referred to as the kernel. We refer the interested reader to [40] and [14] for a detailed exposition on the SPH kernel interpolation theory. An important example of kernel interpolation is the estimation of the density field represented by a set of particles. This is expressed as:  ρ(x) = mb W (x − xb , h), (3) b∈N (x)

In this equation, W is a positive, symmetric weighting function (kernel) that has the dimensions of inverse volume. The equation is to be interpreted in the following sense. The density at a point x, is estimated by a weighted sum of individual particle masses in a neighborhood N (x), about the point. The kernel parameter ζ, and particle smoothing length or smoothing scale [41] h, determine the number of neighbors over which to perform the summation as shown in Fig.1. In this work, we use the truncated Gaussian kernel (ζ = 3):  −q2 1 0≤q≤3 e W (q, h) = d (4) d 0 otherwise π2h is the normalized interaction distance, h is the characteristic smoothing length where q = |x| h and d is the number of spatial dimensions. 2.1. Godunov SPH The discrete equations for GSPH can be derived by using either the kernel interpolation theory or via an action principle [21]. We follow the former approach to elucidate the steps 4

leading to the introduction of the Riemann solver. A fundamental construct in GSPH is to utilize Eq. 3, to obtain the following discrete consistency identities:  mb W (x − xb , h) (5a) 1= ρ(x) b∈N (x)    W (x − xb , h) (5b) mb ∇ 0= ρ(x) b∈N (x)

Using these identities, the SPH approximation of a function can be recast as  < f > (x) ≈ f (x )W (x − x, h)dx Ω  f (x ) W (x − x, h)W (x − xb , h)dx = mb  ρ(x ) Ω

(6)

b∈N (x)

where we have used the first consistency condition in Eqs. 5, to result in a modified SPH convolution integral. Inutsuka [21] has pointed out that this expression reduces to the standard SPH interpolation formula if we make the crude approximation W (x − xb , h) ≈ δ(x − xb , h). From this perspective, the GSPH interpolation is similar to restoring particle consistency in the discrete equations. In like manner, for the equations of motion, equations (1b) and (1c) are convolved and subsequently augmented with the consistency conditions to give   p(x) ¨a = − x (W (x − xb )∇W (x − xa , h) − W (x − xa )∇W (x − xb , h)) dx mb ρ2 (x) b∈N (a)

e˙ a = −

 b∈N (a)

 mb

(7a) p(x) (v − x˙ ∗a ) (W (x − xb )∇W (x − xa , h) − W (x − xa )∇W (x − xb , h)) dx ρ2 (x) (7b)

¨ is the time centred velocity for a particle. The density at x In these equations, x˙ ∗ = x˙ + Δt 2 the position of a particle is estimated as  ρa = mb Wab (ha ) (8) b∈N (a)

which is a direct application of Eq. 3 at the position of the particle. The introduction of the time centered velocity in the GSPH energy equation 7b, represents a “mixing” for the internal energy which is absent for standard SPH. This term is also needed for the discrete conservation of the energy. 5

To arrive at the discrete equations of motion, we need to evaluate the integrals in Eqs. 7. In the original GSPH method, Inutsuka [21] proposed to approximate the integrals in two steps. First, a local coordinate system (s, s⊥ ) is set-up, with the origin at the mid-point of the line joining two interacting particles and unit vector along the line of sight given by 1 b rˆ ab = |xxaa −x . An interpolation for the specific volume V(s) = ρ(s) in this local system is −xb | used obtain specific volume contributions Vab as  √ 1 2 W (x − x , h)W (x − x , h)dx ≈ V (h)W (x − x , 2h) (9) b a a b ab 2 Ω ρ(x) The algebraic manipulations leading to this equation are only valid for the Gaussian kernel. √ The scaling factor ( 2) for the smoothing length also means we must use a larger search radius for each particle. The specific volume contributions can then used to evaluate the integral for an arbitrary function f (x) as  √ f (x) ∗ W (x − xb , h)W (x − xa , h)dx ≈ fab Vab2 (h)W (xa − xb , 2h) (10) 2 Ω ρ(x) ∗ is to be interpreted as the interpolated value of the function at some position s∗ab , where, fab along the line of interaction. Using these expressions, the GSPH equations Eqs. 7 can be written as

 √ √ ¨a = − x mb p∗ab Vab2 (ha )∇Wab ( 2ha ) + Vab2 (hb )∇Wab ( 2hb ) (11a) b∈N (a)

e˙ a = −





√ √ mb p∗ab [v ∗ab − x˙ ∗a ] · Vab2 (ha )∇Wab ( 2ha ) + Vab2 (hb )∇Wab ( 2hb )

(11b)

b∈N (a)

In these equations, the starred quantities (p∗ab , v ∗ab ) are the intermediate state arising from ∗ the solution of a Riemann problem and like fab in Eq. 10, represent some interpolated value of the variables in the local coordinate system, at the position s∗ab . Note that in GSPH, we solve one-dimensional Riemann problems along the line joining two particles with the particle a (respectively b) defining the right (left) state in the local coordinate system. For input to the Riemann solver, the particle velocities are projected upon rˆ ab as ur = v a · rˆ ab ul = v b · rˆ ab

(12a) (12b)

The resulting scalar intermediate velocity u∗ab , is used to define the vector intermediate velocity used in Eq. 11b as   v a + v b u l + ur ∗ ∗ − rˆ ab , (13) v ab = uab rˆ ab + 2 2 where we have assumed s∗ab = 0 for simplicity [21]. Note that the bracketed term is perpendicular to the interaction vector rˆ ab , and is therefore identically zero in Eq. 11b. 6

The Vab terms depend on the kind of interpolation used for the specific volume in the local coordinate system, which can be either a linear or cubic spline interpolation as originally proposed by Inutsuka [21]. Cha and Whitworth [23] and recently, Iwasaki and Inutsuka [13] have suggested the relatively simple approximation ∇W (x−xa )W (x−xb )−∇W (x−xb )W (x−xa ) ≈ ∇(x−xa )δ(x−xb )−∇W (x−xb )δ(x−xa ) (14) which, when used in Eq. 7 readily gives the evolution equations as   1 1 ∗ ¨ xa = − mb pab ∇Wab (ha ) + 2 ∇Wab (hb ) (15a) ρ2a ρb b∈N (a)   1 1 ∗ ∗ ∗ mb pab [v ab − x˙ a ] · ∇Wab (ha ) + 2 ∇Wab (hb ) (15b) e˙ a = − ρ2a ρb b∈N (a)

Note that this formulation admits the use of different SPH kernels and avoids evaluating the specific volume integrals Eq. 9. Moreover, the resulting evolution equations are very similar to the popular Monaghan-Price-Morris (MPM) formulation [14]. The results using this simplified formulation were shown to be comparable with finite volume codes for one and two dimensional problems in hydrodynamics and magneto-hydrodynamics [23, 13]. This suggests a critical role played by the Riemann solver in the equations of motion. This is discussed next. 2.2. The role of the Riemann solver From the derivation of the GSPH equations of motion (Eqs. 11), it is apparent that the starred variables (p∗ , v ∗ ), represent the interpolated pressure and velocity, at some point along the line joining an interacting particle pair. In effect, there is no compulsion to consider them as the intermediate state generated as a solution to the Riemann problem between the two particles. The ansatz however, provides the necessary and sufficient dissipation needed to stabilize the scheme. Monaghan [20] had originally derived the “signal” based artificial viscosity by using a loose analogy with Riemann solvers. In fact, this viscosity can be interpreted as the diffusive contributions from an appropriate Riemann solver (see Sec. 3.2). Hitherto, GSPH implementations [21, 13, 22, 6] have relied on the solution to the exact Riemann problem to determine this intermediate state. The exact solvers are robust but require iterative root finding techniques that are computationally inefficient when compared with non-iterative approximate Riemann solvers. Moreover, exact solutions to the Riemann problem are not available for complex equations of state and/or multi-material problems. This further motivates the need for approximate Riemann solvers. We expect the approximate solution to work for most, if not all scenarios and for general equations of state. Furthermore, the approximate solver should introduce the necessary and sufficient dissipation, thereby keeping GSPH artificial-viscosity-free. We have several choices for the solver which have been used in the CFD community for hyperbolic systems of conservation laws [25, 42, 43, 44, 31, 29, 26, 33, 32]. The approximate Riemann solvers we consider for coupling with the GSPH scheme are described next. 7

3. The Riemann problem and approximate Riemann solvers For the Euler equations of gas dynamics, the one-dimensional Riemann problem is defined by two input states (left, right) abutting an interface, which results in a self-similar solution along lines with constant xt . The solution is composed of different states, delineated by waves which emanate from the interface. The waves are ordered by their speeds λi , so that λ1 < λ2 < λ3 as shown in Fig 2. Of the three waves, the slowest and the fastest




    



Figure 2: Wave structure for the Euler equations. Three waves emanate from an initial discontinuity with left (L) and right (R) states. The middle wave is always a contact discontinuity while the slowest and fastest waves can be either a shockwave or an expansion wave.

waves can be either a shockwave or an expansion wave, while the intermediate wave is called a contact discontinuity. All variables are smooth and discontinuous across an expansion wave and shockwave respectively. The contact discontinuity is special because the pressure and normal component of velocity are constant while the density and energy change discontinuously across it. The contact discontinuity is a linear wave which travels with the local fluid velocity. In the Lagrangian context, since the material is assumed to move with the fluid velocity, we have only two waves and the solution is uniquely determined by obtaining the pressure and velocity in the intermediate state between the fastest and slowest 8

waves. Non-linear equations, relating the states fore and aft of these waves can be obtained and subsequently solved by a root finding algorithm such as the Newton-Raphson method [45, 46, 27]. While highly accurate, these methods are inefficient when we consider the fact that we only seek the intermediate states. The resolution and effective decomposition of the non-linear wave structure that was used to find this state is never used in the solution. Consequently, approximate Riemann solvers were introduced [25, 31, 30, 42, 29, 26] to efficiently provide a reliable estimate for the intermediate states at a fraction of the cost. Approximate solvers are also a necessity when the exact solution to the Riemann problem is not known or impractical for computational purposes. While we do not expect approximate solvers to be robust for all conceivable scenarios, they are expected to provide an efficient approximation for most cases. Their widespread use in the CFD community is testimony to this fact. As was outlined in Sec. 2, the evolution equations for GSPH require, for each interacting particle pair, an interpolated pressure p∗ and projected velocity u∗ , at an imaginary interface along the line joining two particles. The original idea of Inutsuka [21] was to use the exact solution to the Riemann problem that results when the two particles determine the left and right states astride an interface. With the aim to reduce the computational time spent in solving the exact Riemann problem for every interacting particle pair, we explore the possibility of using approximate solvers to determine the intermediate states. We use the solvers reviewed by Rider [26], which are used in Lagrange plus Eulerian remap numerical codes. In addition to these, we also test the HLLC [32, 33] and Ducowicz [31] approximate solvers which were found to be suitable by Cheng and Shu [47] in their ENO-type, Lagrangian Finite Volume Method. The Euler equations in Lagrangian variables, as considered by Rider [26] are given as τt − uξ = 0 u t + pξ = 0 eˆt + (pu)ξ = 0

(16a) (16b) (16c)

where τ = 1/ρ, is the specific volume, eˆ = 12 u2 + e, is the total energy per unit mass and the time derivative is to be understood as a derivative moving with the fluid (material derivative) d ∂ (∗) = ∂t (∗) + u∇(∗). ξ is the mass coordinate, which is related to the spatial coordinate x, dt through the transformation dξ = ρdx. The system is hyperbolic with eigenvalues λ1 = −C, 1 λ2 = 0 and λ3 = C, where C = (γpρ) 2 is the Lagrangian sound speed, which is related to the Eulerian sound speed through C = ρc. We proceed to a brief description of the different approximate Riemann solvers. We would like to mention that among the available Riemann solvers, we have not considered the Engquist-Osher [42] and the flux splitting variants of approximate Riemann solvers [26]. This is because the absence of a sonic point in the Lagrangian equations negates the advantages of using the Engquist-Osher solver and it is not clear how to derive the intermediate pressure and velocity needed for GSPH, from the flux-split AUSM approximate Riemann solver.

9

3.1. Local Lax-Friedrichs approximate Riemann solver (LLXF) The local Lax-Friedrichs Riemann solver [29, 26, 47] is perhaps the most efficient of the approximate Riemann solvers. It can be intuitively understood as a flux contribution that results from an unstable centred approximation and a stabilizing diffusive contribution 1 F ∗ = [F l + F r − ηlr (U r − U l )], 2

(17)

where ηlr is the maximum (Lagrangian) wave speed at the interface. This speed can be taken as an interface constant equal to the maximum sound speed [26], or the maximum eigenvalue for each component in a characteristic decomposition [47]. Efficiency though, comes at a price of being highly diffusive. The intermediate pressure and the product with velocity, computed from Eqs. 17 and Eqs. 16 are 1 [pl + pr − ηlr (ur − ul )] 2 1 (pu)∗ = [pl ul + pr ur − ηlr (ˆ er − eˆl )] 2 p∗ =

(18a) (18b)

2

In these equations eˆ = e + u2 is the total energy per unit mass and we use Roe-averaging to define the averaged sound speed ηlr . The intermediate velocity is taken as u∗ = (pu)∗ /p∗ . 3.2. Roe’s approximate Riemann solver (ROE) Roe [25] proposed to linearise the Euler equations and thereby reduce it’s solution to that of a system of linear equations. The characteristic decomposition of the linearised Lagrangian flux and Jacobian matrices are used to write the flux as  1 F ∗ = [F l + F r − r klr |λklr |(αr − αl )] 2 k

(19)

where the summation is over k characteristic waves (eigenvectors) r, with speeds λ and strengths (wave jumps) αr − αl . An interesting point of this solver is that the form of the diffusive contributions were used by Monaghan [20] as an inspiration to derive the “signal” based artificial viscosity for SPH. The computation of the ROE flux requires a suitable linearisation. Rider [26] proposes an algebraic averaging of the Lagrangian variables (p, τ ) which we adopt in this work. The resulting intermediate states for the pressure and velocity bear a striking resemblance to those of the LLXF approximate Riemann solver.   1 1 ∗ u = u l + ur − (pr − pl ) (20a) 2 Clr 1 p∗ = [pl + pr − Clr (ur − ul )] (20b) 2 where Clr is the averaged Lagrangian velocity at the interface. We refer the reader to [26] for the right and left eigenvectors for the particular linearisation, although we do not need them. 10

3.3. Harten, Lax, van Leer and Einfeldt approximate Riemann solver (HLLE) A theoretical approximate Riemann solver which assumes that two waves emanate from the Riemann problem was proposed by Harten, Lax and van Leer [43] and further developed into a numerical scheme by Einfeldt [48, 30]. Just like the ROE and LLXF approximate Riemann solvers, this solver can also be viewed as a central plus a diffusive term. The intermediate pressure and velocity at the interface is given by F∗ =

S M F l − Sm F r SM Sm + (U r − U l ) SM − Sm SM − Sm

(21)

where, SM and Sm are estimates for the minimum and maximum signal speeds. The intermediate states are given as   1 SM p l u l − S m p r u r S M Sm ∗ u = ∗ + (ˆ er − eˆl ) (22a) p SM − Sm SM − Sm S M p l − Sm p r SM Sm p∗ = + (ur − ul ), (22b) SM − S m S M − Sm which are identical to the LLXF intermediate state for the choice of signal speeds SM = max(Cl , Cr ) and Sm = −SM . For our calculations, we use the following SM = max(ur + Cr , C˜lr ) Sm = min(ul − Cl , −C˜lr ),

(23a) (23b)

where C˜lr denotes a the Roe-averaged, Lagrangian sound speed at the interface. Remark 1. We note that LLXF, ROE and HLLE approximate Riemann solvers are similar considering their general form as a “centred”, unstable term plus a stabilizing “diffusive” contribution. The differences are in the kind of averaging used to define the intermediate states and the diffusion contributions that arise from the differences in the conserved variable. 3.4. Ducowicz approximate Riemann solver (DUCO) Ducowicz [31] introduced a general non-iterative Riemann solver which is characterized by two material dependent parameters, which are the local sound speed and the limit of the shock density ratio for strong shocks. Ducowicz highlights the similarities between artificial viscosity and Godunov schemes as a means to introduce entropy into the isentropic equations. Specifically, the shock Hugoniot curve and it’s relation to the form of the artificial viscosity is exploited to construct the approximate Riemann solver. The solver is expected to be as general as artificial viscosity methods for arbitrary equations of state. The Ducowicz approximate Riemann solver is formally a two shock approximation where transitions through rarefaction waves are treated as transitions through rarefaction shocks. The interface velocity u∗ is obtained by solving the following equation ρ+ A+ |u∗ − u∗min |(u∗ − u∗min ) + ρ− A− |u∗ − u∗max | (u∗ − u∗max ) + p+ − p− = 0 11

(24)

+

where A− are the parameters related to the shock density ratio, which, for ideal gasses is γ+1 2 [47]. u∗min and u∗max are estimates for the minimum and maximum interface speeds that can be obtained from the local Lagrangian sound speed at the interface. The interface pressure p∗ is computed from the velocity as 1 1 1 p∗ = (p+ + p− ) + p+ A+ |u∗ − u∗min |(u∗ − u∗min ) − p− A− |u∗ − u∗max |(u∗ − u∗max ) 2 2 2

(25)

Ducowicz [31] provides a FORTRAN subroutine that can be adapted to the Lagrangian formalism for GSPH. 3.5. Harten, Lax and van Leer scheme with contact (HLLC) The HLLC solver was introduced by Toro et al. [49] as a 3-wave approximate Riemann solver that had better resolution of the contact discontinuity. This is a popular Riemann solver used for example in SPH-ALE schemes ([37, 38, 39]) and recent multi-dimensional extensions for Eulerian FVM schemes developed by Balsara [50] and Balsara et al. [51]. For our purposes, we consider the formulation of the HLLC solver in Lagrangian coordinates, proposed by Lou et al. [32] in their ALE scheme for multi-material, underwater explosions. The solver was also used by Cheng and Shu [47] in their ENO-type, Lagrangian Finite Volume Method. This approximate Riemann solver has the desirable properties of being positivity preserving for scalar quantities, entropy satisfying and to exactly preserve isolated contacts. The expressions for the intermediate pressure and velocity can be written as ⎧ ⎪ pl , if Sl > 0 ⎪ ⎪ ⎪ ⎨ SM [(S − v )M + (ˆ p − pl )] + pˆ, if Sl ≤ 0 < SM l l l M p∗ = Sl S−S (26a) M ⎪ [(S − v )M + (ˆ p − p )] + p ˆ , if S ≤ 0 < S r r r r M r ⎪ Sr −SM ⎪ ⎪ ⎩p , if Sr < 0 r ⎧ ⎪ p l ul , if Sl > 0 ⎪ ⎪ ⎪ ⎨ SM [(S − v )E − p v + pˆS ] + (S − v )ˆ if Sl ≤ 0 < SM l l l l l M M lr p, M (pu)∗ = Sl S−S (26b) M ⎪ [(S − v )E − p v + p ˆ S ] + (S + v )ˆ p , if S ≤ 0 < S r r r r r M M lr M r ⎪ S −S r M ⎪ ⎪ ⎩p u , if Sr < 0 r r These equations require an estimate for the interface velocity vlr , which can be taken as the Roe-averaged velocity on either side of the interface [47]. The quantities vl and vr are velocities relative to the interface velocity, pˆ = ρl (vl − Sl )(vl − SM ) + pl and the signal speeds SM , Sl and Sr are given by ρr vr (Sr − vr ) − ρl vl (Sl − vl ) + pl − pr ρr (Sr − vr ) − ρl (Sl − vl ) Sl = min(vl − cl , −clr ) Sr = max(vr + cr , clr )

SM =

12

(27a) (27b) (27c)

where clr is the Roe-averaged sound speed at the interface. The signal speeds are the wave speed estimates for the three waves considered in the solution to the initial Riemann problem.

Remark 2. An HLLC approximate Riemann solver in Lagrangian coordinates which has zero dissipation for pure convection was proposed by G. Ball [33] for his Free-Lagrange method (see also [52]). We find the results for this version of HLLC (referred herein as HLLC-B) to perform especially well for smooth solutions like the diffusion of a sound-wave (Sec. 5.1) and the two-dimensional accuracy test for convection (Sec. 5.3). 3.6. Numerical efficiency The efficiency of the approximate Riemann solvers is gauged by the time taken to return the intermediate states (p∗ , u∗ ) for representative shock-tube initial data in Table 1. The Test Shocktube Blastwave Sj¨ ogreen’s test

ρl 1 1 1

pl 0.125 1 1

ul 0 0 -2

ur 0 0 2

pl 1.0 1000.0 0.4

pr 0.1 0.01 0.4

Table 1: Sample input to test the efficiency of the approximate Riemann solvers.

three cases correspond to the initial data for the Sod’s shock-tube problem, the blast-wave problem and the Sj¨ogreen’s strong rarefaction test (1-2-3 problem) [27]. The reason for choosing three cases is because for the Newton-Raphson iteration in the exact Riemann solver, the time depends on the convergence tolerance ( = 10−6 ) and maximum number of iterations (niter = 10) which is determined by the input data. For each case, we call the function to return the intermediate states (p∗ , u∗ ) one million times. Results from Ducowicz Roe HLLC HLLE LLXF

Shocktube

Blastwave

Sj¨ ogreen’s test

0.78 0.49 0.57 0.57 0.49

0.63 0.37 0.45 0.45 0.40

1.17 0.72 0.85 0.86 0.77

Table 2: Relative times taken for different Riemann solver functions when returning the intermediate states (p∗ , u∗ ) for three representative inputs. We see that the approximate Riemann solvers are generally more efficient than the iterative exact solver. The Roe’s approximate Riemann solver is the most efficient, being approximately thrice as fast for the blast-wave problem.

this comparison are shown in Table 2, where the average time is displayed as a fraction of the time taken for the reference van Leer exact Riemann solver. We see that the Roe’s approximate Riemann solver is the most efficient for all three representative inputs. The HLLE and LLXF solvers, which differ in our implementation on the choice of signal speeds are equally matched. The Ducowicz and the HLLC Riemann solvers, although non-iterative, 13

require multiple if-else conditions which reduces their efficiency. As a disclaimer, we want to stress that the numbers presented herein are only representative and specific to our implementation. 4. Relation to the SPH artificial viscosity In this section, we establish an equivalence between the dissipative terms in GSPH and the “signal” based artificial viscosity [20, 53, 14] (referred to herein as the Monaghan-PriceMorris (MPM) scheme), under the assumption of a general “centred” plus “diffusive” approximate Riemann solver. To establish the notation used, for a pair of interacting particles (a, b), positioned at xa and xb respectively, we define a local coordinate system (s, s⊥ ) b) b with the origin at (xa +x and with unit vector rˆ ab = |xxaa −x . The GSPH equations solve 2 −xb | a one-dimensional Riemann problem in this local coordinate system with particle velocities projected upon rˆ ab , with particle a (respectively b), defining the right (left) state: ur = v a · rˆ ab ul = v b · rˆ ab The intermediate velocity used in the GSPH energy equation is given by Eq. 13:   v a + v b u l + ur ∗ ∗ v ab = uab rˆ ab + − rˆ ab , 2 2

(28a) (28b)

(29)

where, u∗ab is the intermediate projected velocity returned by the Riemann solver. Note that the bracketed term is orthogonal to the line joining the two particles. In particular, for the kernel gradient expressed as ∇Wab = rˆ ab Fab , we have v ∗ab · ∇Wab = u∗ab Fab . For simplicity, we consider the GSPH equations with constant smoothing lengths although the discussion applies to the case as well. In addition, we make a slight abuse of  variable-h  1 1 2 notation and define Vab = ρ2 + ρ2 for the remainder of this section. With these conditions a b and established notation, the GSPH Eqs. 15 read  ¨ a = −2 mb p∗ab Vab2 ∇Wab (h) (30a) x b∈N (a)

e˙ a = −2



mb p∗ab [u∗ab rˆ ab − x˙ ∗a ] · Vab2 ∇Wab (h),

(30b)

b∈N (a)

We now make the assumption that the general form of the intermediate states p∗ab and u∗ab can be given as a centred plus diffusive contribution: 1 1 u∗ab = u¯ab − (pa − pb ) (31a) 2 ρ¯ab c¯ab 1 (31b) p∗ab = p¯ab − ρ¯ab c¯ab (ua − ub ), 2 where the bar above the variables denotes suitably averaged quantities at the interface. As an example for Roe’s (Sec. 3.2)approximate Riemann solver, we use a simple arithmetic averaging a ¯ab = 12 (aa + ab ). We now proceed to an analysis of the dissipation introduced by these estimates in the GSPH equations of motion. 14

4.1. Dissipation in the Momentum equation Substituting in the momentum Eq. 30a, the second term (the diffusive contribution) in Eq. 31b, we get  GSPH  1 dv a =2 mb ρ¯ab c¯ab (ua − ub )Vab2 ∇Wab dt diss 2 b∈N (a)  = mb ρ¯ab Vab2 c¯ab v ab · rˆ ab ∇Wab b∈N (a)

=

 

ρ¯2ab Vab2

b∈N (a)

 mb c¯ab v ab · rˆ ab ∇Wab ρ¯ab

(32)

Comparing this expression with the signal based based artificial viscosity introduced by Monaghan [20, 14] MPM   mb dv a = α vsig v ab · rˆ ab ∇Wab , (33) dt diss ρab b∈N (a)

we can readily identify the GSPH dissipation as an artificial viscosity with αGSPH = ρ¯2ab Vab2 GSPH and signal velocity vsig = c¯ab . A difference is that while the traditional SPH artificial viscosity is activated only when v ab · rˆ ab ≤ 0, the GSPH viscosity (and indeed, Godunov methods in general) makes no such distinction. Another difference is in the definition of the signal velocity, which, in the MPM scheme is defined as ([14]) MPM vsig =

1 [ca + cb − βv ab · rˆ ab ] , 2

(34)

with the β parameter introduced to prevent particle penetration in high Mach number collisions. This term effectively increases the magnitude of the viscosity proportional to |v ab · rˆ ab |. GSPH does not introduce this term, however, it is compensated by the relatively larger value for αGSPH . Consider the situation for two particles with ρb = f ρa , where f > 0, is a scalar value. The effective value of the GSPH viscosity parameter is then  1 1 1 GSPH 2 α = (ρa + ρb ) + 4 ρ2a ρ2b  1 1 1 2 + = (ρa + f ρa ) 4 ρ2a f 2 ρ2a 1 = 2 (1 + f )2 (1 + f 2 ) (35) 4f For the case of a head-on collision between two equal density particles (f = 1), the value of the parameter is 2. In general the value adapts to the jump in density according to Eq. 35. Thus, the relatively higher values for the effective GSPH viscosity parameter helps to prevent particle penetration. Indeed, Cha and Whitworth [23] have performed subsonic and supersonic flow collision tests with GSPH and observed no penetration. In some extreme 15

cases like the Woodward and Colella blast-wave problem (Sec 5.5), involving the collision and reflection of strong shocks, it might be necessary (depending on the choice of the Riemann solver) to add some additional dissipation to prevent particle penetration. In such cases (which are to be determined by numerical experiment), we can modify the definition of c¯ab in Eq. 31, to include an MPM-like β term. 4.2. Dissipation in the energy equation To analyse the dissipation in the GSPH energy equation 30b, we write it as  e˙ a = −2 mb p∗ab u∗ab Vab2 Fab b∈N (a)

+ va · 2



mb p∗ab Vab2 ∇Wab

b∈N (a)

 Δt ¨a · 2 x + mb p∗ab Vab2 ∇Wab , 2

(36)

b∈N (a)

¨ a . The where we have expanded the definition of the time-centred velocity x˙ ∗a = v a + Δt x 2 energy equation is split into three terms. In general, the dissipation enters through all the terms because of the variables p∗ab and u∗ab . We recognize the third term as a higher-order (because of the presence of the Δt factor) dissipation, proportional to the square of the ¨ 2a . magnitude of particle acceleration and is given by − Δt x 2 The first term in Eq. 36 introduces the product p∗ab u∗ab , which is written as   1 ∗ ∗ pab uab = (pa + pb )(ua + ub ) 4   1 pa + pb 1 − (pa − pb ) − (ua − ub )(pa − pb ) 4 ρ¯ab c¯ab 4   1 − ρ¯ab c¯ab (ua + ub )(ua − ub ) 4

(37)

Dissipation enters through the second and third rows of this equation (differences in the variables). The resulting dissipation for the GSPH energy equation (ignoring the high-order term) is  GSPH   dea 1  2 m b p a + pb =+ ρ¯ab Vab − (ua − ub ) (pa − pb )Fab dt diss 2 ρ¯ab ρ¯ab c¯ab b∈N (a)

1  2 2 mb + ρ¯ab Vab c¯ab (u2a − u2b )Fab 2 ρ¯ab b∈N (a)  mb − ρ¯2ab Vab2 c¯ab ua (ua − ub )Fab , ρ¯ab b∈N (a)

16

(38)

where we have used v a · ∇Wab = ua Fab and the diffusive contribution of p∗ab , to arrive at the last line of Eq. 38. We can group the last two lines of this equation (the viscous contribution to the thermal energy) to get 

dea dt

GSPH diss

  1  2 m b p a + pb =+ ρ¯ab Vab − (ua − ub ) (pa − pb )Fab 2 ρ¯ab ρ¯ab c¯ab b∈N (a)

1  2 2 mb ρ¯ab Vab c¯ab ((v a · rˆ ab ) − (v b · rˆ ab ))2 Fab − 2 ρ¯ab

(39)

b∈N (a)

The first term represents the thermal conduction which acts on pressure differences (pa −pb ), which is in turn proportional to the energy difference through the equation of state. The second term represents the viscous contribution to the thermal energy which is positive definite, since Fab < 0 for all kernels considered in this work. The dependence on energy differences can be made more explicit using the equation of state. For example, using the ideal gas equation, we have pa − pb = (γ − 1)(ρa ea − ρb eb ) = (γ − 1)ρab (fa ea − fb eb ),

(40)

where we have introduced scalar parameters ρa,b = fa,b ρab , to enable factoring out the averaged density. Interestingly, Monaghan [20] had also suggested these scaling factors for the thermal energy. Substituting for the pressure difference in Eq. 39 gives us   GSPH  dea 1  2 2 m b pa + pb =+ (γ − 1)¯ ρab Vab − (ua − ub ) (fa ea − fb eb )Fab dt diss 2 ρ¯ab ρ¯ab c¯ab b∈N (a)

1  2 2 mb − ρ¯ab Vab c¯ab ((v a · rˆ ab ) − (v b · rˆ ab ))2 Fab 2 ρ¯ab b∈N (a)

1  mb GSPH GSPH GSPH vsig (v ab · rˆ ab )2 − (γ − 1)αGSPH vsig,e (fa ea − fb eb ) Fab , α =− 2 ρ¯ab b∈N (a)

(41) where we have used the the following to define the GSPH viscosity parameters αGSPH = ρ2ab Vab2 GSPH vsig = c¯ab

(42a) (42b)

GSPH and introduced the signal velocity for thermal conductivity vsig,e (see Eq. 44 below). On comparison with the dissipation for the thermal energy in the MPM scheme: MPM   dea 1  mb  e αvsig (v ab · rˆ ab )2 − 2αe vsig =− (ea − eb ) Fab , (43) dt diss 2 ρ¯ab

b∈N (a)

17

we find that the GSPH thermal conduction coefficient is given by αeGSPH = γ−1 αGSPH and 2 acts on the scaled differences in thermal energy. Note that the equation implies a nonzero thermal conductivity for GSPH due to the dependence on ρ2ab Vab2 . It is worth noting that thermal conduction is essential to mitigate the spurious pressure “wiggle” at a contact discontinuity in SPH [20, 4, 54, 14]. It is this implicit thermal conduction in Godunov schemes that help suppressing the anomaly when it occurs. We define the equivalent signal velocity for thermal conduction as   p a + pb GSPH vsig,e = − (ua − ub ) ρ¯ab c¯ab 2 = c¯ab − v ab · rˆ ab (44) γ which remains positive in the vicinity of  shocks (v ab · rˆ ab ≤ 0). In the MPM scheme, Price |pa −pb | [54, 14] advocates the use of vsig,e = to smoothen the pressure jump across the ρ¯ab MPM is adaptively adjusted using the Laplacian of contact discontinuity and the value of αe the thermal energy. In GSPH, the conductivity adapts to the viscosity parameter αGSPH . The results in this section establish the equivalence of the GSPH dissipation terms to those of the popular signal-based artificial viscosity used in SPH. Indeed, Monaghan [20] derived the expressions for the signal-based viscosity by considering a loose analogy with Riemann solvers. The arguments however leading to their development were heuristic as pointed out by Monaghan himself. The equivalence is important as it places the hitherto heuristic derivation on a more formal footing and gives us a rational justification for the use of artificial viscosity in SPH. Additionally, the relation to artificial viscosity can be used to interpret the numerical solutions generated by GSPH from within the artificial viscosity framework. A limitation to the derivation carried out in this section is the restriction to the class of approximate Riemann solvers for which the intermediate states can be expressed as a general “centred” plus “diffusive” term (Eq. 31). From the solvers considered in this work, the Ducowicz (Sec. 3.4) and HLLC (Sec. 3.5) Riemann solvers do not fall into this category. 5. Numerical results In this section, we evaluate the performance of GSPH coupled with different approximate Riemann solvers which were introduced in Sec. 3. We are interested in solving the Euler equations of compressible gas-dynamics with an ideal gas equation of state. We can therefore compare the results using the approximate Riemann solvers with the van Leer [45] exact Riemann solver. The problems comprise a test suite that has been used previously [24, 4, 23, 55, 20] to assess the robustness and accuracy of numerical schemes for compressible hydrodynamics. For the GSPH results, we use Eqs. 15, with a second order reconstruction of the primitive variables (ρ, v, p) and the original monotonicity algorithm proposed by Inutsuka [21]. The 18

gradients of the primitive variables are estimated following Iwasaki and Inutsuka [13]  (∇ρ)a = mb ∇a Wab (ha ) (45a) b∈N (xa )

(∇Q)a =

 mb (Qb − Qa )∇a Wab (ha ), ρb

(45b)

b∈N (xa )

for Q = (u, v, p). Although equations (15) permit the use of any arbitrary kernel function, we find that the Gaussian kernel is the most accurate. Unless otherwise mentioned, we use the Gaussian kernel (Eq. 4) in all our calculations. The equations of motion are integrated using xn+1 = xna + Δtx˙ ∗a a v n+1 = v na + Δt¨ xa a n+1 n ea = ea + Δte˙

(46a) (46b) (46c)

with a constant time step. The use of the time-centred velocity x˙ ∗a , ensures discrete energy conservation [21]. The smoothing lengths are determined from the particle positions using two parameters (k ≈ 1, Csmooth = 2), in the following three step process  ρ∗a = mb Wab (xa − xb , Csmooth hna ) (47a) b∈N (a)

1 ma d =k ρ∗ a = mb Wab (xa − xb , hn+1 ) a 

hn+1 a ρn+1 a

(47b) (47c)

b∈N (a)

While displaying the results, we plot the particle weights to avoid any additional smoothing. For the density map for the 2D Riemann problem in Sec. 5.7, we interpolate the density of the particles onto a fixed grid using the Gaussian kernel. 5.1. Diffusion of an acoustic wave A disadvantage of the implicit dissipation in GSPH (Godunov schemes in general), is the inability to artificially turn it off when not required. An example is the propagation of an acoustic wave ([56, 57, 9, 34]). Cullen and Dehnen [57] have shown how SPH without artificial viscosity can propagate this wave in an ideal, dissipation-less manner, preserving the amplitude and phase. As a consequence of the implicit dissipation, we expect a progressive smearing of this wave with GSPH. For the numerical set-up, we distribute Np uniformly spaced particles in the periodic domain x ∈ [0, 1]. The wave is described by a small amplitude (Δρ = 10−6 ) sinusoidal perturbation, with particle properties satisfying the wave equation [34]: ρ = ρ0 + Δρ sin(κx) p = p0 + c20 Δρ sin(κx) u = c0 ρ−1 0 Δρ sin(κx) 19

(48a) (48b) (48c)

In these equations, κ = 2π/λ is the wave number and the wavelength λ is taken equal to the domain length. We take γ = 1.4 and consider ρ0 = γ and p0 = 1.0 as the background state. The sound speed in the undisturbed fluid is c0 = 1. We let the wave propagate for 10 periods and compute the L∞ norm of the error for the velocity with respect to the initial conditions. Table 3 displays these errors when using 400 particles, the Gaussian kernel L∞ × 10

−9

GSPH1 226.7

GSPH2 4.62

LLXF 4.60

ROE HLLE DUCO 4.62 4.60 22.5

HLLC HLLC-B 48.79 4.62

Table 3: Velocity L∞ norm errors (scaled by 10−9 ) for the propagation of a one-dimensional sound-wave. 400 particles have been used in the domain x ∈ [0, 1] and the wave is convected for 10 sound crossing times. The first order scheme (GSPH1 ) is over-diffusive. Among the approximate Riemann solvers, the Ducowicz and HLLC solvers are more diffusive for this test.

and a constant time step of Δt = 10−3 s. GSPH1 and GSPH2 denote the first and second order variants of the GSPH scheme respectively. We find that the HLLC and Ducowicz approximate Riemann solvers are fairly diffusive for this test while the LLXF, HLLE, Roe and HLLC-B solvers result in a decrease in amplitude similar to the reference second-order GSPH scheme. We believe the relatively low dissipation for Roe, LLXF, HLLE and HLLC-B is because the diffusive contributions in these solvers arise from differences in the conserved variables, which is of the order of the perturbation amplitude. 5.2. Accuracy test in 1D We use another one-dimensional problem that was used by Cheng and Shu [47], to test the accuracy of their Lagrangian Finite Volume Method. In an earlier comparison [24], it was found that the GSPH with the exact Riemann solver was approximately second order accurate and comparable in accuracy to the Monaghan, Price and Morris (MPM) scheme [14]. We repeat this test here to gauge the accuracy of the different approximate Riemann solvers when the solution remains smooth. The initial conditions for this test is a fluid in a periodic domain with uniform pressure and sinusoidal density and velocity profiles given by ρ(x) = 2 + sin(2πx),

u(x) = 1 + 0.1sin(2πx),

p(x) = 1

x ∈ [0, 1]

Since there is no analytical solution to this problem, we compare the results of the SPH schemes with a reference solution generated using a TVD MUSCL scheme ([27]), with 10000 cells at the time t = 1s. For the numerical set-up, the particles are distributed uniformly m in the periodic domain x ∈ [0, 1] and their masses set so that ρ = Δx reproduces the initial density profile. We examine the errors between SPH and the reference solution as the number of particles is increased. Since the solution remains smooth, this test can be used to gauge the relative magnitude of dissipation introduced by the different approximate Riemann solvers. Figure 3 displays the L1 norm errors (scaled by 10−3 ) for the density, velocity and pressure profiles for six different Riemann solvers when using the Gaussian kernel. We denote by GSPH1 and GSPH2 , the first and second order variants of the original 20

GSPH scheme using the van Leer exact Riemann solver. The L1 norm, defined for a property f , is normalized by the number of particles as Np 1  SPH L1 = |f − f FVM (xa )| Np a a

L1 × 10−3

25

ρ

25

u

25

24

24

24

23

23

23

22

22

22

21

21

21

20

20

20

2−1

2−1

2−1

2−2 6 2

27 Np

28

2−2 6 2

27 Np

28

2−2 6 2

(49)

GSPH1

p

27 Np

GSPH2 DUCO ROE HLLC HLLE LLXF

28

Figure 3: L1 norm errors (scaled by 10−3 ) for the 1D accuracy test using six different Riemann solvers and the Gaussian kernel. GSPH1 and GSPH2 denote the original GSPH scheme with the exact Riemann solver. All schemes exhibit the expected order of convergence. The HLLE and LLXF solvers are seen to be more diffusive.

We observe the expected rate of convergence of approximately 1 and 2 for the first and second order variants of the GSPH scheme respectively. Linear reconstruction improves the accuracy as expected. The solutions using the approximate Riemann solvers are generally comparable, with the HLLE and LLXF Riemann solvers being slightly more diffusive than the Ducowicz and Roe solvers. 21

Remark 3. Results for the accuracy test at first do not appear to be in consonance with those of the sound-wave propagation (Sec. 5.1), where the Ducowicz and HLLC solvers were found to be very diffusive. We believe this is because of the absence of a wave steepening behaviour due to the low amplitude of perturbation (Δρ = 10−6 ) for the sound-wave problem, which results in a negligible diffusive contribution for the LLXF and HLLE solvers. 5.3. Accuracy test in 2D In addition to the accuracy of the discretizations, particles have a settling tendency in higher dimensions which can break symmetry and corrupt results. It is well known [14] that the SPH equations of motion, derived from a Lagrangian, will attempt to drive particles in pressure equilibrium to a locally regular distribution when the density is not uniform. To check the accuracy of the proposed approximate Riemann solvers in such a scenario, we consider a known periodic solution to the Euler equations for an ideal gas in pressure equilibrium. The problem specifies a gas with a constant pressure and velocity in a doubly periodic domain ([0, 1] × [0, 1]), initialized with a density profile ([55, 58]) ρ(x, y, t) = 1 + 0.2sin(π(x + y))

(50)

We have used the values u = 1 and v = −1 and p = 1. With these values and pressure equilibrium, the particles should simply advect with their initial velocities and return to their original position at t = 1s. The solution remains smooth and as such it is a good test for the spurious motion generated by SPH in the absence of artificial viscosity. This spurious motion can be related to the discrete consistency requirements which GSPH attempts to restore through the volume integrals in Eq. 9. The erroneous motion is of the order accuracy with which these integrals are computed [13]. For the numerical set-up, we distribute particles on a uniform Cartesian grid with different resolutions (40 × 40, 80 × 80, 160 × 160). At the time t = 1, we check the accuracy of the different approximate solvers using the relative L1 norm for the density, calculated as NP 1  |ρ0,a − ρa )|, L = NP a=1 1

(51)

where ρ0,a is the initial density of the particle a. Figure 4 compares the errors for the Ducowicz, Roe, HLLC and HLLC-B approximate Riemann solvers, with the reference second-order GSPH scheme (GSPH2 ). In addition, we display the errors for the MPM scheme (generated using the NDSPMHD code [14]), with standard artificial viscosity parameters. We see that the results for the Ducowicz, Roe and HLLC-B approximate Riemann solvers are virtually identical to the reference GSPH scheme, which is in turn similar to the MPM scheme. The HLLC solver results in the same rate of convergence, albeit with quantitatively larger errors. The LLXF and HLLE solvers perform poorly for this test. Recall that the diffusive terms for these solvers are proportional to the difference in the conserved variables. Given the non-uniform density, this results in a large finite viscosity which eventually corrupts the solution. 22

ρ

24

L1 × 10−3

23

22

21 10 2

GSPH2 MPM DUCO ROE HLLC-Ball HLLC 211

212 Np

213

214

Figure 4: L1 norm errors (scaled by 10−3 ) for the density for the 2D accuracy test for five different Riemann solvers. GSPH2 denotes the original GSPH scheme with the exact Riemann solver. The Ducowicz, Roe and HLLC-B solvers are comparable to the reference GSPH solution. The HLLC approximate Riemann solver exhibits the same convergence rate but quantitatively larger errors.

Remark 4. From the the errors in Fig. 4, we see that SPH does not converge at the expected rate for this problem. This is a result of the spurious motion generated by SPH despite a pressure equilibrium. We expect that a non-conservative formulation like the relativepressure SPH (rpSPH [59]), that relies on pressure differences in the momentum equation will avoid this spurious motion. 5.4. Sj¨ ogreen’s test case Einfeldt et al. [30] found that non-physical solutions were produced by conservative difference schemes when the flow is predominantly kinetic. In their analysis of expansion flows in the presence of solid walls, they showed that no Godunov-type scheme based on a linearised Riemann solver (Roe’s scheme for example) is positively conservative. The problem is characterized by an inequality relating the magnitude of the Mach number and the adi23

abatic index γ, for which the physical solution is positive but linearised Riemann solvers wrongly predict cavitation. An idealized version of this test is called the Sj¨ogreen’s test case. The shock-tube like set-up consists of two regions (left, right), initially separated by a diaphragm, with states (ρl , pl , ul ) = (1, 0.4, −2) and (ρr , pr , ur ) = (1, 0.4, 2). These conditions result in two strong rarefactions that emanate from the initial discontinuity (without creating a vacuum) and a trivial contact discontinuity at the initial jump. In the CFD literature, the test goes by the alias of the Einfeldt rarefaction test or the 1-2-3 problem [27].

ρ

1.2

Exact

1.2

DUCO

1.2

1.0

1.0

1.0

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0.0 −0.4 −0.2

0.0

0.2

0.0 0.4 −0.4 −0.2

0.0

0.2

0.0 0.4 −0.4 −0.2

1.0

1.0

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.0

0.2

0.4

0.0

0.2

0.4

e

1.0

ROE

0.2 −0.4 −0.2

0.0

0.2

0.2 0.4 −0.4 −0.2

0.0

0.2

0.2 0.4 −0.4 −0.2

Figure 5: Numerical results (dots) for the Sj¨ogreen’s test at t = 0.1s using GSPH with the exact (left), Ducowicz (centre) and Roe’s (right) approximate Riemann solver with 200 particles. The density (upper panel) is well re-produced for this test for each of the schemes. All schemes poorly predict the thermal energy at the origin.

For the numerical set-up, we use 200 equal mass particles uniformly distributed in the domain x ∈ [−0.5, 0.5], to give the initial inter-particle spacing Δx = 0.005. We use a constant time step of Δt = 5 × 10−4 s and examine the solution at the time t = 0.1s. We note that it is particularly difficult for a numerical scheme to accurately predict the near-vacuum state at the origin for this test. In particular, the thermal energy at the origin is grossly 24

over-predicted by numerical schemes even if they do not fail by wrongly predicting cavitation. The numerical results (dots) for this test, are shown in Fig. 5 for the van Leer exact, Ducowicz and Roe’s approximate Riemann solvers. The HLLC Riemann solver with second order reconstruction failed for this test with the given resolution. From the figure, we see that the density profile (top panel) is well captured by each of the schemes. Notice the large errors in the thermal energy (bottom panel) at the origin which resembles a non-physical “heating” of the gas. Roe’s approximate Riemann solver produces the most heating for this case. Shen et al. [60] argue that the large thermal energy at the origin for this test is the result of viscosity induced entropy errors, proportional to the numerical dissipation of a scheme, thereby implying that more diffusive schemes would produce larger errors. This can be verified by considering the results for the MPM scheme [14], as we artificially increase the thermal conduction parameter αu in Eq. 43. This is shown in Fig. 6, where we

1.0

e

0.8

0.6

Exact MPM αu = 0

0.4

MPM αu = 0.5 MPM αu = 1.0 0.2 −0.4

−0.3

−0.2

−0.1

0.0 x

0.1

0.2

0.3

0.4

Figure 6: Numerical results (dots) for the Sj¨ogreen’s test at t = 0.1s using the MPM scheme with different values of the thermal conduction parameter αu . The additional dissipation results in larger heating errors at the origin.

25

have used fixed values of αu equal to 0 (blue), 0.5 (green) and 1.0 (red) respectively. Note that artificial viscosity will never be activated for the MPM scheme since v ab · rˆ ab > 0 for this problem. From the discussion on the GSPH dissipation and it’s relation to the SPH artificial viscosity in Sec. 4, we know that GSPH will always include an implicit dissipation through the Riemann solver. This explains the larger heating errors for the GSPH schemes in Fig. 5. Moreover, the relative dissipation introduced by each of the approximate Riemann solvers can be gauged from the results of the 1D accuracy test in Sec. 5.2. Consequently, we expect the slightly more diffusive HLLE and LLXF schemes to exhibit larger errors. This is indeed the case as shown in Fig. 7, which shows the density (top panel) and thermal energy (bottom panel) profiles for the exact, ROE and LLXF Riemann solvers. The LLXF

ρ

1.2

Exact

1.2

ROE

1.2

1.0

1.0

1.0

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0.0 −0.4 −0.2

0.0

0.2

0.0 0.4 −0.4 −0.2

0.0

0.2

0.0 0.4 −0.4 −0.2

1.0

1.0

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.0

0.2

0.4

0.0

0.2

0.4

e

1.0

LLXF

0.2 −0.4 −0.2

0.0

0.2

0.2 0.4 −0.4 −0.2

0.0

0.2

0.2 0.4 −0.4 −0.2

Figure 7: Numerical results (dots) for the Sj¨ogreen’s test at t = 0.1s using GSPH with the exact (left), Roe (centre) and LLXF (right) approximate Riemann solver with 200 particles at t = 0.1s. Note that the more diffusive LLXF scheme results in larger spurious heating at the origin.

solver results in larger heating errors for the thermal energy. While similar results with large errors in the thermal energy were also observed by Cha and Whitworth [23] in the context of GSPH, we would like to point out that the heating is a result of dissipation and is 26

also produced by Eulerian finite volume schemes based on exact and approximate Riemann solvers [55, 27, 28]. A drawback of the GSPH scheme for such problems is the inability to switch off the dissipation as can be done by setting the viscosity parameters to zero in the MPM scheme. 5.5. Woodward and Colella blast-wave problem This is a stringent test introduced by Woodward and Colella [61] to test the accuracy of their Piecewise Parabolic finite volume method (PPM [61]). This test involves multiple shock waves and contact discontinuities and their reflections and interactions. The test was designed to test the accuracy of the overall flow with respect to the resolution with which the discontinuities can be captured. The initial state is a quiescent gas (γ = 1.4) with uniform density in the domain x ∈ [−0.5, 0.5]. The domain boundaries are treated as reflecting solid walls. The initial pressure is given as ⎧ ⎨ 1000 −0.5 ≤ x ≤ −0.4 0.010 −0.4 ≤ x ≤ 0.4 (52) p(x, 0) = ⎩ 100 0.4 ≤ x ≤ 0.5 Note that each jump in the pressure is of the form of a shock-tube problem. The one near the left wall yields a strong shock and contact discontinuity that move into the intermediate low pressure gas. A left-moving rarefaction travels towards the wall. The rarefaction wave will be reflected from the left wall, catch up and weaken the contact discontinuity and shock wave. The shock-tube problem near the right wall has a similar wave pattern but the velocities are reversed. Owing to a milder pressure jump, the waves resulting from this problem are weaker. We set up this problem with 1000 equal mass particles placed uniformly in the domain (Δx = 0.001). The solid wall is implemented by reflecting the particles about either wall and negating their velocities. Since there is no analytical solution to this problem, we compare the SPH results with the Lagrange plus remap version of the Piecewise Parabolic Method (VH-1) [62] using 5000 grid cells and a CFL number of 0.1. Fig. 8 shows the numerical results for the van Leer exact, Ducowicz and Roe approximate Riemann solver at the time t = 0.016s. At this time, for the discontinuity near the left wall, the rarefaction has been reflected and has interacted with the right-moving contact discontinuity. A part of the rarefaction is reflected back and the part that is transmitted is about to catch the shock. The kink in the velocity profile in the FVM solution marks the edge of the rarefaction reflected from the contact discontinuity [61]. On the right, a similar structure develops but the interactions are yet to occur. From the results in Fig. 8 we see that all three schemes capture the behaviour well up until this point. The weakening of the contact discontinuity by the rarefaction has hidden some of the characteristic errors at this point (pressure “wiggling”). On closer examination, we find that the schemes produce a slight over shoot of the density at the contact discontinuity on the right. This is evident in the FVM solution as well and is perhaps an artifice of the Lagrangian nature of the two schemes. In [24], we added a small amount of thermal conduction (originally proposed by Noh [63]) to smear these spikes. This is not done here. The HLLE and LLXF Riemann solvers failed 27

7

Exact

Ducowicz

Roe

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

6 5

ρ

4 3 2 1 0

20 15

u

10 5 0 5 10

Figure 8: Numerical results (dots) for the density (upper panel) and velocity (lower panel), for the Woodward and Colella blast-wave problem at t = 0.016s, for GSPH with the van Leer exact (left), Ducowicz (centre) and Roe (right) approximate Riemann solvers with 1000 particles, compared with the reference PPM solution (solid line). The resolution of the flow structures for the and is comparable with the finite volume reference solution.

for this test when using the second order reconstruction. Fig. 9 show the numerical results at the time t = 0.028s. At this time, the two shocks resulting from the initial jumps have collided resulting in the sharp density spike and an additional contact discontinuity. The right going shock will continue and collide with the left moving contact discontinuity which resulted from the initial jump at x = 0.4. The peak density due to the collision, computed by the reference PPM solution is approximately 28. We see that the van Leer and Roe’s approximate Riemann solver is most accurate in capturing this peak. The Ducowicz solver has underestimated the density peak about 20%. Smearing of a discontinuity over a few smoothing lengths (by viscosity) is the reason why the schemes are unable to capture any sharp spike in the physical variables. We would have to increase the resolution to lower the resolvable length scales to get a better description of this spike. The HLLC approximate 28

30

Exact

Ducowicz

Roe

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

25

ρ

20 15 10 5 0

15

10

u

5

0

5

10

Figure 9: Numerical results (dots) for the density (upper panel) and velocity (lower panel), for the Woodward and Colella blast-wave problem at t = 0.028s, for GSPH with the van Leer exact (left), Ducowicz (centre) and Roe (right) approximate Riemann solvers using 1000 particles, compared with the PPM solution (solid line). The density spike, resulting from the shock collision is underestimated by all schemes.

Riemann solver fails just after this collision. As the flow evolves, the rarefaction catches up with this spike and lowers the value of the peak density [61]. The location of the density peak is found to be in slight error for each of the schemes. This error was not observed in our earlier comparisons of SPH schemes [24] (where we used a first order WAF scheme for the reference results). For the velocity, we find the results are identical for each of the schemes. Numerical scatter in the velocity occurs at exactly the location where there is noise in the PPM solution. As is customary for this problem, we examine the solution at the final time t = 0.038s. At this time, the solution has contact discontinuities at the locations x = 0.1, x = 0.26 and x = 0.3 and shock waves at x = 0.15 and x = 0.37. From the results in Fig. 10, we see that the SPH is in good agreement with the finite volume solution albeit with a slight error in the location of the shocks (x = 0.147 and x = 0.364). The error is of the order of one 29

7

Exact

Ducowicz

Roe

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

0.4 0.2 0.0 0.2 0.4

6 5

ρ

4 3 2 1 0

16 14 12

u

10 8 6 4 2 0 2

Figure 10: Numerical results (dots) for the density (upper panel) and velocity (lower panel), for the Woodward and Colella blast-wave problem at t = 0.038s, for GSPH with the van Leer exact (left), Ducowicz (centre) and Roe (right) approximate Riemann solvers using 1000 particles, compared with the PPM solution (solid line). The shock location is in slight error for each of the schemes. The rest of the flow structures are well captured by the GSPH schemes

smoothing length at this location. In general, we find that GSPH with the Ducowicz and Roe approximate Riemann solvers both produce results comparable with the original GSPH using the van Leer exact Riemann solver [21]. 5.6. Noh’s cylindrical implosion test Our first test in 2D is the cylindrical implosion test proposed by Noh [63]. The problem specifies a cold, ideal gas (γ = 53 ) in the X − Y Cartesian plane with a radial velocity of magnitude 1, directed towards the origin. The exact solution to the problem is known. The initial conditions result in an infinite strength, cylindrical shock that propagates outwards at the speed 3t . The post-shock density is a hot, stationary gas with density 16 while the  density just ahead of the shock front ( x2 + y 2 > 3t ) is given by ρ = 1 + √ 2t 2 . This x +y

30

Figure 11: Radial density profiles for the Noh’s cylindrical implosion test at t = 0.2s for GSPH with van Leer’s exact (left), Ducowicz (center) and Roe’s (right) approximate Riemann solver using 200×200 particles at t = 0.6s. Each of the schemes exhibit the spurious heating at the origin. A higher resolution can be used to reduce the errors in the post-shock region (Fig. 12)

is a difficult test, with numerical methods usually producing an erroneous “wall heating”, characterized by a dip in the density and a corresponding spike in the thermal energy at the origin. Lagrangian schemes are particularly susceptible to this behaviour and as such, this test has been used before ([64, 4, 65]) in the context of SPH. We solve the problem in Cartesian coordinates by distributing nx × nx particles in the quadrant [0, 1] × [0, 1], with symmetry boundary conditions along the planes x, y = 0. To avoid wrongly predicting cavitation, the cold gas is approximated as a gas with thermal energy e = 1 × 10−6 . The density is uniform (ρ = 1) and the particle velocities are directed radially inwards with magnitude 1. The solution is examined at the time t = 0.6s, at which time the shock is at the radial position r ≈ 0.2. Figure 11 shows the radial density profile using the exact (left), Ducowicz (center) and Roe’s (right) approximate Riemann solver using 31

nx = ny = 200 particles. Each of the schemes exhibits the spurious wall heating at the origin for this problem. Investigations into the nature of the spurious heating ([66, 67]) suggest entropy violations that occur near the origin during the transient shock reflection phase. Once the errors are generated, they remain frozen and are simply advected by the particles. Addition of an artificial thermal conduction, to dissipate this excess energy is the typical remedy suggested for this behaviour. We have not used any thermal conduction for the results presented herein. The results with SPH have a tendency to under-predict the post-

Figure 12: Radial density profiles for the Noh’s cylindrical implosion test at t = 0.2s using the GSPH with the van Leer exact Riemann solver at the resolutions of 50 × 50 (blue), 100 × 100 (green) and 200 × 200 (red) particles respectively. The error in the post-shock density reduces when using a higher resolution.

shock values. This is most evident for the solution using the Ducowicz approximate Riemann solver, where the density in the immediate post-shock region peaks, before dropping off to the minimum at the origin. There is also considerable scatter in the profile when compared to the Roe solver. The diffusive HLLE and LLXF solvers fail for this test by wrongly predicting negative values. The HLLC Riemann solver is qualitatively similar to the Roe solver, but with larger errors in the post-shock region. 32

The magnitude of the post-shock values can be made more accurate using a higher resolution. This is shown in Fig. 12, where we have used the exact Riemann solver with 50 × 50 (blue), 100 × 100 (green) and 200 × 200 (red) particles. At higher resolutions, the shock transition region is sharper as expected and the post-shock density approaches the exact value ρ = 16. 5.7. 2D Riemann problem A two-dimensional extension to the Riemann problem can be defined by requiring multiple states to meet at a corner. The resulting solution consists of a complex wave pattern that emerges from the individual jumps. Considering 4 constant states (one in each quadrant) meeting at a corner (origin) in 2 dimensions, Shulz-Rinne [68] categorized 15 unique configurations for the Euler equations, with the assumption of planar elementary waves connecting the states in adjacent quadrants. The elementary waves can be one of either a shock, a rarefaction or a contact discontinuity (slip line). The problems have been recognized as stringent tests for multi-dimensional codes [69, 70, 71, 55, 50] and we consider one (tough) case for SPH. Formally, the test corresponds to configuration 3 ([69, 70]), which involves four planar shocks interacting at a corner. The initial conditions define a constant state in 4 quadrants (NE, NW, SW, SE) with values shown in Table 5.7. This data results in four Variable ρ u v p

NE 1.5 0 0 1.5

NW 0.5323 1.206 0 0.3

SW 0.5323 1.206 1.206 0.029

SE 1.5 0 1.206 0.3

← − ← − ← − ← − shocks S 21 , S 32 , S 34 and S 41 , where the arrow indicates the direction of propagation and ← − the subscripts indicate to which quadrant they belong. For example, S 21 represents a shock that travels from the NE (quadrant 1) to the NW (quadrant 2) quadrant. The solution is symmetric about the line x = y, with a three shock configuration created by the intersection ← − ← − of S 21 and S 32 and a double Mach reflection that leads to a roll up of the slip line along the symmetry axis [69]. We set-up the problem using two kinds of discretizations. For the first, we use a constant particle spacing in the computational domain and adjust the particle mass so that m/Δx2 gives the required density for the particle in that quadrant. We use the same effective resolution of Δx = Δy = 1/400 as employed in [71, 55]. For the results, we focus on the lower left quadrant (SE) which has the most detail. Fig. 13 displays the interpolated density map in the 3rd , or SE quadrant for the PPM (VH-1) (a), van Leer exact (b), Ducowicz (c) and Roe approximate (d) Riemann solvers, at t = 0.2s. The PPM solution was generated using a resolution of Ncells = 400 in each coordinate direction. Qualitatively, each of the schemes are similar, showing good agreement with the PPM solution and symmetry about the line x = y. The slip line roll up at the centre is most prominent for the GSPH with the exact Riemann solver. This roll up is mitigated for the Ducowicz and Roe solutions by the presence of “spot” like structures along the density interfaces. On closer inspection, each of the 33

0.55

0.55

0.50

0.50

0.45

0.45

0.40

0.40

0.35

0.35

0.30

0.30

0.25

0.25

0.20 0.20

0.25

0.30

0.35

0.40

0.45

0.50

0.20 0.20

0.55

0.25

(a) VH-1

0.55

0.55

0.50

0.45

0.45

0.40

0.40

0.35

0.35

0.30

0.30

0.25

0.25

0.20 0.20

0.20 0.20

0.30

0.35

0.40

0.45

0.35

0.40

0.45

0.50

0.55

0.50

0.55

(b) Exact

0.50

0.25

0.30

0.50

0.55

(c) Ducowicz

0.25

0.30

0.35

0.40

0.45

(d) Roe

Figure 13: Density map for the 2D Riemann problem (table 5.7) for the PPM (VH-1), van Leer exact, Ducowicz and Roe approximate Riemann solver at t = 0.2s. The solution is shown in the third, or SE quadrant. The slip line roll up is most prominent for the GSPH with the van Leer exact Riemann solver. The use of constant mass particles has resulted in a “pairing” behaviour between particles of unequal mass.

schemes exhibits a type of pairing instability for this test, whereby groups of particles form regular, hexagonal structures about the symmetry plane x = y as can be seen in Fig. 14. The pattern is consistent with a ring of lighter particles encircling a central cluster of 3 to 4 heavier particles. We want to point out that this behaviour is not limited to GSPH or the use of the approximate Riemann solvers. For example, Fig. 15 shows the results obtained for this problem using the publicly available, NDSPMHD code [14], using nx = ny = 200 particles and the Quintic spline kernel. A similar clumping pattern is observable about the symmetry axis. The instability is therefore more fundamental to SPH rather than an artifice of the GSPH schemes discussed in this work. For the second kind of discretization, we use constant mass particles in the following manner. First, the minimum spacing Δx0 , corresponding to the highest density is determined so that m0 = ρ0 Δx20 for that quadrant. The spacing in the other quadrants is determined as Δx2q = ρ0 /ρq Δx20 , for q = N E, N W, SW, and SE respectively. We have used 34

Figure 14: Particle positions for the 2D Riemann problem using the Roe approximate Riemann solver, with unequal mass particles and at t = 0.2s. A pairing instability is visible whereby particles form regular hexagonal structures about the symmetry plane x = y

Δx0 = 1.75 × 10−3 , which results in a total of 72043 particles. Results using this discretization is shown in Fig. 16. The particle clumping behaviour is eliminated with the use of constant mass particles, but the solution detail is very poor for all the GSPH schemes. This is because the effective resolution in this quadrant is very low given the constraint of equal mass particles. Recall that the reference results were generated on a uniform grid of 400×400 cells, which corresponds to the cell size Δx = 2.5 × 10−3 . In contrast, the GSPH results were generated using an effective resolution (in terms of particle spacing) of Δx3 ≈ 11 × 10−3 in this quadrant. As a result, the slip line roll-up is barely visible for the GSPH schemes. Using a higher resolution yields better results but also leads to a prohibitively large number of particles. With our discretization procedure for example, an effective resolution of ≈ 3.3 × 10−3 in the third quadrant results in a prohibitively large number (≈ 1.5 × 106 ) of particles. Figure 17, shows the result for such a discretization using the Roe approximate Riemann solver. We see that the slip line roll-up is relatively well captured with this reso35

Figure 15: Particle positions for the 2D Riemann problem using the NDSPMHD code, with Np = 200 × 200 particles, the Quintic spline kernel and at time t = 0.15s. The particle clumping, similar to figure (14) is visible about the symmetry plane x = y

lution and we expect results similar to the reference FVM solution were we to use the same effective resolution in the SW quadrant for SPH. This however would result in a very large number of particles, rendering SPH impractical. This is ostensibly a drawback for the SPH method for solving this class of problems. Consequently, there is a need to develop a robust algorithm for SPH that can handle unequal mass particles. 6. Summary and conclusions Lagrangian Approximate Riemann solvers were used with the Godunov Smoothed Particle Hydrodynamics (GSPH) method and applied to test problems in one and two dimensions. The approximate solvers were introduced to improve the efficiency of GSPH schemes which requires the solution to the non-linear Riemann problem for every interacting particle pair. The numerical results in Sec. 5 suggest that the Roe, Ducowicz and HLLC approximate Rie36

0.55

0.55

0.50

0.50

0.45

0.45

0.40

0.40

0.35

0.35

0.30

0.30

0.25

0.25

0.20 0.20

0.25

0.30

0.35

0.40

0.45

0.50

0.20 0.20

0.55

0.25

0.30

(a) VH-1

0.40

0.45

0.50

0.55

0.50

0.55

(b) Exact

0.55

0.55

0.50

0.50

0.45

0.45

0.40

0.40

0.35

0.35

0.30

0.30

0.25

0.20 0.20

0.35

0.25

0.25

0.30

0.35

0.40

0.45

0.50

0.20 0.20

0.55

(c) Ducowicz

0.25

0.30

0.35

0.40

0.45

(d) Roe

Figure 16: Density map for the 2D Riemann problem (table 5.7) for the PPM (VH-1), van Leer exact, Ducowicz and Roe approximate Riemann solver at t = 0.2s. The solution is shown in the third, or SE quadrant. The vortex roll is most prominent for the GSPH with the van Leer exact Riemann solver.

mann solvers are viable alternatives to the exact Riemann solver for the class of problems considered. The use of these approximate Riemann solvers has the added benefit of applying GSPH to multi-material compressible hydrodynamics and/or with a general equation of state. This will be explored in a subsequent work. It is shown how the popular signal-based artificial viscosity terms for SPH can be derived in a consistent manner under a restriction of the class of approximate Riemann solvers. This places the former heuristic derivation ([20]) on a more formal footing and provides a rational justification for the use of artificial viscosity in SPH. The equivalence between the GSPH dissipation and the SPH artificial viscosity can be used to explain some of the features of the GSPH solution. In particular, GSPH has a non-zero viscosity that is useful in suppressing the spurious pressure jump for SPH but at the same time results in larger, dissipation-induced heating errors, as for the Sj¨ogreen’s test in Sec. 5.4. A general SPH clumping instability is observed for the solutions to the two-dimensional Riemann problem (Sec. 5.7) when using unequal mass particles. The constraint of constant mass particles results in an accurate solution but also leads to a prohibitively large num37

(a) VH-1

(b) Roe

Figure 17: Density map for the 2D Riemann problem (table 5.7) for GSPH with the Roe approximate Riemann solver (b) when using an effective resolution of Δx3 ≈ 3.3 × 10−3 in terms of particle spacing. The solution detail is more acceptable at this resolution with the slip line roll-up clearly visible. The drawback is that the constraint of equal mass particles results in a prohibitively large number (≈ 1.5 × 106 ) of particles.

ber of particles. Efficient solutions to such problems, comparable to reference finite volume methods requires a robust SPH algorithm that can handle unequal mass particles. This is left as an avenue for further investigation. ACKNOWLEDGEMENTS We would like to thank the anonymous reviewers for their comments and suggestions that have been very useful in improving the quality of this work. In particular, for the explanation of the clumping behaviour and it’s relation to the use of unequal mass particles.

38

References 1. R.A. Gingold and J.J. Monaghan, . Smoothed particle hydrodynamics: Theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society 1977;181:375–389. 2. J.J. Monaghan and R.A. Gingold, . Shock Simulation by the Particle Method SPH. Journal of Computational Physics 1983;52:374–389. 3. Monaghan, J.. Smoothed Particle Hydrodynamics. Annual Review of Astronomy and Astrophysics 1992;30:543–574. 4. Leonardo Di. Sigalotti and Hender Lopex and Arnaldo Donoso and Elroy Siraa and Jaime Klapp, . A shock capturing SPH scheme based on adaptive kernel estimation. Journal of Computational Physics 2006;212:124–149. 5. Agertz et al., . Fundamental differences between SPH and grid methods. Monthly Notices of the Royal Astronomical Society 2007;380:963–978. 6. G. Murante and S. Borgani and R. Brunino and S-H. Cha, . Hydrodynamic Simulations with the Godunov SPH. Monthly Notices of the Royal Astronomical Society 2011;417:136–153. 7. M. Steinmetz and E. M¨ uller, . On the capabilities and limits of smooth particle hydrodynamics. Astronomy and Astrophysics 1993;268:391–410. 8. Stephan Rosswog, . Astrophysical smooth particle hydrodynamis. New Astronomy Reviews 2009;53:78– 104. 9. Volker Springel, . Smoothed Particle Hydrodynamics in Astrophysics. Annual reviews of Astronomy and Astrophysics 2010;48:391–430. 10. S. Børve and M. Omang and J. Trulsen, . Regularized Smoothed Particle Hydrodynamics: A New Approach To Simulating Magnetohydrodynamic Shocks. The Astrophysical Journal 2001;561:82–93. 11. D. J. Price and J. J. Monaghan, . Smoothed Particle Magnetohydrodynamics - I. Algorithm and tests in one dimension. Monthly Notices of the Royal Astronomical Society 2004;348:123–138. 12. S. Børve and M. Omang and J. Trulsen, . Multidimensional MHD Shock Tests of Regularized Smoothed Particle Hydrodynamics. The Astrophysical Journal 2006;652:1306–1317. 13. Kazunari Iwasaki and Shu-ichiro Inutsuka, . Smoothed Particle Magnetohydrodynamics with a Riemann solver and the method of characteristics. Monthly Notices of the Royal Astronomical Society 2011;418:1668–1688. 14. Price, D.. Smoothed particle hydrodynamics and magnetohydrodynamics. Journal of Computational Physics 2012;231:759–794. 15. J. Reisner and J. Serencsa and S. Shkoller, . A space-time smooth artificial viscosity method for nonlinear conservation laws. Journal of Computational Physics 2013;235:912–933. 16. Alexander Kurganov and Yu Liu, . A new adaptive artificial viscosity method for hyperbolic systems of conservation laws. Journal of Computational Physics 2012;231:8114–8132. 17. I. G. Cameron, . An Analysis of the Errors Caused by Using Artificial Viscosity Terms to Represent Steady-State Shock-Waves. Journal of Computational Physics 1966;1:1–20. 18. John. W. White, . A New Form of Artificial Viscosity. Journal of Computational Physics 1973;11:573– 590. 19. Mark L. Wilkins, . Use of Artificial Viscosity in Multidimensional Fluid Dynamic Calculations. Journal of Computational Physics 1980;36:281–303. 20. Monaghan, J.. SPH and Riemann Solvers. Journal of Computational Physics 1997;136:298–307. 21. Shu-Ichiro Inutsuka, . Reformulation of Smoothed Particle Hydrodynamics with Riemann Solver. Journal of Computational Physics 2002;179:238–267. 22. S-H. Cha and Shu-ichiro Inutsuka and Sergei Nayakshin, . Kelvin-Helmholtz instabilities with Godunov SPH. Monthly Notices of the Royal Astronomical Society 2010;403:1165–1174. 23. S-H. Cha and A. P. Whitworth, . Implementations and tests of Godunov-type particle hydrodynamics. Monthly Notices of The Royal Astronomical Society 2003;340:73–90. 24. Kunal Puri and Prabhu Ramachandran, . A comparison of SPH schemes for the compressible Euler equations. Journal of Computational Physics 2014;256:308–333.

39

25. P. L. Roe, . Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics 1981;43:357–372. 26. Willian J. Rider, . A Review of Approximate Riemann Solvers with Godunov’s Method In Lagrangian Coordinates. Computers & Fluids 1994;23(2):397–493. 27. E. F. Toro, . Riemann solvers and numerical methods for fluid dynamics. Springer; 2009. 28. R. LeVeque, . Finite Volume Methods for Hyperbolic Problems. Cambridge University Press; 2002. 29. Chi-Wang Shu and Stanley Osher, . Efficient Implementation of Essentially Non-oscillatory ShockCapturing Schemes. Journal of Computational Physics 1987;77:439–471. 30. B. Einfeldt and C.D, Munz and P.L. Roe and B. Sj¨ ogreen, . On Godunov-Type Methods near Low Densities. Journal of Computational Physics 1991;92:273–295. 31. John K. Dukowicz, . A General, Non-Iterative Riemann Solver for Godunov’s Method. Journal of Computational Physics 1985;61:119–137. 32. Hong Luo and Joseph D. Baum and Rainald L¨ ohner, . On the computation of multi-material flows using ALE formulation. Journal of Computational Physics 2004;194:304–328. 33. Ball, G.J.. A Free-Lagrange method for unsteady compressible flow: simulation of a confined cylindrical blast wave. Shock Waves 1996;5:311–325. 34. Fedir V. Sirotkin and Jack J. Yoh, . A Smoothed Particle Hydrodynamics method with approximate Riemann solvers for simulating strong explosions. Computers & Fluids 2013;88:418–429. 35. J. P. Vila, . On Particle Weighted Methods and Smooth Particle Hydrodynamics. Mathematical Models and Methods in Applied Sciences 99;09. 36. D. Hietel and K. Steiner and J. Struckmeier, . A Finite-Volume Particle Method for Compressible Flows. Math Models Methods Appl Sci 2000;10:1363–1382. 37. Benidict D. Rogers and Robert A. Dalrymple and Peter K. Stansby, . Simulation of caisson breakwater movement using 2-D SPH. Journal of Hydraulic Research 2010;48:135–141. 38. Jean-Christophe Marongiu and Francis Leboeuf and Jo¨elle Caro and Etienne Parkinson, . Free surface flows simulations in Pelton turbines using and hybrid SPH-ALE method. Journal of Hydraulic Research 2010;48:40–49. 39. Phoevos K. Koukouvinis and John S. Anagnostopoulos and Dimitris E. Papantonis, . An improved MUSCL treatment for the SPH-ALE method: comparison with the standard SPH method for the jet impingement case. International Journal For Numerical Methods in Engineering 2013;71:1152–1177. 40. Monaghan, J.. Smoothed Particle Hydrodynamics. Reports on Progress in Physics 2005;68:1703–1759. 41. Dehnen, W., Aly, H.. Improving convergence in smoothed particle hydrodynamics simulations without pairing instability. Monthly Notices of the Royal Astronomical Society 2012;1082:1068–1082. 42. Bjorn Engquist and Stanley Osher, . One-Sided Difference Approximations for Nonlinear Conservation Laws. Mathematics of Computation 1981;36(154):–. 43. Ami Harten and Peter D. Lax and Bram van Leer, . On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws. SIAM Review 1983;25(1):35–61. 44. S. F. Davis, . Simplified Second-Order Godunov-Type Methods. SIAM Journal of Scientific and Statistical Computation 1988;9(3):445–473. 45. B. Van Leer, . Towards the Ultimate Conservative Difference Scheme. Journal of Computational Physics 1997;20:229–248. 46. Phillip Colella and Paul R. Woodward, . The Piecewise Parabolic Method (PPM) for gas-dynamical simulations. Journal of Computational Physics 1984;54:174–201. 47. Juan Cheng and Chi-Wang Shu, . A high order ENO conservative Lagrangian type scheme for the compressible Euler equations. Journal of Computational Physics 2014;257:143–168. 48. Bernd Einfeldt, . On Godunov-type Methods for Gas Dynamics. SIAM Journal on Numerical Analysis 1988;25(2):294–318. 49. E. F. Toro and M. Spurce and W. Speares, . Restoration of the Contact Surface in the HLL-Riemann Solver. Shock Waves 1994;4(22):25–34. 50. Dinshaw S. Balsara, . A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows. Journal of Computational Physics 2012;231:7476–7503.

40

51. Dinshaw S. Balsara and Michael Dumbser and Remi Abgrall, . Multidimensional HLLC Riemann solver for unstructured meshes - With appliction to Euler and MHD flows. Journal of Computational Physics 2012;261:172–208. 52. Bailey, D.A.. Lagrange-Remap Methods for the Euler Equations for Single and Multi Gas Flows. Tech. Rep.; The University of Reading; 2003. 53. J.P. Morris and J.J. Monaghan, . A switch to reduce SPH Viscosity. Journal of Computational Physics 1997;136:41–50. 54. Price, D.. Modelling Discontinuities and Kelvin-Helmholtz instabilities in SPH. Journal of Computational Physics 2008;227:10040–10057. 55. Richard Liska and Burton Wendroff, . Comparison of Several Difference Schemes on 1D and 2D Test Problems for the Euler Equations. SIAM Journal of Scientific and Statistical Computation 2003;25(3):995–1017. 56. Frederic A. Raiso, . Particle Methods in Astrophysical Fluid Dynamics. Progress of Theoretical Physics 2000;138:609–621. 57. Lee Cullen and Walter Dehnen, . Inviscid smoothed particle hydrodynamics. Monthly Notices of the Royal Astronomical Society 2010;408:669–683. 58. Guang-Shan Jiang and Chi-Wang Shu, . Efficient Implementation of Weighted ENO Schemes. Journal of Computational Physics 1996;126:202–228. 59. Tom Abel, . rpSPH: a novel Smoothed Particle Hydrodynamics Algorithm. Monthly Notices of the Royal Astronomical Society 2010;413:271–285. 60. Zhijun Shen and Wei Yan and Guixia Lv, . Behaviour of viscous solutions in Lagrangian formulation. Journal of Computational Physics 2010;229:4522–4533. 61. P. Colella and P. R. Woodward, . The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations. Journal of Computational Physics 1984;54:174–201. 62. John Blondin, . VH-1 The Virginia Numerical Bull Session ideal hydrodynamics PPMLR. 1990. URL: http://wonka.physics.ncsu.edu/pub/VH-1/index.html. 63. W. H. Noh, . Errors for Calculations of Strong Shocks Using an Artificial Viscosity and Artificial Heal Flux. Journal of Computational Physics 1978;72:78–120. 64. Leigh Brookshaw, . Smooth Particle Hydrodynamics in Cylindrical Coordinates. The Australian and New Zealand Industrial and Applied Mathematics Journal 2003;44:114–139. 65. M. Omang and S. Børve and J. Trulesn, . SPH in cylindrical coordinates. Journal of Computational Physics 2006;213:391–412. 66. Ralph Menikoff, . Errors When Shock Waves Interact Due to Numerical Shock Width. SIAM Journal of Scientific and Statistical Computation 1994;15(5):1227–1242. 67. Willian J. Rider, . Revisiting Wall Heating. Journal of Computational Physics 2000;162:395–410. 68. Carsten W. Shulz-Rinne, . Classification of the Riemann problem for two-dimensional gas dynamics. SIAM Journal on Mathematical Analysis 1993;24(1):76–88. 69. Carsten W. Shulz-Rinne and James P. Collins and Harland M. Glaz, . Numerical Solution of the Riemann Problem for Two Dimensional Gas Dynamics. SIAM Journal on Scientific Computing 1993;14(6):1394–1414. 70. Peter D. Lax and Xu-Dong Liu, . Solution of Two-Dimensional Riemann Problems of Gas Dynamics by Positive Schemes. SIAM Journal on Scientific Computing 1998;19(2):319–340. 71. Alexander Kurganov and Eitan Tadmor, . Solution of Two-Dimensional Riemann Problems for Gas Dynamics without Riemann Problem Solvers. Numer Methods Partial Differential Equations 2000;18:548– 608.

41