Reformulation of Smoothed Particle Hydrodynamics
with Riemann Solver
1

Shu-ichiro Inutsuka

Department of Physics, Kyoto University, Kyoto, 606-8502, Japan

E-mail: inutsuka@tap.scphys.kyoto-u.ac.jp

Smoothed Particle Hydrodynamics is reformulated in terms of the convolution of the original hydrodynamics equations, and the new evolution
equations for the particles are derived. The same evolution equation of
motion is also derived using a new action principle. The force acting oneach particle is determined by solving the Riemann problem. The use of
the Riemann Solver strengthens the method, making it accurate for the
study of phenomena with strong shocks. The prescription for the variablesmoothing length is shown. These techniques are implemented in strict
conservation form. The results of a few test problems are also shown.

1. INTRODUCTION

Smoothed particle hydrodynamics (SPH, [14, 7]) is a fully Lagrangian particle method to describe
the hydrodynamical phenomena. The Lagrangian particle methods are especially suited to hydrodynamical problems that have large empty regions and moving boundaries. Those problems naturally
arise in engineering science as well as geophysics and astrophysics. A variety of astrophysical problems
have been studied by SPH because of its simplicity in programming the two- and three-dimensional
codes and its versatility of incorporating the various physical e ects such as self-gravity, radiative cooling, and chemical reactions. A broad discussion of the method can be found in a review by Monaghan[18]. However, the inconsistency of the method is emphasized by Dilts [3, 4] who modi ed the methodby means of the moving-least-square basis functions to obtain accurate derivatives regardless of thepositions of the SPH particles. Another concern of the method is its poor description of the strong
shocks. In the two- or three-dimensional calculation of colliding gases, particles often penetrate into
the opposite side. This unphysical e ect can be partially eliminated by the so-called XSPH prescription [17] that does not introduce the (required) additional dissipation but results in the additional
dispersion of the waves. Therefore it is very desirable to construct a method that is simple and able
to describe the strong shock phenomena accurately.
In this paper, a new method for handling shocks in particle hydrodynamics is constructed. The
force acting on each particle is determined by solving the Riemann problem (RP). This use of the socalled Riemann Solver is introduced as an simple analogy of the grid-based method [6]. The previous
attempts to introduce Riemann Solver into the particle method failed to give the method an exact
1 to appear in Journal of Computational Physics, Vol. 179, pp.238-267, 2002

1

2

SHU-ICHIRO INUTSUKA

conservation form [9, 10]. This paper describes how to include the exact Riemann Solver into thestrictly conservative particle method.
The alternative approach to including small but sucient dissipation into the numerical solution isto nd a good limiter to switch a dissipative method and a less-dissipative method [5, 21]. In principle,
however, the switch is always an option to any numerical scheme including the present one, and its
discussion is beyond the scope of this paper.
Section 2 provides a description of the method where we derive the evolution equations for theparticles in terms of the convolution of the original hydrodynamic equations with the so-called kernel
function. The same evolution equations are derived from an action principle which is di erent fromthe previous ones [19, 20], which sheds light on the \hidden" approximation in the expression for
the velocity eld in SPH formalism. The detailed explanation for the implementation is described in
Section 3. The usage of the Riemann Solver is analogous to the grid-based second-order Godunov
scheme [24]. A variable smoothing length is also considered. Numerical examples involving strong
shocks are presented in Section 4. Section 5 is for summary.

2. THE METHOD

We consider the following set of equations for non-radiating inviscid uid:d = ;r  v ;dt(1)

dv = ; 1 rP;
dt


(2)

du = ; P r  v;
dt


(3)

P = ( ; 1)u;(4)

where u denotes the speci c internal energy, denotes the ratio of speci c heats, and other symbols
have their usual meanings.
In the standard SPH method, each particle has its own mass, and density at each particle's location
is simply assigned by

i =

X
j

mj W (xi ; xj ; h);(5)where subscripts denote particle labels, mj is the mass of the j -th particle, and W (x; h) is a sphericallysymmetric kernel function which is normalized to be unity if integrated in space:

W (x; h) = W (;x; h);Z

W (x; h)dx = 1;

(6)
(7)

and h is the parameter of the kernel function. For later convenience we also de ne the e ective width
he of the kernel function W (x) byh2eZ

 2 x2 W (x; h)dx:(8)

SPH WITH RIEMANN SOLVER

3

Although there are many possible forms for W (x; h), we use the following Gaussian kernel throughout
in this paper:dW (x; h) = hp1  e;x2 =h2 :(9)In the above equation, d is the number of dimensions, and h is the so-called smoothing length. Notethat he = h for this Gaussian kernel. We assume h is constant in space in Sections 2.1-3.2. The
spatially variable smoothing length is considered in Section 3.3 and the subsequent Sections.
The standard procedure to derive sets of evolution equations for SPH has been previously described[18]. In the following sections, we show how new evolution equations are derived from the originalhydrodynamics equations. We then show that our new equation of motion can be derived from theaction principle for an appropriately de ned Lagrangian function.

2.1. Key Concept

In general any method of computational uid dynamics inevitably brings errors into the solutions.
Therefore the essential feature of the method can be characterized by how to introduce the errors into
the solutions. For example, in many kinds of grid-based methods, the truncation of the Taylor-series
expansion is the origin of the errors in the numerical solutions. In the SPH formalism we consider the
convolution of the physical function f (x) with the kernel function,Z

hf i(x)  f (x0 )W (x0 ; x; h)dx0 :(10)

We use the symbol hi to denote the convolution. Taylor-series expansion of f around x = xi in Eq.(10)shows that the di erence of hf i(xi ) and f (xi ) is the second-order in he :2

hf i(xi ) = f (xi ) + h4e r2 f (x) + O(h4e ):(11)

Thus if we use hf i as an approximate solution for f , the errors of the second-order in he are introducedby the convolution with the symmetric function W .
We also haveZ0@f
h i(x) = @f (x ) W (x0 ; x; h)dx0@x0@ W (x0 ; x; h)dx0= ; f (x0 ) @x0Z
@ f (x0 )W (x0 ; x; h)dx0 ;
= @x

@x

Z

(12)

where we used the integration by part and Eq.(6). That is, the kernel convolution of rf is just rhf i.
If we de ne the density distribution by the summation of the kernel functions at particle positions:(x) 
then we have the following identities:
1=
0=

X

mj W (x ; xj ; h);(13)

mj W (x ; x ; h);
j
(x)
j

(14)

j

X

X
j





mj r W (x;(xx) j ; h) :

(15)

4

SHU-ICHIRO INUTSUKA

These equations are the key equations in the present method. Using Eq.(14) the value of the kernel
convolution can be cast into the following expression:Zf (x0 )W (x0 ; xi ; h)dx0
Z X0=
mj f((xx0 )) W (x0 ; xi ; h)W (x0 ; xj ; h)dx0jX=fi;j ;

fi  hf i(xi ) =jwhere we de nefi;j Z0

mj f((xx0 )) W (x0 ; xi ; h)W (x0 ; xj ; h)dx0 :

(16)

(17)

Note that mi fi;j is symmetric with respect to i and j . In this way the value (fi ) of the physicalvariable f at the i-th particle position is expressed as the summation of the contributions (fi;j ) fromthe surrounding particles. This expression shows the spirit of the present method, although the actual
detailed evolution equations are presented in the following sections.
In contrast the standard SPH adopts the following form for fi;j :fi;j  mj f((xxj )) W (xi ; xj ; h):j(18)This corresponds to the crude approximation, W (xj ; x0 ; h)  (x0 ; xj ) in Eq.(17), where (x) is
the Dirac's  funcion. The standard SPH also implicitly assumes the following approximate equation:1

X
j

mj W (x ; x ; h);(xj )

j(19)

instead of Eq.(14), although this is very poor approximation.
In the MLSPH scheme by Dilts, an equation analogous to Eq.(14) is provided by the movingleast-square basis functions, although he needed somewhat heavy operations to construct his basis
functions. Note that our method does not restrict the functional form of the kernel function to satisfy
Eqs.(13),(14).

2.2. Equation of Motion

To derive the appropriate evolution equation for particles we take the convolution of the equationof motion (EoM), Eq.(2),Z
dv(x) W (x ; x0 ; h)dx = ; Z 1 rP (x) W (x ; x0 ; h)dx:
(20)
dt
(x)
The right-hand side of this equation can be transformed into the following expression:


Z 
P
(
x
)
P
(
x
)
0
0
;
r (x) W (x ; x ; h) + 2 (x) [r(x)] W (x ; x ; h) dx
9
8
2
3
Z <
=
X
mj W (x ; xj ; h)5 W (x ; x0 ; h); dx= : P((xx)) r W (x ; x0 ; h) ; P2((xx)) 4r
j
=

X
j

mj

Z

P
0
0
2 frW (x ; x ; h)W (x ; xj ; h) ; W (x ; x ; h)rW (x ; xj ; h)g dx(21)

SPH WITH RIEMANN SOLVER

5

where we integrated by part and used the identity Eq.(14).
Thus if we adopt the following equation for the evolution of the particle positions,x i 

Z

dv (x) W (x ; x ; h)dx = ; Z 1 rP (x) W (x ; x ; h)dx;
i
i
dt
(x)

(22)

we obtain the evolution equation that is consistent and spatially second-order accurate to the originalhydrodynamical equation of motion:

mi x i = ;mi

X
j

mj

Z





@
P (x) @
2 (x) @ xi ; @ xj W (x ; xi ; h)W (x ; xj ; h)dx;(23)

where the overdot means a time-derivative. The antisymmetric appearance of i and j on the righthand side guarantees the conservation of the linear and angular momentum of the particle system.
From this equation we can deduce some properties of the EoM of the particle system. For example,
if the pressure distribution is constant in space, the acceleration vanishes exactly in Eq. (23), whichis evident in Eq. (22). The evolution equation of MLSPH also has a similar property although thestandard SPH does not [3].
In this way the main approximation introduced in the present method is simply expressed in Eq.
(22), while the standard SPH formalism does not have such simple explanation.

2.3. Energy Equation

We can follow the evolution of barotropic uid with P = P () only by the EoM derived in theprevious section. However an additional evolution equation for energy is required to describe theevolution of non-barotropic uid. The energy equation is also needed to handle shocks. The energy
equation guarantees the conservation of the total energy of the system in the absence of other sources
or losses of energy, such as radiative heating or cooling. On the other hand, the strict conservation
property in the numerical calculation guarantees the accurate description of strong shocks. This isbecause the structure of shock wave is determined by the Rankine-Hugoniot relation which is the
direct consequence of the physical conservation laws. Thus we must derive the energy equation thathas strict conservation property and a convenient form to include additional physical dissipation.
As in Eq. (20), we multiply both sides of Eq. (3) by W (x ; x0 ; h), and integrate in space (withrespect to x),Zdu(x) W (x ; x0 ; h)dx = ; Z P (x) [r  v] W (x ; x0 ; h)dx:
dt
(x)

(24)

The right-hand side of this equation can be transformed into the following expression:Z; P((xx)) [r  v] W (x ; x0 ; h)dx
ZZ
= ; (1x) [r  P v] W (x ; x0 ; h)dx + (1x) [v  rP ] W (x ; x0 ; h)dx

(25)

At this point, we again make an approximation according to the following equation (see also Eqs.(36),(38)):Z

1 [v  rP ] W (x ; x ; h)dx = Z 1 [x_  rP ] W (x ; x ; h)dx + O(h2 ):ii
(x)(x) i(26)

6

SHU-ICHIRO INUTSUKA

Then we haveZ

du(x) W (x ; x ; h)dxidtZZ ; (1x) [r  P v] W (x ; xi ; h)dx + (1x) [x_ i  rP ] W (x ; xi ; h)dxZ= P [v ; x_ i ]  r W (x ; xi ; h) dx
ZX=
mj P2 [v ; x_ i ]  [rW (x ; xi ; h)W (x ; xj ; h)

u_i 

j

;W (x ; xi ; h)rW (x ; xj ; h)] dx;(27)

where we have again used the identities (14), (15). This form of the energy equation will be further
modi ed to account for the physical dissipation in Section 3.2, and its conservation property will beshown in Section 3.4.

2.4. Action Principle

In this section we show another feature of the equation of motion in the present method. In the casewithout dissipation, the equation of motion for the particles in our method can be derived from anaction principle. Therefore the time evolution of the particle system can be formulated in a canonicaltransformation, that might further enable the possible improvement of the time integration.
For the moment, we consider barotropic uid in which pressure is a function of density alone,
P () = K . In this case, the equation of motion (Eq.[2]) for the uid can be derived by minimizing theaction S which is the time integral of the Lagrangian function L expressed in Lagrangian coordinates:Z

S = Ldt;
Z 

Z

(28)


L = Ldx = 21 v2 ; u dx:(29)
If we use Eulerian coordinates, the appropriate Lagrangian density L must include constraint equations
along with Lagrange multipliers [13, 22], which introduce additional complication. In contrast, theLagrangian density shown above is simple and analogous to the classical particle mechanics form, owing
to the use of Lagrangian coordinates. In addition we can derive the Hamiltonian and the evolution
becomes the canonical transformation which is free from dissipation. Although this formalism is forthe continuous media, we can use this formalism to derive the evolution equation for the system ofparticles.
First we consider the density eld expressed in the following equation:
(x) =

Xj

mj W (x ; xj ; h):(30)This continuous distribution of density is de ned only by (the nite number of) the positions ofparticles, xi (i = 1; 2; 3; :::; Np). Time-derivative of the above equation gives@(x) = X m x_  @ W (x ; x ; h) = ; X m x_ rW (x ; x ; h) = ;r  F ;
j j
j j
j
m
j
@t
@ xj
j
j

(31)

7

SPH WITH RIEMANN SOLVER

where we de ned
F m(x) Xjmj x_ j W (x ; xj ; h):(32)Comparing Eq. (31) with the ordinary continuity equation, we realize that F m de nes the mass ux,
F m (x) = (x)v(x)(33)

From this equation we can deduce the de nition of velocity eld given by the following equation:Pm x_ W (x ; xj ; h)
F
(
x
)
m
v(x) = (x) = Pj j j:(34)
j mj W (x ; xj ; h)Note that at x = xi the above equation givesXv(x ) = x_ + 1
m [x_ ; x_ ]W (x ; x ; h):
i

i

i j

j

j

i

i

(35)

j

The di erence of v(xi ) and x_ i is second order in h. In a one-dimensional case, for example, this canbe shown by using an appropriate smooth function fx that satis es fx(xi ) = x_ i (i = 1; 2; :::; Np) asvx (xi ) ; x_ i


2 fx
X
1
@f
1
@
x
2
3
mj @x [xj ; xi ] + 2 @x2 [xj ; xi ] + O([xj ; xi ] ) W (xi ; xj ; h)
= 
i

j





2
1 @ 2 fx  + O(h4 ):x @= 1 h2 @f+@x @x 2 @x2
i(36)

Eqs. (30) and (34) can be used for v(x) and (x) in Eq. (29) to obtain the Lagrangian functionthat is de ned only by the positions of particles.

L=

8h
P
Z >
<
j j j
P
>
:
j j

i2

2

3 9

>
=
m x_ W (x ; xj ; h)
X
K
4
5
;
m
dx:
j W (x ; xj ; h)>
2 m W (x ; xj ; h)
;1 j
;

(37)

This is the exact Lagrangian for the \ uid" of which density and velocity are constrained by Eqs. (30)and (34).
Instead of adopting this complicated Lagrangian function, we now make an approximation for thevelocity eld by Taylor-series expansion of the velocity eld:v(x) = v(xi ) + [x ; xi ]  rv(xi ) + O(jx ; xi j2 )X
= x_ i + 1 mj [x_ j ; x_ i ]W (xi ; xj ; h) + [x ; xi ]  rv(xi ) + O(jx ; xi j2 ):i

From this we have

(38)

j

Z

v(x)W (x ; xi ; h)dx = x_ i + O(h2 ):With this equation we can simplify the kinetic term of Lagrangian function as follows:Z
1 v2 dx = 1 Z F  vdx
m
2
2

(39)

8

SHU-ICHIRO INUTSUKA

= 21

X
j

Z

mj x_ j  vW (x ; xj ; h)dx

X= 12 mj x_ 2j + O(h2 ):j(40)

Now we can de ne a simple Lagrangian function that is second-order accurate in h as

Lnew 

X
i





Z

mi 12 x_ 2i ; u(x)W (x ; xi ; h)dx :(41)
RWe can derive the evolution equation for the particles by minimizing the action S = Lnew dt . In
fact the Euler-Lagrange's equation,d @Lnew ; @Lnew = 0;
dt @ x_ i
@ xi

(42)

gives
x i = ;

X
j

mj

Z





P (x) @@
2 (x) @ xi ; @ xj W (x ; xi ; h)W (x ; xj ; h)dx:(43)

The manipulation to derive Eq. (43) from Eq. (42) is explained in Appendix. This equation is thesame as Eq. (23).
This is the exact evolution equation for the system of which Lagrangian function is de ned by
Eq. (41). The space-symmetry of Lnew guarantees the conservation of linear momentum and angularmomentum, which is also obvious in the anti-symmetric appearance of i and j in Eq. (43).
If we start with more simpli ed Lagrangian function,


LSPH = mi 21 x_ 2i ; u(i ) ;
i
X

(44)

we would obtain the standard SPH equation [19],"#x i = ; mj P2i + P2j @@x W (xi ; xj ; h):
i
i
j
j
X

(45)

This Lagrangian function LSPH for the standard SPH can be obtained if we make a crude approximation, W (x ; xi ; h)  (x ; xi ) in Eq.(41). Thus the di erence of the evolution equations, Eqs.
(43) and (45) is related to the di erence of the degrees of approximations to the original Lagrangianfunction of the uid. In this context, the standard SPH is also assuming the approximation for thevelocity eld Eq.(39) as in the present method.
Dilts [3] reported that the SPH approximation is derived by means of the Galerkin approximation
followed by the kernel approximation. In contrast to the Galerkin approximation, the above action
principle has deep physical consequences such as the conservation of the linear and angular momentum of the particle system. In addition Hamiltonian structure of the method possibly gives further
sophistication to the time-integration scheme (see, e.g, the symplectic integrator for the astronomical
self-gravitating system [25]), although it is beyond the scope of this paper.

3. IMPLEMENTATION

9

SPH WITH RIEMANN SOLVER

In this section we describe the numerical implementation of the method in detail. In Section 3.1
we describe how to evaluate the integrals in the basic evolution equations. Section 3.2 is for the
introduction of Riemann Solver into the method. Section 3.3 shows the prescription for the variablesmoothing length. In Section 3.4 we show that the present method remains a fully conservative scheme
even after the discretization in space and time.

3.1. Convolution

In this subsection we describe how to evaluate the spatial integrals in Eq.(23) and (27). First
we have to know the distribution of ;2 (x) to calculate the integrals. That is, we have to constructappropriate interpolation (or extrapolation) function for ;2 (x) around each pair (i and j ) of particles.
For convenience, we de ne the s-axis which is along the vector xi ; xj and has the origin at(xi ; xj )=2. We use s? to symbolically denote the components of the axes that are perpendicular
to the s-axis, and de ne the s-coordinate system. The unit vector in the s-direction is ei;j  (xi ;xj )=jxi ; xj j: We set si;j  si ; sj = jxi ; xj j, where si and sj denote the s-components of thepositions, xi and xj , respectively.
We de ne the speci c volume and its gradient as follows:

V (x) = (1x) ;rV (x) = ; 21(x) r(x) = ; 21(x)

(46)
X
j

mj rW (x ; xj ):(47)

We will make approximate function for V (x) in the s-coordinates. In the following subsections, twosets of equations with the di erent order of accuracy are shown.
3.1.1.

asLinear Interpolation

The most simple choice for the approximate function is the linear interpolation which is expressed

V (s) = (1s) = Ci;j s + Di;j ;

(48)

where
) ; V (xj ) ;Ci;j = V (xisi;jD = V (xi ) + V (xj ) ;
i;j

(49)

2

From this we obtained

V 2 (s) = 21(s) = [Ci;j s + Di;j ]2 ;With this interpolation, we can calculate the integral as follows:
Z

(50)

p

;2 (x)W (x ; xi ; h)W (x ; xj ; h)dx = Vi;j2 (h)W (xi ; xj ; 2h);(51)

2 + D2 :Vi;j2 (h)  14 h2 Ci;j
i;j

(52)

where we de ned

10

SHU-ICHIRO INUTSUKA

The above calculation corresponds to an approximation based on the linear expansion of 1=2 in theperpendicular direction to the vector xi ; xj , that is,

;2 (s) = V 2 (s) + s?  rV 2 (s) + O(js? j2 ):


(53)



Note that the integral of s?  rV 2 (s) in Eq. (51) vanishes identically owing to the symmetry of thekernel.
Next, we consider the following integral for arbitrary function f (x):Zf (x) W (x ; x ; h)W (x ; x ; h)dx:(54)ij2 (x)
When we have some interpolation for the function f (x) and ;2 (x) , we can de ne the weighted
 by the following equation:
average fi;jZ

f (x) W (x ; x ; h)W (x ; x ; h)dx = f  Z 1 W (x ; x ; h)W (x ; x ; h)dx:
i
j
i
j
i;j
2 (x)
2 (x)(55)

A straight-forward calculation of the above equation with the linear approximation for the function
f (x) ( s[fi ; fj ]=si;j + [fi + fj ]=2) and the speci c choice (Eq.[50]) of the interpolation of ;2 (x)gives = fi ; fj s + fi + fj ;
fi;j
(56)
si;j i;j
2
where
2

si;j = h C2Vi;j2Di;j :

(57)

i;j

 is the value of the (linearly approximated) function f at the position si;j .
That is, fi;j
The evolution equations now become the following:p
mj P  Vi;j2 (h) @@x W (xi ; xj ; 2h);(58)p
mj ([P v ] ; P  x_ i )Vi;j2 (h) @@x W (xi ; xj ; 2h);(59)x i = ;2

u_i = ;2

X
j

X
j

i

i

where the formal meanings of the P  and [P v] are the values at the position si;j of the linearlyinterpolated functions. In section 3.2, however, we will change the values of these quantities byconsidering the physical dissipation.
3.1.2.

Cubic Spline Interpolation

Another convenient method for the approximation is based on the cubic spline interpolation of ;1along the s-axis :(60)
V (s) = (1s) = Ai;j s3 + Bi;j s2 + Ci;j s + Di;j ;where
Vi ; Vj ) + (Vi0 + Vj0 ) ;Ai;j = ;2 ((si;j )3 (si;j )2

11

SPH WITH RIEMANN SOLVER(V 0 ; V 0 )Bi;j = 21 is j ;i;jCi;j = 32 (Vi ;s Vj ) ; 41 (Vi0 + Vj0 );i;jDi;j = 12 (Vi + Vj ) ; 18 (Vi0 ; Vj0 )si;j ;Vi = V (xi );Vj = V (xj );Vi0 = ei;j  rV (xi );Vj0 = ei;j  rV (xj ):

(61)

Then,

V 2 (s) = 21(s) = Ai;j s3 + Bi;j s2 + Ci;j s + Di;j 2 ;




(62)

In this case Eq. (51) becomesZp;2 (x)W (x ; xi ; h)W (x ; xj ; h)dx = Vi;j2 (h)W (xi ; xj ; 2h);(63)

15 h6 A2 + 3 h4 (2A C + B 2 ) + 1 h2 (2B D + C 2 ) + D2 :Vi;j2 (h)  64i;j i;ji;j i;j
i;j
i;j
i;j
i;j
16
4

(64)

where we de ned

Eq. (57) becomes
15 h6 Ai;j Bi;j + 3 h4 (Ai;j Di;j + Bi;j Ci;j ) + 1 h2 Ci;j Di;j
8
2
:
V2

si;j = 32

i;j

(65)

To avoid undershooting and overshooting of the interpolating function, we need to use linear interpolation (Eq.[48]) in the case Vi0 Vj0 < 0.

3.2. The Usage of the Riemann Solver

The description of shock waves requires (physical) dissipative process which is not considered in thefundamental equations (Eqs.[1]-[4]). The Godunov's scheme and its second-order sequel (MUSCL, [24])and third-order sequel (PPM, [2]) use the exact Riemann Solver to include the (possibly) minimum
and sucient amount of dissipation into the method. At present these methods remain state-of-theart grid-based methods for computational uid dynamics. In this paper we introduce the exact andspatially second-order Riemann Solver into the particle method.
The Godunov's scheme uses the result of Riemann problem (RP) at each cell interface in thecalculation of numerical ux (see, e.g., [24]). Likewise, we want to use the result of RP at the vicinityof the middle point of the i-th particle and the j -th particle. This is achieved by modifying the valuesof P  and [P v] in Eqs. (58),(59). The nite di erence expression for this is the following:

x_ i = ; X m P  Z 1
@
@
j i;j
2 (x) @ xi ; @ xj W (x ; xi ; h)W (x ; xj ; h)dx
t

j
= ;2

X
j

p
mj P  Vi;j2 (h) @@x W (xi ; xj ; 2h);
i

(66)

12

SHU-ICHIRO INUTSUKAui = ; X m P  (v  ; x_  ) Z 1@ ; @ W (x ; x ; h)W (x ; x ; h)dx
j i;j i;j
i
j
i
2
t
 (x) @ xi @ xj
j
= 2

X
j

p
mj (P  v ; P  x_ i )Vi;j2 (h) @@x W (xi ; xj ; 2h):

(67)

i

where  stands for the nite di erence of each variable, and x_ i is the time-centered velocity of thei-th particle, and Pi;j and vi;j are the results of RP between the i-th particle and the j -th particle.
How to de ne and calculate these variables are explained below.
a)

f(s)

b)
fR

fi
fL

f*i,j
fj

fj

Ci dt

Cj dt

sj

fi

si

s* 0

s

sj

si

s*

s

FIG. 1.
a) A schematic picture of the distribution of the function f (s) in the vicinity of the i-th particle andthe j -th. The value of the weighted average f  are determined by Eqs. (56), and its location s by Eq. (57) or (65).
b) The setup of Riemann problem. We de ne the piecewise linear distribution of physical variable f (s) in both the i-th
region and the j -th region. The initial values of each side of the one-dimensional Riemann problem are the averagevalues in each domain of dependence.
i;j

i;j

Figure (1) shows schematic picture of the distribution of the function f (s) in the vicinity of the i-th
particle and the j -th. According to the original (grid-based) Godunov's scheme, it is natural to de ne
the interface of the two regions for the RP at si;j . To make the spatially second-order method like
MUSCL, we de ne the piecewise linear distribution of physical variable f (s) in both the i-th region
and the j -th region. The gradient of f (s) in each region can be simply assigned by the gradient at
each particle's position. The initial values on each side of one-dimensional Riemann problem is theaverage values of each domain of dependence, which is expressed as follows:


 



 + Cs;i t ; si ;sR = i + @
@s i i;j
2

 

@P
PR = Pi + @s si;j + Cs;i 2t ; si ;
i

 

@
v
vR = vi + @s si;j + Cs;i 2t ; si ;
i

 

@
L = j + @s si;j ; Cs;j 2t ; sj ;
j
 



t
@P

PL = Pj + @s si;j ; Cs;j 2 ; sj ;
j

 

@v
vL = vj + @s si;j ; Cs;j 2t ; sj ;j(68)SPH WITH RIEMANN SOLVER

13

vi = vi  ei;j ;vj = vj  ei;j :

(69)

where

and Cs;i denotes the sound speed at xi . Thus we can solve the RP at the interface si;j . The detailed
explanation for the solution of RP can be found in [24, 2], and we will not repeat it here.
 and the projected velocity vi;j.
After solving the RP at si;j , we have the resulting pressure Pi;j
Three-dimensional (or two-dimensional) velocity vi;j is calculated as follows: ei;j + [vi;j ; vi;j ei;j ];vi;j = vi;j(70)

vi;j = vi [1=2 + ] + v j [1=2 ; ];(71)vi;j = vi [1=2 + ] + vj [1=2 ; ];s = si;j :
i;j

(72)

where

(73)

 and vi;j are used in the calculation of the right-hand side of Eqs. (66) and (67). Note that
These Pi;jEqs.(71), (72), (73) are not needed in the actual calculation, because the bracketed second term on
the right-hand side of Eq. (70) is perpendicular to ei;j and vanishes identically in Eq.(67).
In the higher-order grid-based methods (e.g., MUSCL and PPM) we need to impose monotonicity
constraint on the gradients of physical variable to obtain the stable description of the discontinuity.
This is also true in the present method. Our experience shows that the monotonicity constraint only
on the velocity eld suces to make the stable method. Thus in the actual calculation of our method,
we mimic the monotonicity constraints by setting@v  =  @v  = 0 if
@s i
@s j



@v    @v  < 0:@s i @s j(74)

The numerical scheme should be rst-order in space at the surface of the shock wave. The shock
surface tends to include a few particles in the actual calculation. Therefore, if the velocity di erenceof a certain particle pair corresponds to the sound speed divided by a small number, the pair isconsidered to be in the shock surface, and we need to use the rst-order Riemann Solver for this pair.
We implement this condition by setting








@f = @f = 0 if C
(75)
shock ei;j  (v j ; v i ) > min(Cs;i ; Cs;j );
@s i
@s j
where f = ; P; v, and Cshock is a numerical constant corresponding to the number of particles atthe shock surface. We adopt Cshock = 3 throughout in this paper.

3.3. Variable Smoothing Length

In the previous subsections, we assume that the smoothing length, h, is constant in space. In actual
calculation we need to change h to enlarge the dynamic range of the spatial resolution. Therefore wehave to extend the present method for the variable smoothing length.
The formal derivation of the evolution equation with variable smoothing length is the same as
Section 2.2, except that we should start with the de nition of density as

(x) =

Xj

mj W (x ; xj ; h[x]):

(76)

14

SHU-ICHIRO INUTSUKA

This de nition of density corresponds to the so-called \gather" formulation of SPH [8]. The resulting
evolution equations are the following:x i =

u_ i =Xj

mj P 

Zi;j1

2 (x) frW (x ; xi ; h[x])W (x ; xj ; h[x]);W (x ; xi ; h[x])rW (x ; xj ; h[x])g dx;(77)

mj Pi;j (vi;j ; vi ) 21(x) frW (x ; xi ; h[x])W (x ; xj ; h[x])j;W (x ; xi ; h[x])rW (x ; xj ; h[x])g dx:

(78)

Z

X

These equations are essentially the same as Eqs. (66) and (67), although the analytic integration
for these equation is not possible even with the polynomial approximation for ;2 (x). Therefore we
prefer to use simple approximation for the integral as in the following form:
x_ i = ; X m P  V 2 (h ) @ W (x ; x ; p2h )
j
i
j
i
i;j i
t
@ xi
j


p



+Vi;j2 (hj ) @@x W (xi ; xj ; 2hj ) ;

(79)

i

ui = ; X m [P  v ; P  x_  ] V 2 (h ) @ W (x ; x ; p2h )
j
i
j
i
i
i;j i
t
@ xi
j


p



+Vi;j2 (hj ) @@x W (xi ; xj ; 2hj ) ;i(80)

where, in spirit, we used hi for the half of the integration space which include xi , and used hj for theother half.
The above approximations assume that the h(x) should not vary much within the neighborhoodof each particle. One possible way to determine the smoothing length with this constraint is thefollowing:1=di;
hi =  m





i

(81)

where

i =

X
j

mj W (xi ; xj ; hi ); hi = hi Csmooth :(82)

 is more smooth than  itself if Csmooth > 1 . Numerical experiments shows that  ' 1 with
Csmooth ' 2 works ne in Section 4. The e ective number of neighbors around each particle depends

on the ratio of the smoothing length and the mean separation of particles at that particle. For
example, we can usually ignore the contribution from the j -th particle to the i-th particle if theirdistance jxi ; xj j is larger than 3hi , because exp(;32 )  1:234  10;4. Thus, in one-dimensional
calculations, the number of neighbors for calculating Eq.(82) is 6Csmooth excluding the i-th particle2
itself. The number of neighbors becomes about 282Csmoothin two-dimensional calculations, and33
about 113 Csmooth in three dimensional calculations.

15

SPH WITH RIEMANN SOLVER

3.4. Conservation Property

The nal discretized form for the equation of motion for the i-th particle becomes the following:x_ i = ;tXjpmj Pi;j Vi;j2 (hi ) @ W (xi ; xj ; 2hi )


@ xi

p

+Vi;j2 (hj ) @@x W (xi ; xj ; 2hj ):



(83)

i

This calculation guarantees the conservation of the total momentum of the system,Ximi x_ i = 0;(84)

because the terms in the square bracket in Eq. (83) are anti-symmetric in i and j .
The nal form for the equation of energy for the i-th particle becomes the following:ui = ;tXjp
mj P  (v  ; x_  ) V 2 (hi ) @ W (xi ; xj ; 2hi )


i;j

i;j

i

@ xi

i;j

p
+V 2 (h ) @ W (x ; x ; 2h )
j

i;j

@ xi

i



j

j

where we de ned the time-centered velocity of the i-th particle,
x_ i = x_ i + 1 x_ i :2This expression guarantees the conservation of the total energy of the system, because


X

=

i



mi 12 x_ 2i + ui

X
i
X







=
mi x_ i x_ i + 12 x_ i + ui
i
= ;

XX
i

j





p
mi mj Pi;j vi;j Vi;j2 (hi ) @ W (xi ; xj ; 2hi )


@ xi

p
+V 2 (h ) @ W (x ; x ; 2h )
i;j

= 0:

(86)



mi 12 [x_ i + x_ i ]2 + [ui + ui ] ; 21 x_ 2i ; ui


(85)

j

@ xi

i

j



j

(87)

Thus the present scheme conserves energy exactly. This is in contrast to the ordinary energy equationof the standard SPH which is only accurate in the rst-order in time (see, e.g., [1, 18]). We also notethat the conservation of energy does not require constant smoothing length, which isPobvious in Eq.
(87). This is because the expression of the total energy in the numerical calculation i mi ( 21 x2i + ui )does not explicitly depends on the choice of smoothing length. In other words, a sudden change ofsmoothing length at any timestep does not change the numerical value of the total energy.

3.5. Overall Procedure

In this section we summarize the actual procedures in the sequence of executions. The main loopof the time integration corresponds to Steps 1-4.

16

SHU-ICHIRO INUTSUKA

Step 0. Problem Setup

We rst setup the problem in the computer program. We appropriately place the particles torepresent the density distribution that corresponds to the initial condition of the problem to solve.
This may require some relaxation technique to nd the appropriate positions of the particles [18]. We
also determine i , ri , and hi to start the main loop of the time integration. Various constants andthe initial variables are calculated in this step. For example we determine the initial timestep t.

Step 1. Gradient Calculation

We calculate the gradient of physical variables P; v for the use in the Riemann Solver.

Step 2. Source Term Calculation

We calculate the RHS of Eqs.(83),(85) either by Eq.(52) or Eq.(64). The subroutine for the
Riemann Solver is called once for every pair of particles to calculate P  and v in Eqs.(83),(85).

Step 3. Time Evolution
We update xi , x_ i , and ui according to Eqs.(86),(83),(85) ,xi (t + t) = xi (t) + x_ i t;x_ i (t + t) = x_ i (t) + x_ i ;ui (t + t) = ui (t) + ui :(88)

Step 4. Density Updation

According to the updated positions of particles, we update density distribution. The smoothinglength of each particle is also updated. The timestep t for the next integration is also determined.
We turn to Step 1 for the next time integration.

4. NUMERICAL EXAMPLES

The present method was tested on a variety of 1D, 2D, and 3D problems, a few of which aredescribed below. Other sets of test calculations will be described elsewhere.
In determining the amount of the timestep t, we have to consider the Courant condition that isin spirit similar to the Courant condition for the grid-based Lagrangian methods,
t = CCFL mini(mi 1=di,)Cs;i :(89)

The numerical experiments show that we can safely use CCFL  0:5 in most of hydrodynamical
problems. Note that we don't need to consider the (e ective) di usion timescale that is related to thearti cial viscosity adopted by the standard SPH and other particle methods. Thus t in the present
method can be much larger than that in the other SPH methods.
The following examples are calculated with the Fortran program in which the single precision real
number is used. To accelerate the computation we use the data structure based on link lists [16].
Five cycles of iterations in the Riemann Solver is sucient for the following test problems.

4.1. Shock Tube

In order to compare the capabilities of the present method and standard SPH we test our method
against the shock tube problems, of which the exact solutions are available. The density variation in
the rst problem is not so large that we can use both the variable smoothing length and the constantsmoothing length. The initial parameters of the problem are the following:L = 1 ;PL = 1 ;vx;L = 0 ;R = 0:5;
PR = 0:2;
vx;R = 0;

(90)

17

SPH WITH RIEMANN SOLVER

3

u

2.5
2
1.5
1
.5
.8

vx

.6
.4
.2

ρ

P

0
-.2
1.2
1
.8
.6
.4
.2
0
1.2
1
.8
.6
.4
.2
0
-.4

-.2

0

.2

.4

x
FIG. 2.
Results of the shock tube problem in which the Mach number is 1.526. The right-hand side of the initial
discontinuity includes 40 equal-mass particles. The open circles plot the snapshot at t = 0:2 by the present method withthe cubic spline approximation for the convolution Eq.(64) and the variable smoothing length ( = 1, Csmooth = 2:0).
The solid lines correspond to the analytic solution.

18

SHU-ICHIRO INUTSUKA

where the subscript \L" denotes the variable on the left-hand side of the initial discontinuity, andR denotes the right-hand side. The ratio of speci c heats is = 7=5. The Mach number of theresultant shock wave is 1:526. The value of the post-shock pressure is P  = 0:5099. Figure 2 plots thesnapshot at t = 0:2 by the present method with the cubic spline approximation for the convolution
Eq.(64) and the variable smoothing length ( = 1, Csmooth = 2:0). The method with the linear
approximation for the convolution Eq.(52) gives very similar result. Figure 3 plots the corresponding
result of the standard SPH where we used the \standard" arti cial viscosity ( = 1; = 2) describedin the review paper by Monaghan [18]. In both cases, the number of equal-mass particles is 80 (40)on the left(right)-hand side of initial discontinuity (i.e., xL = 0:005, xR = 0:01 ). The solid lines
correspond to the analytic solution. In this problem, these two methods give similar results except
for the contact discontinuity where the standard SPH produced a \wiggle" in pressure and speci cinternal energy. This is due to the inconsistency of EoM of the standard SPH (45). The present
method produced the smooth distribution of internal energy at the contact discontinuity, and hence
provide the almost constant pressure distribution at the density discontinuity.
The total energy of the system (;0:4 < x < 0:4) is de ned as
Z 0:4

1 u dx = 1:2:
2
;0:4

(91)

The initial numerical value of the total energy of the particle system was the following:Ximi ui = 1:20000017(92)

where the last few digits have no signi cant meaning in this Fortran single precision calculation. The
nal (t = 0:2) value of the total energy of the particle system was found to beXimi ( 21 x_ 2i + ui ) = 1:20003033(93)

The relative error remains suciently small (E=E < 10;4) when we change t in the other runs.
Thus the error of the total energy is only due to the round-o error of the single precision calculation,and not due to the truncation errors in the numerical modeling of the evolution. This result guarantees
the strict conservation property of the present scheme (see Section 3.4).
The present method spent 0.17 seconds for the total 52 timesteps (t  0:004) to compute with
the Hewlett-Packard workstation C240 (PA-RISC 8200/236MHz), which corresponds to 3:27  10;3
seconds per step and 2:72  10;5 seconds per particle. In contrast the standard SPH method used 956
timesteps (t  0:0002) for the stable evolution in this problem, because the arti cial viscosity limits
t. Indeed numerical experiments shows larger t causes unphysical oscillations in the solution. As
a result the standard SPH method spent 1.35 seconds to nish the calculation. It corresponds to
1:4  10;3 seconds per step and 1:2  10;5 seconds per particle. This means that the present method
is about two times slower than the standard SPH for the required operations per particle but thepresent method tends to spend much smaller (total) computation time than the standard SPH does.
Figure 4 plot the result of the present method where we used the constant smoothing length h = 0:01.
The rarefaction wave on the left hand side is described less-accurately than in Fig. 2 as expected,
because the smoothing length in high density region in Fig. 2 is smaller than the constant smoothinglength in Fig. 4. The method with the constant smoothing length has no advantage over the methodwith the variable smoothing length even in this kind of shock tube problem where the density variation
is not so large.
Figure 5 plot the result of the present method with the rst-order Riemann Solver where we set@=@s = @P=@s = @v=@s = 0 in the Riemann Solver. The sharp pro les at the shock front and the19

SPH WITH RIEMANN SOLVER

head and tail of the rarefaction wave are smeared out as in the grid-based rst-order Godunov method.
In this calculation we don't need to calculate the gradients of P and vx . However the di erence of thecomputation times with the rst- and second-order method is negligible. There is no reason to use
this rst-order method.

3

u

2.5
2
1.5
1
.5
.8

vx

.6
.4
.2

ρ

P

0
-.2
1.2
1
.8
.6
.4
.2
0
1.2
1
.8
.6
.4
.2
0
-.4

-.2

0

.2

.4

x
FIG. 3.
Results of the shock tube problem in which the Mach number is 1.526. The open circles plot thesnapshot at t = 0:2 by the standard SPH. The solid lines correspond to the analytic solution.

20

SHU-ICHIRO INUTSUKA

3

u

2.5
2
1.5
1
.5
.8

vx

.6
.4
.2

ρ

P

0
-.2
1.2
1
.8
.6
.4
.2
0
1.2
1
.8
.6
.4
.2
0
-.4

-.2

0

.2

.4

x

FIG. 4.

Same as Fig. 2 except that we used the constant smoothing length h = 0:01.

21

SPH WITH RIEMANN SOLVER

3

u

2.5
2
1.5
1
.5
.8

vx

.6
.4
.2

ρ

P

0
-.2
1.2
1
.8
.6
.4
.2
0
1.2
1
.8
.6
.4
.2
0
-.4

-.2

0

.2

.4

x

FIG. 5.

Same as Fig. 2 except that we used the rst-order Riemann Solver.

22

SHU-ICHIRO INUTSUKA

In Fig. 6, the results of the standard Sod's shock tube [23] are plotted. The initial parameters of
the problem are the following:L = 1 ; R = 0:125;
PL = 1 ; PR = 0:1;
vx;L = 0 ; vx;R = 0:(94)

The ratio of speci c heats is = 7=5. The solid lines correspond to the analytic solution. The numberof the equal-mass particles in the x-direction is 40 on the right-hand side of the initial discontinuity.
We used the variable smoothing length ( = 1, Csmooth = 2:0).

4.2. Extreme Blast Wave

To explore the capability of the present method, we test our method against an extremely strong
shock tube problem. The initial parameters of the problem are the following:L = 1 ; R = 1;
PL = 3000 ; PR = 10;7;
vx;L = 0 ; vx;R = 0;(95)

where the initial pressure of the gas on the left-hand side is 3  1010 times that of the right-hand side.
The ratio of speci c heats is = 5=3. The Mach number is as large as 105. We used 100 particles
on each side of the initial discontinuity. Figure 7 plots the result of the present method where the
cubic spline approximation for the convolution Eq.(64) and the variable smoothing length is used( = 1, Csmooth = 2:0). Even in this severe problem, the present method gave stable and accurate
results and is free from penetration problems. The pressure distribution shows a small wiggle at thecontact surface. This is due to the approximation for the convolution. Figure 8 plots the result of
the present method with the linear interpolation Eq.(52). The amplitude of the pressure wiggling at
the contact surface is slightly larger than that of the cubic spline result. In order to eliminate thisunphysical wiggling, we need to develop more accurate approximation for the numerical convolution
than Eq.(64). Numerical experiments shows that the standard SPH cannot produce acceptable resultat least with  = 1,  1, and  2. The possible ne-tuning of the parameters in the standard
SPH may enable the calculation in this case. Note that the present method can describe this extreme
phenomena with the same parameters ( = 1, Csmooth = 2:0) used in the other test problems, without
any tuning of the method.
The accurate calculation with the single precision real numbers in the Fortran program is due to theconservative formulation with the internal speci c energy u instead of the total energy e (= v2 =2 + u)
or E (= v2 =2 + u) which are usually used in the grid-based conservative numerical schemes (seeSection 3.4). If we adopt e or E as the main variable for the integration, we have to do the subtraction
to obtain u (= e ; v2 =2), which brings huge error into the numerical value of u in the extremely
supersonic motion (i.e., when u  e  v2 =2). In this sense our method has potential advantage overthe grid-based higher-order Godunov methods, at least, in describing extremely supersonic ows.

4.3. Wind-Tunnel

To test the present method against multi-dimensional problem, we adopt the following wind-tunnel
problem [12], that can be in principle directly compared with the laboratory experiment. The geometry
is a two-dimensional channel with a 15 wedge on the lower wall. A 15 expansion corner is also
included. The in ow Mach number is two. This kind of problem is poorly suited for the particle
method because a rigid-wall boundary condition must be set up at the wall of the tunnel. We present
this test problem only to demonstrate the capability of our scheme.

23

SPH WITH RIEMANN SOLVER

4

u

3.5
3
2.5
2
1.5

vx

1
.5

ρ

P

0
1.2
1
.8
.6
.4
.2
0
1.2
1
.8
.6
.4
.2
0
-.4

-.2

0

.2

.4

x
FIG. 6.
Results of the Sod's shock tube problem by the present method with the cubic spline approximation
for the convolution Eq.(64) and the variable smoothing length ( = 1, Csmooth = 2:0). The solid lines correspond to
the analytic solution.

24

SHU-ICHIRO INUTSUKA

Mach Number = 105
5000

u

4000
3000
2000
1000
0
40

vx

30
20
10
0

P

3000
2000
1000
0
5

ρ

4
3
2
1
0
-.4

-.2

0

.2

.4

x
FIG. 7.
Results of 1D calculation for the extremely strong shock tube problem of which Mach number is about
105 . We used 100 particles on each side of the initial discontinuity. The open circles plot the result of the present
method where the cubic spline approximation for the convolution Eq.(64) and the variable smoothing length are used( = 1, Csmooth = 2:0). A small wiggle in the pressure distribution at the contact surface is due to the approximation
for the convolution.

25

SPH WITH RIEMANN SOLVER

Mach Number = 105
5000

u

4000
3000
2000
1000
0
40

vx

30
20
10
0

P

3000
2000
1000
0
5

ρ

4
3
2
1
0
-.4

-.2

0

.2

.4

x
FIG. 8.
Same as Fig.7 except that we used the linear interpolation in the convolution Eq.(52). The amplitudeof the pressure wiggling at the contact surface is slightly larger than that of the cubic spline result.

26

SHU-ICHIRO INUTSUKA

Fifteen Degree Wedge Channel Flow

Y

1.0
0.5
0.0
-1.0

0.0

1.0

2.0

X
FIG. 9.
Results of calculation on \ fteen degree wedge channel ow" problem. Initially 3296 particles areowing inside the channel. Mach-number contours from 0.92 to 1.97 with an increment of 0.05 are plotted.

y

.5
0
-.5
-1

-.5

0

.5
x

1

1.5

2

FIG. 10.
Same as Fig.9 but positions of particles are plotted. The crosses denote the positions of \ghost
particles" that are introduced to mimic the rigid wall boundary condition.

27

SPH WITH RIEMANN SOLVER

Y

1.0

Fifteen Degree Wedge Channel Flow

0.5
0.0
-1.0

0.0

1.0

2.0

X
FIG. 11.
Results of calculation on \ fteen degree wedge channel ow" problem with the rst-order Riemann
Solver. Initially 3296 particles are owing inside the channel. The contour levels are the same as Fig.9.

Y

1.0

Fifteen Degree Wedge Channel Flow

0.5
0.0
-1.0

0.0

1.0

2.0

X
FIG. 12.
Results of calculation on \ fteen degree wedge channel ow" problem with the standard SPH. The
standard arti cial viscosity with = 1; = 2 are used. Initially 3296 particles are owing inside the channel. The
contour levels are the same as Fig.9.

Y

1.0

Fifteen Degree Wedge Channel Flow

0.5
0.0
-1.0

0.0

1.0

2.0

X
FIG. 13.
Results of calculation on \ fteen degree wedge channel ow" problem with the standard SPH. The
standard arti cial viscosity with = 2; = 4 are used. Initially 3296 particles are owing inside the channel. The
contour levels are the same as Fig.9.

28

SHU-ICHIRO INUTSUKA

Figure 9 shows our method's result. Mach-number contours from 0.92 to 1.97 with an incrementof 0.05 are plotted. Figure 10 shows particle positions. In this calculation, we realized the rigid-wallboundary condition by placing \ghost particles" in the wall. The in owing particles outside of theleft-hand side boundary were prepared at appropriate timesteps. The variable smoothing length isused ( = 1, Csmooth = 1:0). Initially, 3296 particles were owing inside the channel, which may becompared to the 3296 grids nite di erence calculation in [12]. Our results were satisfactory. For
comparison, Figure 11 shows the result for the rst-order Riemann Solver. The sharp features are
smeared out.
Figure 12 shows results of calculation with the standard SPH, where the standard arti cial viscositywith = 1; = 2 are used. Note that the contour levels are the same as Fig. 9. Owing to the moderate
strength of the shocks in this problem, the standard SPH could produce reasonable result except for
the region just after the Mach stem on the upper wall. A small hollow semicircle seen at x  0:5,
y  1 corresponds to the small value of Mach number (= 0:806) that is caused by an exceedingly
large value of the post shock pressure. The calculations with larger values of the arti cial viscosity
parameters tend to remedy this pressure overshooting. Figure 13 shows the result of the standard
SPH with = 2; = 4 where the minimum value of Mach number behind the Mach stem is 0.957.
In this case, however, the sharp features are smeared out. Even in this two dimensional problem ofmoderate Mach number, the present method shows its potential ability in analyzing supersonic ows.
It is clear that a grid-based method must be used in the actual study of this kind of steady rigid-wallboundary problem. This result is presented only to demonstrate the capability of this particle scheme.
The present method becomes more useful when we study the hydrodynamical problems without rigid
boundary that are common in astrophysics and space sciences.

5. SUMMARY

Smoothed Particle Hydrodynamics is reformulated by the formal convolution of the original hydrodynamics equations, and by a new action principle, in which the second order (in h) approximation
is used for the kinetic term of the Lagrangian function. The force acting on each particle is determined by solving the Riemann problem for each particle pair. The prescription for variable smoothinglength is also shown. These techniques are implemented in the strict conservation form. Numericalexamples involving an extremely strong shock are shown. The other test calculations will be described
elsewhere.
Although the method with a spatially constant smoothing length is formulated in a rigorous manner,
the method with the variable smoothing length is based on a crude approximation (Eqs. [79],[80]).
A more re ned technique for the variable smoothing length is to be studied. A better approximation
than the cubic spline interpolation in the numerical convolution is required to eliminate the \wiggle"
at the contact discontinuity in the blast wave problems (see Section 4.2).
This paper presented a piece of concepts that are not discussed in detail in the literature. Those are
the convolution of the original uid equation (Eqs. [20],[24]), the de nition and approximation for thevelocity eld (Eqs. [34], [39]), and the modi cation of the force due to dissipation (Eqs. [66], [67]).
Incorporating these concept, the nal evolution equation was cast into a form similar to the standard
SPH equation. Therefore those concepts may enable the rigorous examination of the accuracy andstability of the method and may enable further modi cation.
The numerical examples in Section 4 show that the present particle method based on Riemann
Solver can handle severe problems with strong shocks, those might include the description of explosion/implosion and supersonic jet phenomena. In this respect, further modi cation of the SPH
method in modeling relativistic ows is promising with the help of relativistic Riemann Solver [15].
In adittion, the Lagrangian particle methods have advantage over Eulerian grid-based methods, indescribing chemically reacting (multi) uid and radiatively heating/cooling uid [11]. This is because

SPH WITH RIEMANN SOLVER

29

we can simply assign the chemical composition and entropy to each particle as a uid element. This
direction has a wide area of applications.

ACKNOWLEDGMENT

The author thanks the anonymous referee for valuable comments. The author also thanks Toru Tsuribe, Yusuke
Imaeda, and Shoken M. Miyama for useful discussions.

APPENDIX
A.1. DERIVATION OF THE EQUATION OF MOTION

In this Appendix we show the derivation of the equation of motion from the Euler-Lagrange's
equation (Eq.[42]).

@L = ; X m Z @u W (x ; x ; h)dx
k
k
@ xi
@ xi
k
Z
;mi u @@x W (x ; xi ; h)dx:
i

(A.1)

The rst term on the RHS of this equation becomes the following:Z@u W (x ; x ; h)dx
k
@
xi
k
Z
X
= ; mk P2 @@
xi W (x ; xk ; h)dxkZXxi ; h) W (x ; x ; h)dx:= ; mk P2 mi @W (x@;
k
x
i
k

;

X

mk

Before we manipulate the second term on the RHS of Eq. (A.1), we note that
@W (x ; xi ; h) = ; @W (x ; xi ; h) :@ xi
@xWe also obtainZZ
f (x) @W (x@;x xi ; h) dx = ; @f@(xx) W (x ; xi ; h)dx;
by the integration by parts and W (x) ! 0 (jxj ! 1) .
Using the above identities, we can transform the second term on the RHS of Eq. (A.1) as

(A.2)

(A.3)
(A.4)

Z

;mi u @@x W (x ; xi ; h)dxiZ@
= +mi u @ x W (x ; xi ; h)dx
Z@u W (x ; x ; h)dxi@xZ
= ;mi P2 @@x W (x ; xi ; h)dx
ZX
= ;mi P2 mj @W (x@;x xj ; h) W (x ; xi ; h)dxjZXxj ; h) W (x ; x ; h)dx:= mi P2 mj @W (x@;
i
x
j
j
= ;mi

(A.5)

30

SHU-ICHIRO INUTSUKA

As a result, Euler-Lagrange's equation gives the following:Z

P @ W (x ; x ; h)W (x ; x ; h)dxij
2 @ xijZX+ mj P2 @@x W (x ; xi ; h)W (x ; xj ; h)dx:
jjx i = ;X

mj

REFERENCES(A.6)

1. W. Benz, Smooth Particle Hydrodynamics - a Review, in Proc Numerical Modelling of Nonlinear Stellar Pulsations:
Problems and Prospects, Les Arcs, France, March 1989, edited by J. Robert Buchler, (Kluwer, Dordrecht/Boston,1990), P.269.
2. P. Colella and P. Woodward, Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations, J. Computat.
Phys. 54, 174 (1984)
3. G. A. Dilts, Moving-Least-Squares-Particle Hydrodynamics {I. Consistency and Stability, Int. J. Numer. Meth.
Engng. 44, 1115 (1999)
4. G. A. Dilts, Moving-Least-Squares-Particle Hydrodynamics {II. Conservation and Boundaries, Int. J. Numer. Meth.
Engng. 48, 1503 (2000)
5. D. A. Fulk and D. W. Quinn, Hybrid Formulations of Smoothed Particle Hydrodynamics, Int. J. Impact Engng.
17, 329 (1995)
6. S. K. Godunov, Di erence Methods for the Numerical Calculation of the Equations of Fluid Dynamics, Mat. Sb.
47, 271 (1959)
7. R. A. Gingold and J. J. Monaghan, Smoothed Particle Hydrodynamics { Theory and Application to Non-SphericalStars, Mon. Not. R. Astron. Soc. 181, 375 (1977)
8. L. Hernquist & N. Katz, TREESPH - A Uni cation of SPH with the Hierarchical Tree Method, Astrophys. J. Supp.
70, 419 (1989)
9. S. Inutsuka, Godunov-type SPH, J. Italian, Astron. Soc. 65, 1027 (1994)
10. S. Inutsuka, A Numerical Scheme for Collisions of Stars: Godunov-type Particle Hydrodynamics, in Proc. IAU
Symp. 174, Dynamical Evolution of Star Clusters, ed. J. Makino & P. Hut, (Kluwer, Tokyo, 1995), p. 373
11. S. Inutsuka, Advanced Methods in Particle Hydrodynamics, in Proc. Numerical Astrophysics 1998, eds. S. M.
Miyama, K. Tomisaka, & T. Hanawa, (Kluwer, Boston, 1999) , p. 367.
12. D. W. Levy, K. G. Powell, and B. van Leer, Use of a Rotated Riemann Solver for the Two-Dimensional Euler
Equations, J. Computat. Phys. 106, 201 (1993)
13. C. C. Lin, Hydrodynamics of Helium II, Proc. Int. Sch. Phys. XXI, 93(1963)
14. L. Lucy, A Numerical Approach to the Testing of the Fission Hypothesis, Astron. J. 82, 1013 (1977)
15. J. M. Mart and E. Muller, Extension of the Piecewise Parabolic Method to One-Dimensional Relativistic Hydrodynamics, J. Computat. Phys. 123, 1 (1996)
16. J. J. Monaghan, Particle Methods for Hydrodynamics, Phys. Rep. 3(2), 71 (1985)
17. J. J. Monaghan, On the Problem of Penetration in Particle Methods, J. Computat. Phys., 82, 1 (1989)
18. J. J. Monaghan, Smoothed Particle Hydrodynamics, Annu. Rev. Astron. Astrophys. 30, 543 (1992)
19. R. A. Gingold and J. J. Monaghan, Kernel Estimate as a Basis for General Particle Hydrodynamics, J. Computat.
Phys. 46, 429 (1982)
20. R. Nelson and J. Papaloizou, Variable Smoothing Lengths and Energy Conservation in Smoothed Particle Hydrodynamics, Mon. Not. R. Astron. Soc. 270, 1 (1994)
21. J. P. Morris and J. J. Monaghan, A Switch to Reduce SPH Viscosity, J. Computat. Phys. 136, 41 (1997)
22. J. Serrin, Mathematical Principles of Classical Fluid Mechanics, Handbuch der Physik VIII-I, 125 (1959)
23. G. A. Sod, A Survey of Several Finite Di erence Methods for Systems of Nonlinear Hyperbolic Conservation Laws,
J. Computat. Phys. 27, 1 (1978)
24. B. van Leer, Towards the Ultimate Conservative Di erence Scheme. V. A Second-Order Sequel to Godunov's Method,
J. Computat. Phys. 32, 101 (1979)
25. J. Wisdom and M. M. Holman, Symplectic Maps for the N-Body Problem { Stability Analysis, Astron. J. 104, 2022
(1992)
26. P. Woodward and P. Colella, The Numerical Simulation of Two-Dimensional Fluid Flow with Strong Shocks, J.
Computat. Phys. 54, 115 (1984)