A smoothed particle hydrodynamics numerical scheme with a consistent diffusion term for the continuity equation

Accepted Manuscript

A smoothed particle hydrodynamics numerical scheme with a consistent diffusion term for the continuity equation Mashy D. Green, Renato Vacondio, Joaquim Peiro´ PII: DOI: Reference:

S0045-7930(18)30900-9 https://doi.org/10.1016/j.compfluid.2018.11.020 CAF 4062

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

31 March 2018 7 September 2018 21 November 2018

Please cite this article as: Mashy D. Green, Renato Vacondio, Joaquim Peiro, ´ A smoothed particle hydrodynamics numerical scheme with a consistent diffusion term for the continuity equation, Computers and Fluids (2018), doi: https://doi.org/10.1016/j.compfluid.2018.11.020

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.

ACCEPTED MANUSCRIPT

Highlights • Parameter-free numerical diffusive term based on a Riemann solver.

CR IP T

• Scheme guarantees consistency both inside the fluid and close to the free surface.

• The van Albada flux limiter is the best performer in standard SPH bench-

AC

CE

PT

ED

M

AN US

marks.

1

ACCEPTED MANUSCRIPT

A smoothed particle hydrodynamics numerical scheme with a consistent diffusion term for the continuity equation

CR IP T

Mashy D. Greena,∗, Renato Vacondiob , Joaquim Peir´ oa a Department b Department

of Aeronautics, Imperial College London, London SW7 2AZ, UK of Engineering and Architecture, Parma University, Parco Area delle Scienze 181/A, 43124 Parma, Italy

AN US

Abstract

A novel formulation for the diffusive term in the continuity equation is proposed to improve the stability of smoothed particle hydrodynamics weakly compressible scheme avoiding the introduction of empirical parameters. Densities at particle-particle interface have been computed by means of a first-order consis-

M

tent total variational diminishing reconstruction and a one-dimensional Roe’s approximate Riemann solver is applied to add the correct amount of diffusion. Results of numerical tests also demonstrate that the proposed method is able

ED

to guarantee consistency both inside the fluid and close to the free surface. Furthermore, a numerical analysis of several flux limiter functions has been conducted, finding that the choice of this function is a critical point to guaran-

PT

tee the accuracy of the method. It has been assessed, through the monitoring of the internal energy, that the van Albada limiter is more effective in dissipating

CE

spurious density fluctuations. Keywords: Smoothed particle hydrodynamics; Free-surface flow; Consistency; Riemann solvers; Total variation diminishing reconstruction; Artificial

AC

dissipation.

∗ Corresponding

author Email address: [email protected] (Mashy D. Green)

Preprint submitted to Elsevier

November 21, 2018

ACCEPTED MANUSCRIPT

1. Introduction The smoothed particle hydrodynamics (SPH) method is a Lagrangian meshless numerical scheme originally introduced by Gingold and Monaghan [1] and by

CR IP T

Lucy [2] to solve astrophysical problems. Since then the method has progressed

significantly as described, for instance, in the various reviews by Monaghan

[3, 4, 5, 6] and the textbook by Viouleau [7]. Several variants of the original approach have been proposed [8, 9] enabling the simulation of different engineer-

ing applications which were considered too complex for classical computational fluid dynamics (CFD) Eulerian approaches [10, 11, 12, 13].

AN US

The most common SPH numerical schemes are based on the so-called weakly compressible approach. The SPH spatial discretization via a kernel does not account for the convective nature of the governing equations and, as a result, numerical simulations exhibit spurious numerical oscillations. To damp these oscillations and achieve stability, a diffusive term is often added to the momen-

M

tum equation [4]. This additional term, often referred to as artificial viscosity, incorporates a tuning parameter which requires calibration. Such an approach allows to reproduce the kinematic of violent flows with fragmented free-surfaces.

ED

On the other hand, it is also well known that the standard SPH scheme with artificial viscosity is affected by spurious oscillations in the pressure field [14].

PT

To address this issue, Molteni and Colagrossi [15] introduced an additional diffusive term in the continuity equation which effectively improved the accuracy of the pressure field, however it was over-dissipative and could not preserve the

CE

hydrostatic properties of the fluid. Antuono et al. [16] proposed a variant of the diffusive term introduced in [15] which is consistent at the free surface and thus reproduces more accurately the hydrostatic components of the pressure. The

AC

same formulation was adopted also by Marrone et al. [17] to simulate violent impact flows and in [18, 19, 20] for flow past different bodies with and without free surface. It is worth pointing out that these SPH methods are not well-balanced in the sense of the C-property discussed in [21] and therefore cannot guarantee an exact balance, within the precision of computer arithmetic, between the flux

3

ACCEPTED MANUSCRIPT

and gravity terms. In all the previously mentioned works the diffusion term in the continuity equation is controlled by the so-called δ empirical parameter which is unique for all particles and it has to be adjusted by the user. Usually

CR IP T

the δ parameter is assumed equal to 0.1, but there is little theoretical justification for this choice, beyond its general range based on linear stability analysis [22], and seldom a sensitivity analysis to different values of δ is presented. An

alternative approach for stabilizing the SPH scheme without the addition of an

artificial viscosity borrows ideas from finite-volume discretizations of hyperbolic equations where the numerical fluxes on the control volume are evaluated us-

AN US

ing an exact or approximate Riemann solver. This is often referred to as the

Godunov SPH (GSPH) approach initially introduced in [23, 8] where the Riemann solver is applied to a virtual interface between two particles. A number of reconstruction-based SPH schemes have been proposed in the literature. For instance, references [8, 24, 25, 26] present second-order accurate versions whilst reconstructions of higher order have been proposed in references [27, 28]. In

M

simulating free-surface flows, some GSPH formulations could be too dissipative and thus inaccurate when the free-surface deformation is large.

Ferrari et al.

ED

[29] applied an approximate Riemann solver only in the continuity equation, reducing the amount of diffusion of the original GSPH approach and avoiding the introduction of any empirical parameter. However, similarly to Molteni & Co-

PT

lagrossi [15], the formulation diverges at the free surface and thus it is not able to preserve the hydrostatic solution. A comprehensive comparison of different

CE

diffusion terms in the continuity equation has been conducted in [30]. We propose to re-interpret the formulation of Antuono et al. [16] as an

approximate Riemann solver with first-order reconstruction of the density at the

AC

particle-particle interface. This leads to a novel formulation of the diffusion term in the continuity equation which is consistent at the free surface and does not require any empirical parameter that has to be tuned. The scheme is obtained by employing a total variation diminishing (TVD) reconstruction and adopting Roe’s approximate Riemann solver aimed at introducing the correct amount of diffusion. Similarly to the scheme proposed by Antuono et al. [16], the herein 4

ACCEPTED MANUSCRIPT

proposed formulation guarantees the preservation of the hydrostatic pressure component. In addition, we have conducted an assessment of several flux limiter functions since the amount of added numerical dissipation depends on the choice

CR IP T

of limiter and this is a critical aspect in ensuring that the method is consistent at the free surface and that the volume and energy are accurately conserved by the method.

The contents of the paper are organized as follows. Section 2 introduces the

governing equations, their spatial and temporal discretisation and a detailed description on the continuity dissipation terms. This is followed by the intro-

AN US

duction of the Godunov SPH (GSPH) formulation for the continuity equation, the re-interpretation of the formulation in Antuono et al. [16] in the GPSH framework and lastly the new approximate Riemann solver formulation for the continuity equation. Section 3 presents several verification and validation benchmark test cases where the different dissipation schemes are numerically assessed

future work.

ED

2. Formulation

M

and compared. Finally, Section 4 presents our conclusions and suggestions for

We introduce the governing equations and the basic SPH formulation first,

PT

followed by a description of the δ-SPH framework, the Godunov formulation of the continuity equation, and a discussion of their equivalence.

CE

2.1. Governing equations The Lagrangian form of the compressible isentropic Navier-Stokes equations

AC

is

Dρ Dt Dv Dt Dx Dt

=

−ρ∇ · v = M

=

1 ∇σ + fe = F ρ

=

v

5

(1)

ACCEPTED MANUSCRIPT

where ρ is the density, x is the position vector, v is the velocity vector, and fe is an external force per unit mass, such as gravity.

CR IP T

The stress tensor σ is given by  
2 σ = −pI + µ − (∇ · v) I + ∇v + ∇vT 3

and the pressure, p, is calculated using a suitable equation of state to close the system. Here we use p(ρ) = c20 (ρ − ρ0 )

(2)

where ρ0 is a reference density and c0 is the reference speed of sound, chosen to

AN US

correspond to a nominal low Mach number, for instance Ma< 0.1. 2.2. Basic SPH

A SPH discretization assumes that fluid variables are associated with a set of N particles with coordinates xa ; a = 1, . . . , N [7]. The value of a variable β

M

at a particle a, denoted by βa , is approximated by using a kernel interpolant as βa =

X b

Vb βb W (kxa − xb k, h)

(3)

ED

where W is the SPH kernel function and the summation extends over the set of neighbour particles, b, within a distance h of the particle a, which has an

PT

associated volume Va . Introducing this spatial discretization in the governing equations (1) and the equation of state (2) leads to the semi-discrete SPH for-

AC

CE

mulation

Dρa Dt Dva Dt Dxa Dt pa

=

−ρa ∇ · va = Ma

(4)

=

1 ∇(−pa ) + fe = Fa ρa

(5)

=

va

(6)

=

c20 (ρa − ρ0 )

(7)

The terms of the right-hand side of the continuity and the momentum equa-

6

ACCEPTED MANUSCRIPT

tions can be expressed in conservative form, see e.g. [7], as X

Vb (va − vb ) · ∇a Wab + Da

(8)

1 X Vb (pa + pb ) ∇a Wab + Πa + fe ρa

(9)

= ρa

b

Fa

= −

b

CR IP T

Ma

where Vb = mb /ρb is the volume associated to particle b, Wab = W (kxa − xb k, h) is the SPH kernel, and ∇a Wab is the spatial gradient of kernel evaluated at the

location of particle a. For the simulations described here we have selected the Wendland kernel [31] given by q αd (2q + 1)(1 − )4 ; hd 2

and its derivative ∇W (x − x0 , h) = where q =

1 h

0≤q≤2

AN US

W (x − x0 , h) =

   αd q 3 x − x 0 ; −5q 1 − hd+1 2 kx − x0 k

0≤q≤2

(10)

(11)

kx − x0 k, d is the number of dimensions of the problem, and 7 4π

M

the values of the coefficient αd are

and

21 16π

in two and three dimensions,

respectively. The variables Da and Πa represent diffusion terms. The term

ED

Da is a numerical dissipation added to the continuity equation which will be discussed in detail in Section 2.3. The diffusion term Πa can act as either the physical viscosity, or as a numerical dissipation commonly referred to as

PT

artificial viscosity. We consider only inviscid flows and interpret it exclusively as an artificial viscosity. The expression for the term Πa that we have adopted

AC

CE

is that of reference [3], which reads X αhc0 µab   ∇a Wab , (va − vb ) · (xa − xb ) < 0 2mb  (ρa + ρb ) b Πa =   0, (va − vb ) · (xa − xb ) ≥ 0

(12)

where α is a problem-dependent tuning parameter and µab is defined as µab =

(va − vb ) · (xa − xb ) kxa − xb k2

(13)

The solid boundaries are described using the fixed ghost particle formulation described in reference [32]. In this description, a number of layers, N , of fixed 7

ACCEPTED MANUSCRIPT

ghost particles are placed on the outer side of a solid boundary to assure that fluid particles always have a full kernel support. We exclusively use three layers of ghost particles and therefor fix the smoothing length of the kernel to h =

CR IP T

1.5∆x, where ∆x is the particle spacing. To obtain a continuous flow field, the ghost particles are assigned pressure values obtained by a force balance on either side of the solid wall as P P f pf Wwf + f ρf (g + as ) · (xw − xf ) Wwf P pw = f Wwf

(14)

where the subscripts f and w denote the fluid and the ghost particles, respec-

AN US

tively, and as is the specified wall acceleration, i.e. as = 0 when the wall is fixed. The ghost particles have a constant mass, taken to be equal to the initial mass of a fluid particle, and density is obtained from the equation of state (7) as ρw = ρ0 +

pw c20

(15)

M

This treatment ensures that the fluid pressure field adjacent to the solid walls is smooth. Assuming only a free-slip boundary condition, artificial viscosity (equation (12)) is switched off in the fluid-ghost particle interactions so that no

ED

velocity terms appear in the momentum conservation equation (9). The velocity of the ghost particles in the continuity equation is then set to the specified wall

PT

velocity

vw = vs

(16)

noting that vs = 0 when the wall is fixed.

CE

Lastly, the semi-discrete system (4–7) is integrated in time using a symplectic

AC

predictor-corrector time-stepping scheme [33] given by n+ 21

ρa

n+ 12

va

n+ 21

xa

ρn+1 a

∆t n M 2 a ∆t n = van + F 2 a ∆t n = xna + v 2 a    n+ 12 2− Ma = ρna ;  = −∆t 2+ ρa = ρna +

8

(17)

ACCEPTED MANUSCRIPT

n+ 12

van+1 = van + ∆tFa xn+1 = xna + a


∆t n va + van+1 2

To ensure stability, the time-step is calculated as ∆t = c∆t min(∆tf , ∆tcv )

CR IP T

where the pressure is calculated with the equation of state (7) at each half-step.

(18)

where c∆t is a coefficient in the range 0.1 ≤ c∆t ≤ 0.3, and a

∆tcv = min a

p

h/kFa k,

AN US

∆tf = min

h h(va − vb ) · (xa − xb ) c0 + maxb kxab k2

Here ∆tf is based on the force per unit mass and ∆tcv combines the convective and viscous time-step restrictions.

M

2.3. Density dissipation in the continuity equation

According to the δ-SPH formulation proposed in reference [22], the dissipa-

ED

tive term, Da , in the continuity equation (8) can be written as Da = −2c0 h

X b

ψba

(xb − xa ) · ∇a Wab Vb k(xa − xb )k2

(19)

PT

where ψba has the dimension of density and can vary according to different formulations. For instance, the term proposed by Molteni et al. [15] is equivalent

CE

to defining ψba as

ψba = δ (ρb − ρa )

(20)

AC

where δ is a tuning coefficient, usually taken to be δ = 0.1. This term has proven to be effective when dealing with short-duration advection-dominated flows. However, reference [22] has shown that this term is over-dissipative and does not preserve the hydrostatic component of the pressure or conserves the fluid volume.

9

ACCEPTED MANUSCRIPT

To address these shortcomings, Antuono et al. [16] proposed an alternative

(21)

CR IP T

expression for the term ψba as  
1 C C ψba = δ (ρb − ρa ) − ∇ρa + ∇ρb · (xb − xa ) 2

C where ∇ρC a and ∇ρb are the corrected density gradients for particles a and

b, respectively, and the same value of δ is typically used. In [16] the density

gradients in equation (21) are obtained by a SPH-corrected interpolation [34] ∇ρC a = with La =

X b

(ρb − ρa )La ∇a Wab

AN US

as

X −1 (xb − xa ) ⊗ ∇a Wab b

(22)

(23)

where the symbol ⊗ denotes a vector cross product.

2.4. Godunov SPH formulation of the continuity equation

M

An alternative form of writing the semi-discrete SPH equations is given by the Godunov SPH method (GSPH). In this method, the continuity equation is X Dρa ¯ ab ) · ∇a Wab Vb = 2ρa (va − v Dt

ED

written as [25]

(24)

b

¯ ab is the velocity obtained from the approximate Riemann solver. Knowwhere v

CE

PT

ing that ∇a Wab can be rewritten as ∇a Wab =

1 xa − xb ∂W h kxa − xb k ∂q

(25)

AC

we can rewrite equation (24) as Dρ ρa X p ∂W p = −2 (va − v¯ab ) Vb Dt h ∂q

(26)

b

where vap = va ·

xb −xa kxb −xa k

p ¯ ab · and v¯ab =v

xb −xa kxb −xa k .

p Adopting Roe’s Riemann solver, the velocity v¯ab can be calculated as [24] p v¯ab =

2 (prb − pra ) 1 p (va + vbp ) − 2 c0 (ρra + ρrb ) 10

(27)

ACCEPTED MANUSCRIPT

where the superscript r denotes the reconstructed value of the variables (pressure, density and velocity) at the interface between particles a and b. In this work the reconstructed values are computed by means of a TVD reconstruction, 1 C βar = βa + ξ(∇βab ) · (xb − xa ) 2

CR IP T

which reads (28)

where β is an arbitrary value to be reconstructed (i.e. the density, pressure or C velocity component), ξ is the limiter function [35] and ∇βab =

C ∇βa ·(xa −xb ) . βb −βa

Substituting the pressures with the density following the equation of state (7), equation (27) can be rewritten as 1 p 2 (ρr − ρra ) (va + vbp ) − c0 rb 2 (ρa + ρrb )

AN US

p v¯ab =

(29)

Finally, substituting equation (29) in equation (26) we can write the continuity equation as

Dρ ρa X p ∂W 2ρa c0 X 2 (ρrb − ρra ) ∂W =− (va − vbp ) Vb − Vb Dt h ∂q h (ρra + ρrb ) ∂q

(30)

b

M

b

2.5. Antuono δ-SPH formulation in a Godunov SPH framework

ED

We will demonstrate that the δ-SPH formulation of equation (19) is equivalent to the adoption of Roe’s Riemann solver in the continuity equation. Substituting equation (25) in the continuity equation (8) and in equation(19)

PT

of the δ-SPH formulation, we can rewrite the continuity equation as X Dρa ρa X p ψba ∂W ∂W =− (va − vbp ) Vb − 2c0 Vb Dt h ∂q kxa − xb k ∂q

CE

b

(31)

b

Realising that the first terms on the RHS of equations (30) and (31) are

AC

identical allows us to write ψba as ψba =

kxa − xb k 2 (ρrb − ρra ) ρa r h (ρa + ρrb )

(32)

which makes equation (30) equivalent to equation (31). This means that the formulation presented in Antuono et al. [16] might be seen as a particular type of Riemann solver, applied only in the continuity equation.

11

ACCEPTED MANUSCRIPT

Furthermore, taking reconstructed values of the density using equation (28)

where Bab is defined as

Bab =

kxa − xb k 2 ρa h (ρra + ρrb )

(33)

CR IP T

and substituting it in equation (32), we can now write  
1 C C ψba = Bab (ρb − ρa ) − ξab ∇ρa + ξba ∇ρb (xb − xa ) 2

(34)

The structure of equations (21) and (33) is very similar, the only differences are the presence of the limiter functions ξab and ξba used in the TVD reconstruction,

AN US

the δ coefficient which does not appear in equation (33), and its substitution by the Bab term.

With the aim of further demonstrating that the formulation of Antuono et p al. [16] is equivalent to a particular Riemann Solver the v¯ab velocity (equation

27) can be re-defined as

1 p h 2 (prb − pra ) (va + vbp ) − δ 2 kxa − xb k c0 (ρra + ρrb )

M

p v¯ab =

(35)

Comparing equation (35) with equation (27) it can be noticed that the new

ED

p v¯ab expression contains the coefficient

h kxa −xb k

which varies between 0.5 (when

kxa − xb k = 2h) to ∞ (when kxa − xb k = 0+ ) and the δ coefficient. Substituting equation (35) in equation (26) and using the equation of state

PT

(7), the continuity equation can be expressed as

CE

X ρa X p 2 (ρrb − ρra ) ∂W Dρ ∂W 1 =− Vb − δ2ρa c0 Vb (36) (va − vbp ) Dt h ∂q kxa − xb k (ρra + ρrb ) ∂q b

b

Imposing that equation (36) has to be identical to equation (31) and substitut-

AC

ing equation (28) for the reconstructed densities, the expression of ψab can be obtained as  
1 L ψba = Cab δ (ρb − ρa ) − ξab ∇ρL + ξ ∇ρ (x − x ) ba b a a b 2

where Cab is defined as

Cab =

2 ρa (ρra + ρrb ) 12

(37)

ACCEPTED MANUSCRIPT

Please note that, at least for a weakly-compressible flow Cab ∼ = 1 and so the only differences between the Riemann solver obtained from equation (35) and

CR IP T

the formulation of Antuono et al. [16], i.e. equation (21), are the limiters.

2.6. Time integration analysis

We have conducted a dimensional analysis to identify if an additional time-

step restriction has to be introduced to ensure the stability of the numerical scheme of equation (30). We have shown in the previous section that our GSPH

AN US

formulation (30) is equivalent to the δ-SPH formulation of Antuono et al. [16] if

Bab replaces δ in equation (33). Starting from this and following the procedure proposed by Antuono et al. [22], the stability restriction on the time-step is   h ∆tD < C∆t (38) Bab c0 Substituting equation (34) in equation (38) and recalling that the flow is as2 ρa

M

sumed to be weakly compressible, thus

ED

∆tD < C∆t



(ρra +ρrb )

' 1, we obtain

h2 kxa − xb kc0



PT

Knowing that kxa − xb k ≤ ξh (here ξ = 2), it follows that   h ∆tD < C∆t ξc0

(39)

(40)

which is equivalent to the customary CFL condition (18) for explicit SPH

CE

schemes. This means that no additional time-step restriction has to be adopted for the present algorithm.

AC

3. Verification and validation of the proposed formulation Three validation and verification tests are discussed to establish the agree-

ment between the δ-SPH scheme and the Riemann solver based scheme proposed. In the following we will assess the various dissipative schemes in the continuity equations defined as follows:

13

ACCEPTED MANUSCRIPT

• The Antuono formulation of the δ-SPH scheme [16] (ADSPH) as described in equations (8), (19) and (21); • the Molteni density dissipation scheme [15] (MDSPH) as described in This scheme is also referred to as the

CR IP T

equations (8), (19) and (20).

Molteni δ-SPH scheme for convenience due to fitting within the same general framework;

• the full Riemann solver formulation of the continuity equation (FRS) as defined in equation (30);

AN US

• the Riemann solver scheme as described in equations (36) with limiter function assumed always equal to 1 (ARS), which should be equivalent to the formulation of Antuono et al. [16].

The first test is a hydrostatic long duration test where the equivalence of the ADSPH and ARS schemes is validated, and continues with studying the effects

M

of various limiters in the FRS scheme on energy conservation. The second test validates the scheme in the dynamic case of a jet impact, and a convergence study is conducted. The final test is a verification benchmark test of a 2D dam

ED

break where the FRS scheme is compared against the MDSPH scheme. 3.1. Hydrostatic Dissipation

PT

A fluid rectangular domain of width L = 1m and height H = 0.5m is defined for this case with a particle spacing of ∆x = 0.005m. The fluid is assumed to be

CE

inviscid with density ρ0 = 1 000kg/m3 , resulting in 20 000 fluid particles. The p speed of sound is chosen as c0 = 20 |gH| with g = −9.81m/s2 the acceleration

due to gravity in the z-axis direction. No artificial viscosity is used, i.e. α = 0 in

AC

equation (12), unless indicated otherwise. The parameter in the δ-SPH schemes is chosen as δ = 0.1 wherever applicable. Periodic boundaries in the x-axis direction are employed to reduce influences from the choice of solid boundary treatment. The duration of each of the simulations is set to 50s. This first hydrostatic case compares the different dissipation schemes. These are the ADSPH, ARS and FRS. The FRS scheme is using the van Albada limiter 14

ACCEPTED MANUSCRIPT

[36] for reasons explained in Section 3.2. For reference, a case using only the artificial viscosity (AV), as defined in equation (12), with a value of α = 0.01 is used. Smaller values of α were numerically determined to be unstable under this

CR IP T

configuration. This test verifies the equivalence of the formulation of the ARS scheme to the standard ADSPH formulation, and ensures that the full Riemann

solver formulation is comparable. In the following, both the pressure field and

energy conservation are examined. The method for calculating the total energy and its components is described in Appendix A.

The three schemes that provide dissipation to the continuity equation are all

AN US

similarly capable of obtaining an accurate pressure profile at the end of the 50s

simulation, as seen in Figure 1, while the case using only artificial dissipation contains some noise. The reason for this becomes clear after examining the energy conservation. Despite the total energy conservation of the AV scheme reaching a steady level, shown in Figure 2(a), the internal energy conservation, shown in Figure 2(d), shows it continues to have strong oscillations throughout

M

the simulation, which are cancelled out in the total energy by the oscillations in the potential energy, shown in Figure 2(b). It is also worth noting that the

ED

kinetic energy during the first instances of the simulation is much lower. While this might seem like a beneficial feature, it actually restricts the particles from rearranging themselves into a low energy state. Focusing on the other schemes,

PT

the ADSPH and ARS formulations are indeed extremely similar, validating their equivalence. The FRS formulation shows a large initial drop in potential energy

CE

(and hence total energy), but stays roughly constant after the initial drop whilst being effective in dissipating internal energy and not over-dissipating kinetic energy, which are the two desired features for an effective numerical dissipation

AC

scheme in SPH. 3.2. Assessment of Limiter’s Performance in Hydrostatic Simulation This test consists of the same properties described in Sec. 3.1. Here, a

comparison between the FRS scheme using a variety of limiters is conducted. The limiters chosen are all TVD limiters, ranging between minmod (ξmm ) (the 15

1

0.5

0.5

0 0

0.2

0.4

0.6

0.8

(a)

1

0 0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

(b)

1

M

1

0.4

0.6

0.8

1

0.5

0 0

0.2

PT

0.2

ED

0.5

0 0

AN US

1

CR IP T

ACCEPTED MANUSCRIPT

(c)

0.4

(d)

CE

Figure 1: Normalised pressure of all fluid particles against their height at t = 50s, compared to

the analytic hydrostatic pressure profile: (a) The full Riemann solver formulation (FRS); (b) the Antuono scheme written in terms of a Riemann solver (ARS); (c) the standard Antuono

AC

scheme (ADSPH); (d) the Monaghan artificial viscosity scheme (AV).

16

10-5

5

0 -0.5

-5

-1

-10

-1.5

-15 0

10

20

30

40

(a) 10

5

1

20

PT

10

ED

3 2

10

20

30

40

50

40

50

(b)

-5

4

0 0

-2 0

50

M

5

10-4

AN US

0

CR IP T

ACCEPTED MANUSCRIPT

30

10

-5

4 3 2 1

40

50

10

20

30

(d)

CE

(c)

0 0

Figure 2: Temporal evolution of the normalised energy components of the different dissipation

AC

schemes: (a) total energy; (b) potential energy; (c) kinetic energy; (d) internal energy.

17

ACCEPTED MANUSCRIPT

most dissipative) to superbee (ξsb ) (the least dissipative) [37], which include the van Albada (ξva ) [36] and the van Leer (ξvl ) [38] limiters. The expressions of these limiters are, respectively, given by =

max[0, min(1, β)];

ξsb (β)

=

max[0, min(2β, 1), min(β, 2)];

ξva (β)

=

ξvl (β)

=

β2 + β ; β2 + 1 β + |β| ; 1 + |β|

lim ξmm (β) = 1

β→∞

lim ξva (β) = 1

β→∞

lim ξvl (β) = 2

β→∞

(41)

CR IP T

ξmm (β)

lim ξsb (β) = 2

β→∞

(42)

(43)

(44)

AN US

For reference the ARS scheme is included in the analysis too.

Figure 3 shows the pressure of each of the particles against its height compared to the analytic hydrostatic profile. All the limiters provide a good hydrostatic profile close to the free surface. However, the midmod limiter shows a noticeable drop in pressure towards the bottom of the tank, and the superbee

M

shows a significant increase. Only the van Albada and van Leer limiters provide an accurate hydrostatic profile throughout the domain. The same behavior is observed in the conservation of energy in Figure 4, where the potential (and

ED

hence total) energy shows the minmod limiter to be over-dissipative and the superbee under-dissipative. Focusing on the more accurate limiters in Figure 5 there are noticeable differences between the van Leer and van Albada limiters.

PT

The kinetic energy is very similar, however, the potential energy shows a large initial drop in the van Albada case, while the van Leer shows an increase. Over-

CE

all, the case using the van Leer limiter is closer to zero, while the case with the van Albada limiter is slightly less oscillatory after the initial drop. However, looking the the internal energy, it is clear that the van Albada limiter has su-

AC

perior capabilities of dissipating the spurious density fluctuations, which is the key point of the continuity equation artificial dissipation schemes. Due to this reason, the van Albada limiter is the limiter of choice for all further investigations in this work and is the strong recommendation of the authors for the use in the proposed scheme.

18

1

0.5

0.5

0 0

0.2

0.4

0.6

0.8

(a)

1

0 0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

(b)

1

1.5

M

0.8 0.6 0.4

1

0.4

0.6

0.8

1

PT

0.2

ED

0.5

0.2 0 0

AN US

1

CR IP T

ACCEPTED MANUSCRIPT

(c)

0 0

0.2

0.4

(d)

CE

Figure 3: Normalised pressure of all fluid particles against their height at t = 50s using the full Riemann solver formulation (FRS), compared to the analytic hydrostatic pressure profile: (a) using the van Albada limiter; (b) using the van Leer limiter; (c) using the minmod limiter;

AC

(d) using the superbee limiter.

19

ACCEPTED MANUSCRIPT

10-3

8 6 4

0 -2 -4 0

5

10

15

20

25

30

35

(a) 8 6 4 2 0 -2 -4 0

10

15

20

25

45

50

30

35

40

45

50

30

35

40

45

50

30

35

40

45

50

M

5

40

AN US

10-3

CR IP T

2

(b)

5

10-5

ED

4 3

PT

2 1

0 0

10

15

20

CE

5

AC

10

25

(c)

10-5

5

0

-5 0

5

10

15

20

25

(d) 20 Figure 4: Temporal evolution of the normalised energy components of the full Riemann solver formulation (FRS) with all the different limiters tested against the Antuono scheme given in terms of a Riemann solver (ARS): (a) total energy; (b) potential energy; (c) kinetic energy; (d) internal energy.

10-4

10-4 0.5

0.5

0

0

-0.5

-0.5

-1

-1

AN US

1

-1.5

-1.5 0

10

20

30

40

(a)

M

5

4

2 1

20

10

20

30

40

50

30

40

50

10-5

4 3 2 1

40

50

0 0

10

20

30

PT

10

ED

3

0 0

-2 0

50

(b)

10-5

5

CR IP T

ACCEPTED MANUSCRIPT

Figure 5: Temporal evolution of the normalised energy components of the full Riemann solver

CE

formulation (FRS) with only the van Leer and van Albada limiters against the Antuono scheme given in terms of a Riemann solver: (a) total energy; (b) potential energy; (c) kinetic energy;

AC

(d) internal energy.

21

ACCEPTED MANUSCRIPT

3.3. Verification: Jet Impact To verify the formulation for a dynamic test case, a two-dimensional jet impact case is studied.

CR IP T

The problem consists of a column of water impacting perpendicular to a

flat plate and has a steady-state analytical solution [39]. The component of the velocity on the plate, u, normalized by the jet velocity, U , and the corresponding

AN US

pressure, P , can be obtained by solving the system x 1 1 + q + sin−1 q = ln H 2 1−q p 1 − q2 − 1 u = q 1 2 P = ρU (1 − u2 ) 2

(45)

where x is a coordinate on the plate, H is the width of the water column and q is an auxiliary variable. Due to the lack of an implementation of inflow and outflow boundary conditions in the current version of the code, the simulation is

M

conducted as follows. The initial height of the column of water is chosen so that it will not fully collapse on the plate at the end of the simulation, and the plate

ED

is set to be wide enough so that the edges have minimal influence on the fluid domain in the vicinity of the impact. The dimension, L, of the computational domain is taken to be H/L = 0.4. The particle spacing is set to H/∆x = 80 and

PT

the velocity, u, is normalised by the inflow velocity U . The numerical speed of sound is set to c0 = 30U and the artificial viscosity parameter is set to α = 0.1.

CE

The simulation is run for a duration of tU/L = 10. As it reaches the steady state only at tU/L = 2, the time interval tU/L = [0, 2] can be neglected in the analysis. We monitor the time history of the values of the pressure at the

AC

stagnation point P0 = [0, 0] and plot the pressure along the plate in the interval

x = [−0.8L, 0] at a given time. The FRS scheme using the van Albada limiter

[36] is compared against the MDSPH, ADSPH and the ARS schemes. The different dissipation schemes all seem to be similarly capable of matching

the analytical solution whilst avoiding any significant spurious oscillations in the pressure field. The ADSPH and ARS cases both show slightly higher oscillations 22

ACCEPTED MANUSCRIPT

(both temporal and spatial) than the MDSPH and FRS cases. Further, the FRS case shows a slightly higher value at P0 while the MDSPH case shows a slightly lower one. However, these changes are all within 2% of the analytical values

CR IP T

and are likely related to the use of artificial viscosity and the implementation of the boundary conditions as much as to the dissipation schemes themselves.

Finally, in order to verify the FRS scheme using the van Albada limiter, a convergence study with varying particle spacing of H/∆x = 10 · 2n for n = 1, 2, 3, 4, and using a much lower artificial viscosity value of α = 0.01 (in order to reduce the impact the α parameter on the varying resolution) is

AN US

provided. All the test use the same fixed ratio of h = 1.5∆x and three layers

of ghost particles. Figure 7 shows the pressure field for the highest resolution adopted (H/∆x = 160) where no spurious pressure oscillations can be appreciated. The simulations clearly converge towards the analytic solution as the resolution increases, however due to the nature of the SPH method it is difficult to obtain a clear picture of the rate of convergence. Nevertheless, the following

M

material represents our best attempts to conduct such a study. Figure 8(a) shows the pressure along the plate at tU/L = 10, and Figure 8(b) shows accu-

ED

mulated error along with the calculated order of convergence (COC), which is slightly better than a linear convergence of order one. The accumulated error was evaluated by integrating (over space) the L2 norm of the difference be-

PT

tween the analytical and computed pressures. To better understand the rate of convergence of the solution over time, the pressure at P0 over the interval

CE

tU/L = [8, 10] is provided in Figure 8(c) and accumulated error of these results (integrated over time) in Figure 8(d). Here, there is still a clear convergence towards the analytical solution, but the pressure oscillations in time are still

AC

present in all resolutions up to H/∆x = 160 where they are greatly reduced.

This gives a sharp slope between H/∆x = 80 and H/∆x = 160, which shows that the COC is not as linear as might be concluded from Figure 8(b).

23

ACCEPTED MANUSCRIPT

CR IP T

0.55

0.45 2

3

4

AN US

0.5

5

6

7

8

9

10

-0.4

-0.3

-0.2

-0.1

0

(a)

0.6

M

0.5 0.4

0.2 0.1

PT

0

ED

0.3

CE

-0.1 -0.8

-0.7

-0.6

-0.5

(b)

AC

Figure 6: Comparison of SPH simulations using different dissipation schemes with the analytical solution: (a) normalised pressure at P0 over time in the interval tU/L = [2, 10]; (b)

Normalised pressure over the plate at time tU/L = 10.

24

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

CE

Figure 7: SPH simulation with particles coloured by non-dimensional pressure. Simulation parameters correspond to the highest resolution simulation of H/∆x = 160 using the van

AC

Albada limiter and artificial viscosity with α = 0.01.

25

0.6

0.4

10-2

0.2

0 -0.6

-0.4

-0.2

(a)

0

10

0.6

M

0.55 0.5

ED

0.45

4

40

80

160

80

160

(b)

0.65

0.4 2

10-3 20

AN US

-0.8

CR IP T

ACCEPTED MANUSCRIPT

6

8

10

10-2

20

40

(d)

PT

(c)

-1

Figure 8: Convergence tests using a particle spacing of H/∆x = 10 · 2n for n = 1, 2, 3, 4. (a)

CE

Normalised pressure over the plate at time tU/L = 10; (b) accumulated error between the pressure over the plate at tU/L = 10 in the SPH simulations and the steady state analytical solution, along with the computed order of convergence (COC); (c) normalised pressure at P0 over time in the interval tU/L = [8, 10]; (d) accumulated error between the pressure at

AC

P0 over time in the SPH simulations and the steady state analytical solution, along with the COC.

26

ACCEPTED MANUSCRIPT

3.4. Dam break A two-dimensional dam break simulation is chosen as a validation benchmark, which is commonly studied in the SPH literature [32, 17, 22]. The prob-

CR IP T

lem consists of an open top two-dimensional tank of width Lw = 5.336H with a fluid initially occupying a domain of width L = 2H and height H = 600mm

at the lower corner of the tank. Upon the start of the simulations, the fluid column is allowed to collapse and then plunges towards the far wall where it impacts the wall and collapses into the fluid as an overturning wave.

For the simulations, the fluid is initialised with hydrostatic pressure and the

AN US

particle spacing is set to H/∆x = 100, resulting in 20 000 fluid particles. The fluid is considered inviscid with density ρ0 = 1 000kg/m3 . The speed of sound is p chosen as c0 = 10 |gH| with g = −9.81m/s2 the acceleration due to gravity in

the z-axis direction. No artificial viscosity is used, i.e. α = 0 in equation (12),

and the parameter in the δ-SPH scheme is chosen as δ = 0.1.

Here, we compare the FRS scheme using the van Albada limiter [36] against

M

the MDSPH scheme. For this test, results for the ADSPH and ARS schemes are not shown because equivalence between them has already been demonstrated

ED

for previous test cases. Further, the ADSPH scheme is known to produce results very similar to MDSPH [22]. Snapshots of the particles coloured according to pressure at various times are shown in figure 9. Energy conservation components

PT

(excluding thermal) are compared in figure 10. Both schemes are capable of obtaining a smooth pressure field with low spatial oscillations and show highly

CE

similar energy evolution properties. Due to its simplicity the MDSPH scheme works well for short duration, advection dominated phenomena, however this case goes to show the general applicability of the FRS scheme, whilst avoiding

AC

the need of tuning the δ parameter.

27

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

Figure 9: Pressure field for the dam break simulation at times (top to bottom) 0s, 0.75s, 1.5s,

2s and 2.25s: (left) the Molteni δ-SPH scheme (MDSPH); (right) the full Riemann solver scheme using the van Albada limiter (FRS).

28

CR IP T

ACCEPTED MANUSCRIPT

0 0

-0.2 -0.4

-0.5

-0.6

-1 0

2

4

(a)

5

2

4

6

10-3

4

M

0.6

ED

0.4

2

4

3 2 1 0 6

PT

0 0

-1 0

6

(b)

0.8

0.2

AN US

-0.8

(c)

-1 0

2

4

(d)

CE

Figure 10: Energy evolution over time of SPH simulations comparing the δ-SPH Molteni scheme (MDSPH) and the full Riemann solver scheme using the van Albada limiter (FRS): (a) normalised total energy; (b) normalised potential energy; (c) normalised kinetic energy;

AC

(d) normalised internal energy.

29

6

ACCEPTED MANUSCRIPT

4. Conclusions We have presented a method for stabilising a weakly-compressible SPH formulation in the form of a Riemann solver that adds dissipation to the continuity

CR IP T

equation. The method is consistent at the free surface, preserves volume and reproduces hydrostatic pressures with good accuracy. It could be interpreted as a variant of the method proposed by Antuono et al. [16] without the tuning

parameter, δ, and with the amount of numerical dissipation depending on the choice of TVD flux limiter in the Riemann solver. We implemented the van Albada, van Leer, minmod and superbee limiters. Their numerical performance

AN US

was assessed by monitoring the evolution of the energy in well-established benchmark test problems for SPH simulations. For the hydrostatic simulation, the van Albada limiter produces the least dissipative results whilst preserving volume and consistency near the free surface. It is important to note that this conclusion is limited to the Wendland kernel and the other specific features chosen.

M

Further analysis of the limiters in relation to different kernels is left for future work. For the dynamic cases (jet impact and dam break flows) the resulting SPH scheme performs on a par with the schemes proposed by Antuono et al. [16] and

ED

Molteni and Colagrossi [15] in dissipating spurious pressure oscillations. Due to the absence of a tuning parameter, the proposed scheme could help reduce the

PT

complications associated with variable resolution and multiphase formulations where different tuning parameters may be required. Exploring this potential benefit is left for future work. Other ongoing work includes the thorough verifi-

CE

cation and validation in complex 3D cases, an assessment of the numerical costs and optimisation of the algorithms for GPU implementation.

AC

Acknowledgments This work was partially supported by the Italian Ministry of Education, Uni-

versities and Research under the Scientific Independence of Young Researchers project, grant number RBSI14R1GP, CUP code D92I15000190001. This work

30

ACCEPTED MANUSCRIPT

was also partially supported by the PRISM platform (Underpinning Technologies for Finite Element Simulation) under EPSRC grant EP/L000407/1.

CR IP T

Appendix A Given the particle nature of the SPH method, the energy of the system can

be derived from the Hamiltonian of a set of particles [7]. The total energy offers a unique way of assessing the quality of a SPH scheme by measuring its

energy dissipation. For instance, in a still water case an ideal scheme should perfectly conserve energy. Even though this is rarely the case due to various

AN US

choices of formulation, boundary conditions and numerical dissipation schemes,

monitoring the temporal evolution of the energy during the simulation provides a way of reliably quantifying potential negative effects of artificial dissipation and analyzing the performance of the various schemes.

The Hamiltonian H, a measure of the energy of a system, is given as the

M

sum of kinetic and potential energy

H = Ekin + Epot

(46)

ED

which is conserved for an isolated system in the absence of physical dissipation.

PT

The kinetic energy for a set of Np discrete particles is given as Ekin =

Np X

a=1

ma

va2 2

(47)

CE

where va2 is the dot product of the velocity, i.e. va · va . The potential energy can be further decomposed into its internal Eint and external Eext , components, where the internal energy depends on the interactions between the discrete

AC

particles, and the external energy depends on the interactions with an external system such as gravitational pull. The external potential energy due to gravity is

g Eext =

Np X

a=1

31

ma g · x

(48)

ACCEPTED MANUSCRIPT

The internal potential energy can be derived from the first law of thermodynamics δEint = δQ + δW

(49)

CR IP T

where δEint is the change in the internal energy of the system, δQ is the heat transfer into the system, and δW is the work done on the system. Energy

conservation equations and thermal effects are not included in the framework, so for the purposes of this work, δQ is dropped. The work done on the system

is given by the change in the volume of the system, referred to as the pressurevolume work, and is given by δW = −P dV , where P is the pressure of the

AN US

system and dV the change in volume, so that

δEint = −P dV

(50)

From the continuity equation and taking the volume as V = m/ρ, the change in the internal energy can be obtained as

M

δEint =

mP dρ ρ2

(51)

where the pressure is calculated from the equation of state (7). The internal

PT

ED

energy is then calculated by integrating δEint leading to    Z ρf Z ρf mp 1 ρ Eint (ρf ) = dρ + E (ρ ) = m p − 1 dρ + Eint (ρi ) int i 0 2 ρ2 ρ0 ρi ρi ρ (52) Evaluating the integral from the initial density ρi = ρ0 and the final density

CE

ρf = ρt , where t is the current time of the simulation, the internal energy reads     m p0 ρ ρ0 Eint (ρ) = ln + − 1 + Eint (ρ0 ) (53) ρ0 ρ0 ρ

AC

Finally, for a system of Np discrete particles the internal energy can be evaluated as Eint (ρ) =

    Np X ma p0 ρa ρ0 ln + − 1 + Eint (ρ0 ) ρ0 ρ0 ρa a=1

(54)

The total energy of the system Etot is then the sum of the kinetic energy Ekin g (47), the external (potential) energy Eext (48) and the internal energy Eint (54).

32

ACCEPTED MANUSCRIPT

For the purpose of post-processing, only the conservation of energy, and not its absolute value, are of interest, giving the final expression for the total energy as

where t0 is the initial time of the simulation.

(55)

CR IP T

g g Etot (t) = Ekin (t) + Eext (t) + Eint (t) − (Ekin (t0 ) + Eext (t0 ) + Eint (t0 ))

Other methods for calculating the total energy of the system can be found

in the SPH literature, for example [40]. However, they require the addition of the energy conservation equation to the scheme itself and cannot be performed

AN US

solely at the post-processing stage.

References

[1] R. A. Gingold and J. J. Monaghan, “Smoothed particle hydrodynamic: theory and application to non-spherical stars,” Monthly Notices of the Royal Astronomical Society, vol. 181, pp. 375–389, 1977.

M

[2] L. B. Lucy, “A numerical approach to testing the fission hypothesis,” The Astronomical Journal, vol. 82, no. 12, pp. 1013–1924, 1977.

ED

[3] J. J. Monaghan, “Smoothed particle hydrodynamics,” Annual Review Astronomics and Astrophysics, vol. 30, pp. 543–574, 1992.

PT

[4] ——, “Simulating free surface flows with SPH,” Journal of Computational Physics, vol. 110, pp. 399–406, 1994.

CE

[5] ——, “Smoothed particle hydrodynamics,” Reports on Progress in Physics, vol. 68, pp. 1703–1759, 2005.

AC

[6] ——, “Smoothed particle hydrodynamics and its diverse applications,” Annual Review of Fluid Mechanics, vol. 44, pp. 323–346, 2012.

[7] D. Violeau, Fluid mechanics and the SPH method: Theory and applications. Oxford University Press, 2012.

33

ACCEPTED MANUSCRIPT

[8] J. P. Vila, “On particle weighted methods and smoothed particle hydrodynamics,” Mathematical Models and Methods in Applied Sciences, vol. 9, no. 2, pp. 161–209, 1999.

CR IP T

[9] R. Xu, P. Stansby, and D. Laurence, “Accuracy and stability in incompress-

ible SPH (ISPH) based on the projection method and a new approach,” Journal of Computational Physics, vol. 228, no. 18, pp. 6703–6725, 2009.

[10] J. Nu˜ nez-Ramirez, J.-C. Marongiu, M. Brun, and A. Combescure, “A partitioned approach for the coupling of SPH and FE methods for transient

AN US

nonlinear FSI problems with incompatible time-steps,” International Journal for Numerical Methods in Engineering, vol. 109, no. 10, pp. 1391–1417, 2017.

[11] G. Fourtakas and B. D. Rogers, “Modelling multi-phase liquid-sediment scour and resuspension induced by rapid flows using Smoothed Particle Hy-

M

drodynamics (SPH) accelerated with a Graphics Processing Unit (GPU),” Advances in Water Resources, vol. 92, pp. 186–199, 2016.

ED

[12] M. Shadloo, G. Oger, and D. Le Touz´e, “Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and challenges,” Computers & Fluids, vol. 136, pp. 11–34,

PT

2016.

[13] C. Altomare, J. M. Dom´ınguez, A. J. C. Crespo, J. Gonz´ alez-Cao,

CE

T. Suzuki, M. G´ omez-Gesteira, and P. Troch, “Long-crested wave generation and absorption for SPH-based DualSPHysics model,” Coastal En-

AC

gineering, vol. 127, pp. 37–54, 2017.

[14] D. Le Touz´e, A. Colagrossi, G. Colicchio, and M. Greco, “A critical investigation of smoothed particle hydrodynamics applied to problems with freesurfaces,” International Journal for Numerical Methods in Fluids, vol. 73, no. 7, pp. 660–691, 2013.

34

ACCEPTED MANUSCRIPT

[15] D. Molteni and A. Colagrossi, “A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH,” Computer Physics Communications, vol. 180, no. 6, pp. 861–872, 2009.

CR IP T

[16] M. Antuono, A. Colagrossi, S. Marrone, and D. Molteni, “Free-surface flows solved by means of SPH schemes with numerical diffusive terms,” Computer Physics Communications, vol. 181, no. 3, pp. 532–549, 2010.

[17] S. Marrone, M. Antuono, A. Colagrossi, G. Colicchio, D. Le Touz´e, and

G. Graziani, “δ-SPH model for simulating violent impact flows,” Computer

1542, 2011.

AN US

Methods in Applied Mechanics and Engineering, vol. 200, no. 13, pp. 1526–

[18] R. Vacondio, B. D. Rogers, P. K. Stansby, P. Mignosa, and J. Feldman, “Variable resolution for SPH: A dynamic particle coalescing and splitting scheme,” Computer Methods in Applied Mechanics and Engineering, vol.

M

256, no. 0, pp. 132–148, 2013.

[19] P. N. Sun, A. Colagrossi, S. Marrone, M. Antuono, and A. M. Zhang,

ED

“Multi-resolution δ + –SPH with tensile instability control: Towards high Reynolds number flows,” Computer Physics Communications, vol. 224, pp. 63–80, 2018.

PT

[20] P. N. Sun, A. Colagrossi, S. Marrone, and A. M. Zhang, “The δ + –SPH model: Simple procedures for a further improvement of the SPH scheme,”

CE

Computer Methods in Applied Mechanics and Engineering, vol. 315, pp. 25–49, 2017.

AC

[21] A. Bermudez and M. V´ azquez-Cend´ on, “Upwind methods for hyperbolic conservation laws with source terms,” Computer & Fluids, vol. 23, pp. 1049–1071, 1994.

[22] M. Antuono, A. Colagrossi, and S. Marrone, “Numerical diffusive terms in weakly-compressible SPH schemes,” Computer Physics Communications, vol. 183, no. 12, pp. 2570–2580, 2012. 35

ACCEPTED MANUSCRIPT

[23] S. Inutsuka, “Godunov-type SPH,” Journal of the Italian Astronomical Society, vol. 65, pp. 1027–1031, 1994. [24] K. Puri and P. Ramachandran, “Approximate Riemann solvers for the

CR IP T

Godunov SPH (GSPH),” Journal of Computational Physics, vol. 270, pp. 432–458, 2014.

[25] A. N. Parshikov and S. A. Medin, “Smoothed particle hydrodynamics using interparticle contact algorithms,” Journal of Computational Physics, vol. 180, no. 1, pp. 358–382, 2002.

AN US

[26] L. Chiron, G. Oger, M. de Leffe, and D. L. Touz´e, “Analysis and improve-

ments of Adaptive Particle Refinement (APR), through CPU time, accuracy and robustness considerations,” Journal of Computational Physics, vol. 354, pp. 552–575, 2018.

[27] D. Avesani, M. Dumbser, and A. Bellin, “A new class of moving-least-

M

squares WENO SPH schemes,” Journal of Computational Physics, vol. 270, pp. 278–299, 2014.

ED

[28] X. Nogueira, L. Ram´ırez, S. Clain, R. Loub`ere, L. Cueto-Felgueroso, and I. Colominas, “High-accurate SPH method with multidimensional optimal order detection limiting,” Computer Methods in Applied Mechanics and

PT

Engineering, vol. 310, pp. 134–155, 2016. [29] A. Ferrari, M. Dumbser, E. F. Toro, and A. Armanini, “A new 3D parallel

CE

SPH scheme for free-surface flows,” Computers & Fluids, vol. 38, no. 6, pp. 1203–1217, 2009.

AC

[30] J. L. Cercos-Pita, R. A. Dalrymple, and A. Herault, “Diffusive terms for the conservation of mass equation in SPH,” Applied Mathematical Modelling, vol. 40, no. 19, pp. 8722 – 8736, 2016.

[31] H. Wendland, “Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree,” Advances in Computational Mathematics, vol. 4, no. 1, pp. 389–396, 1995. 36

ACCEPTED MANUSCRIPT

[32] S. Adami, X. Y. Hu, and N. A. Adams, “A generalized wall boundary condition for smoothed particle hydrodynamics,” Journal of Computational Physics, vol. 231, no. 21, pp. 7057–7075, 2012.

CR IP T

[33] M. G´ omez-Gesteira, B. D. Rogers, A. J. C. Crespo, R. A. Dalrymple,

M. Narayanaswamy, and J. M. Dom´ınguez, “SPHysics - development of a free-surface fluid solver. Part 1: Theory and formulations,” Computers & Geosciences, vol. 48, pp. 289–299, 2012, SPHysics is available at: www.sphysics.org.

AN US

[34] P. W. Randles and L. D. Libersky, “Smoothed particle hydrodynamics:

some recent improvements and applications,” Computer Methods in Applied Mechanics and Engineering, vol. 139, no. 1-4, pp. 375–408, 1996. [35] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems. bridge University Press, 2002.

Cam-

M

[36] G. D. van Albada, B. van Leer, and W. W. Roberts Jr., “A comparative study of computational methods in cosmic gas dynamics,” in Upwind and Springer, 1997, pp. 95–103.

ED

High-Resolution Schemes.

[37] P. L. Roe, “Characteristic-based schemes for the Euler equations,” Annual

PT

Review of Fluid Mechanics, vol. 18, no. 1, pp. 337–365, 1986. [38] B. van Leer, “Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme,” Jour-

CE

nal of Computational Physics, vol. 14, no. 4, pp. 361–370, 1974.

AC

[39] G. I. Taylor, “Oblique impact of a jet on a plane surface,” Philosophical Transactions for the Royal Society of London. Series A, Mathematical and Physical Sciences, pp. 96–100, 1966.

[40] A. H´erault, G. Bilotta, and R. A. Dalrymple, “Achieving the best accuracy in a SPH implementation,” in Proceedings of the 9th International SPHERIC Workshop, Paris, France, June 3–5 2014, pp. 134–139.

37