JOURNAL OF COMPUTATIONAL PHYSICS 52, 374-389 (1983)

Shock Simulation by the Particle Method SPH

J.J. MONAGHAN

Mathematics Department,
Monash University, Clayton, Victoria, 3168, Australia

AND

R. A. GINGOLD

Mt. Stromlo and Siding Spring Obs.

Woden, Private Bag. A.C.T., Australia
Received September 2, 1982; revised March 15, 1983

The particle method SPH is applied to one-dimensional shock tube problems by iricor-
porating an artificial viscosity into the equations of motion. When the artificial viscosity is
either a bulk viscosity or the Von Neumann-Richtmyer viscosity, in a form analogous to that
for finite differences, the results show either excessive oscillation or excessive smearing of the
shock front. The reason for the excessive particle oscillation is that, in the standard form, the
artificial viscosity cannot dampen irregular motion on the scale of the particle separation
since that scale is usually less than the resolution of the interpolating kernel. We propose a
new form of artificial viscosity which eliminates this problem. The resulting shock simulation
has negligible oscillation and satisfactorily sharp discontinuities. Results with a gaussian inter-
polating kernel (with second-order errors) are shown to be greatly inferior to those with a
super gaussian kernel (with fourth-order errors).

1. INTRODUCTION

When applying particle methods to the simulation of shocks it is convenient to use

the artificial viscosities proposed and tested for finite difference methods (Roache,
1975). These artificial viscosities can be considered as forms of viscous pressure, and
they can be easily included in the equations of motion.
In order to calculate the viscous pressure it is first necessary to calculate the
divergence of the velocity, V . v. For particle methods which use a grid, V . v can be
calculated by finite differences. For a particle method such as SPH (smoothed
particle hydrodynamics) (Gingold and Monaghan, 1977, 1982; Monaghan, 1982)
V . v can be calculated using an interpolating kernel with the particles as inter-
polating points. The length scales V . v are therefore bounded below either by the
resolution of the grid or the resolution of the interpolating kernel. However, in order
to maintain reasonable accuracy, the resolution length scale must be much greater
than the typical separation of the particles, and irregular motion on this latter scale is
374
0021.9991/83 $3.00
Copyright 0 1983 by Academic Press. Inc.
All rights of reproduction in any form reserved.
 THEPARTICLEMETHOD SPH 375

then only weakly affected by the artificial viscosity. The scale separation creates
problems in shock transitions for two reasons. The first is that in the highdensity
side of the shock the scale separation is exacerbated. The second is that, in the shock
transition, the particles mimic molecules and, in the absence of an appropriate
artificial viscosity, the shock is simulated by inducing irregular particle motions on
length scales corresponding to the particle separation. Artificial viscosities calculated
in the usual way are very inefficient at damping this irregular motion. They can be
made to do so by increasing the artificial viscosity, but the result is excessive
broadening of the shock front. Experiments confirming these ideas are described in
Section 3 of this paper.
It is clear that a particle method requires a damping term which acts more directly
on the relative motion of particles. In Section 4 we describe a new artificial viscosity
which is equivalent to a bulk viscosity but which acts effectively on irregular motion
on the short scale. Experiments described in Section 5 confirm that this artificial
viscosity produces negligible oscillation and good resolution of the shock front and
contact discontinuity.
Recent theory and experiment with finite difference methods (Van Leer, 1979)
show that methods which have third-order accuracy, when used correctly, give
excellent results for shock transitions. In this paper we compare results using a
gaussian kernel, which has first-order accuracy, with those using a super gaussian
kernel which interpolates with third-order accuracy. The latter gives excellent results
which rival those found by Sod (1978) for the best of a variety of finite difference
methods.
In the following section we establish the equations of motion from the exact
equations using the ideas and formalism of Monaghan (1982).

2. EQUATIONS OF MOTION

The SPH equations, like those for other particle methods, were first derived by
intuitive methods (Lucy, 1977; Gingold and Monaghan, 1977). Recently (Monaghan,
1982) it has been shown how the SPH equations can be derived from the exact
equations of motion. The formalism used for this purpose can also be applied to
spectral and finite difference methods which are then seen to differ from particle
methods by (a) using a different interpolating kernel, and (b) not using (in general)
interpolating points which move with the fluid. The main points of the analysis are
the following:
(i) All interpolation methods involving linear operations on the function inter-
polated (e.g., spectral or local polynomial interpolation) can be written in the form

(‘4(r))=(DA(TI) W(r,r’,h)dr’, (2-l)

376 MONAGHAN AND GINGOLD

where A(r) is the function interpolated, D is the domain, and W is a kernel with the
following properties:

(4 ID
W(r, r’, h) dr’ + 1 as h -+ 0. (2.2)

(b) If ,4(r) is a continuous function

Mr)WW as h -+ 0. (2.3)

The kernel therefore mimics a delta function and it does so more closely as h + 0.
(ii) If the interpolation points are distributed in space with number density
n(r) then (A) can be approximated by

x.i -iA9 Wr, rjy h), (2.4)

where, for any function B, Bj = B(rj). If the domain is one dimensional, and the
points are equi-separated,the choice of an appropriate kernel allows all the standard
interpolation formula to be recovered from (2.4). The expression (2.4) for (A) is the
general SPH interpolation formula.
(iii) Equations for numerical work can be constructed by multiplying each
term of the exact equations by the kernel and integrating over the domain where a
solution is required. Integration by parts with, if necessary, approximations for
nonlinear terms gives the equations for numerical work.
The Momentum Equation
We first consider inviscid flow for which the momentum equation is
dv
= f Vp. (2.5)
z-
Following (iii) we find

J,$ W(r,r’,h)dr’=-jDiVpW(r,r’,h)dr’. (2.6)

It is sufftcientiy general for most purposes to assumethat

W(r,r’,h)E W(ir-r’I,h). (2.7)

A consequenceof this assumption is that for any function $

(V$> = V(4) + j $(r’) W(lr - r’ I, h) n da, (2.8)

S
 THE PARTICLE METHOD SPH 377

where the surface integral is over the surface S of the domain D. For the problems we
consider all the variables vanish on the surface and we can take

W) = V@)* (2.9)
In general, however, surface terms appear in the equations of motion.
Since the RHS of (2.6) involves the nonlinear combination Vp/p we need further
approximations to evaluate the integral. An obvious approximation is to choose

= -&v(P), (2.10)

but if exact conservation of linear and angular momentum is required it is better to

note that

VP
-=v $ +$vp, (2.11)
P 0

so that (2.6) can be approximated by

(f)=--q-)-$+(p). (2.12)

No special interpolation method has yet been specified. We now specialize to particle
methods. We assign a mass m to each interpolation point and call it a particle, and
we move the particle with the acceleration an element of fluid would experience at the
position of the particle.
If we assign a mass to each particle the density p is given by

mn(r). (2.13)

Referring to (2.4) we find (dropping angle brackets for convenience) that (2.12)
becomes, at particle i,

(2.14)

where Wij = W(jr, - rjj, h) and Vi denotes the gradient taken with respect to the
coordinate ri.
An artificial viscosity in the form of a viscous pressure q can be included in the
equation of motion by replacing p in (2.14) by p + q. The new system of equations
conserves total linear and angular momentum exactly.
It can be shown (Monaghan, 1982) that if the same procedure is applied to the
378 MONAGHAN AND GINGOLD

continuity equation it leads to an equation for the time change of the interpolated
density

@(r))=m 5 W(lr-rjITh) (2.15)

j=I

which is satisfied automatically locally (that is, at each particle). Total mass is, of
course, conserved exactly.
The Energy Equation
The energy equation for particle methods can be derived directly from the exact
energy equation using the procedure described earlier (Monaghan, 1982). However,
in Section 4 we shall introduce an artificial viscosity which is not easily related to a
viscous term in the exact equations. It is therefore preferable to construct an energy
equation directly from (2.14) and use our knowledge of the exact energy equation to
guide the identification of the terms in the energy equation. From (2.14) we find

mCVi .~=-m2CCuijvi-ViW,,
i 1 i

uij ..- Pi+4i I Pj + Sj

(2.17)
Pt pi”

Interchanging the labels in the double summation we can write (2.16) in the form

(2.18)

With the assumptions we have made, W, = Wji, uij = u,~ and Vi W, = -Vj Wji.
Accordingly (2.18) can be written in the form

~mvi.~=--$$70ijvij.ViWij, (2.19)
J

where vii = vi - vj. The LHS of (2.19) is the rate of change of the total kinetic
energy. We provisionally identify

FT uijvij . vi w,

as the rate of change of thermal energy per unit mass at the position of particle i. If
our choice is correct (2.20) should be the particle representation of

-%.v,
P
 THE PARTICLE METHOD SPH 379

when q = 0. The correspondence can be established by noting that

-$o.“=-t[~~v.~“)-“.vp}+v. (T) -“.V(f)]. (2.21)

This somewhat surprising combination of terms is suggested by the relation (2.11)

which leads to (2.14). If the RHS is converted into the particle form using approx-
imations equivalent to those used in deriving (2.14) we recover (2.20).
If ui is the thermal energy per unit mass the thermal energy equation is

(2.22)

where the second term on the RHS of (2.22) is the viscous dissipation.
We could use the energy equation in this form since it would lead to exact energy
conservation. It is, however, preferable for our purposes to replace the first term on
the RHS of (2.22) by its exact equivalent so that (2.22) can be rewritten so as to give
the rate of change of entropy per unit mass s:

(2.23)

The advantage of this form is that a relatively crude algorithm for the integration of
(2.23) will give satisfactory results over the bulk of the flow. We have, however, used
both (2.22) and (2.23) in our calculations and the change in the results is negligible.

The Artificial Viscosity

For the present we consider two viscous pressures. These are (a) the Von
Neumann-Richtmyer viscous pressure defined by

q = apP(V * zg2, v-v<o,

(2.24)
= 0, v * v > 0,

where a is a constant, and (b) a bulk viscosity for which

q = -aphcV . v, v*v<o,
(2.25)
= 0, v*v>o,

where a is a constant and c is the speed of sound.

Both viscous pressures require the calculation of V . v. In the particle represen-
tation we can use either

vj * vi wij
mC (2.26)
j Pj ’
380 MONAGHANANDGINGOLD

or
!I.\- vji . vi wij5 (2.27)
Pi 7
where the latter has been derived from the relation

v *v = f [V * (pv) -v . VP]. (2.28)

The expression (2.27) is more accurate than (2.26); in particular it gives the correct
result when v is constant. In the calculations to be described in Section 3 we use
(2.27), but the results are very similar to those found when (2.26) is used.
As we noted in Section 1 the estimates of V . v are only weakly affected by
irregular variations of v on a scale & since the other terms in the summations in
(2.26) or (2.27) vary slowly on a scale -& and only the average of the irregular
variations in v contributes to V . v. This average is small.

3. NUMERICAL RESULTS-STANDARD ARTIFICIAL VISCOSITY

We consider a shock tube problem for a perfect gas analogous to that considered
by Sod (1978). The initial conditions are

x < 0; p= 1, p=l
(3.1)
x > 0; p = 0.25, p = 0.2154.

In the absenceof dissipation the energy equation becomes

p=A(s)pY, (3.2)
with y = 1.4, and A(s) = 1 for those particles with x < 0 initially, and A(s) = 1.25 for
those particles with x > 0 initially.
For the particle simulation we distribute N’ particles uniformly in 0 to 0.6 and 4N’
uniformly in -0.6 to 0.0. Typically we take N’ = 80 and h twice the particle
separation in the low-density region. No special boundary conditions are applied at
the end points since, for the time elapse we consider, effects at the end points do not
have time to propagate to the shock front and contact discontinuity which are near
x = 0.
The time step dt is limited by the Courant condition. For the calculations to be
discussedhere 6t = Min(0.3h/ci), where ci is the, speedof sound ar particle i.
The momentum equation was integrated using the leapfrog algorithm. The thermal
energy equation (2.23) was integrated by writing it in the form

d-i
-= AQ(r-1)~
dt P ’
 THE PARTICLE METHOD SPH 381

where Q is the RHS of (2.23), and then updating A with the algorithm

,‘I(“+‘) d’“‘(1 + &(y - l)p’“‘Q/p’“‘). (3.4)

This algorithm is only first-order accurate in the time but, as shown by the results in
Section 5, it gives satisfactory accuracy once the correct artificial viscosity is used.
The interpolated pressure is calculated from (3.2) by first calculating @A) then
setting

(P> = (&x,PY l. (3.5)

For these calculations we use two interpolating kernels. The first is the gaussian

W(u,h)=-L -G/h=

hfi ’

which gives an interpolation error of O(h2). The second is the super gaussian

W(u, h) = -h;-n e-“lh2 (f+‘h’),

which gives interpolation errors of O(h4). For a given h the super gaussian kernel has
a better resolution than the gaussian kernel (3.6). When the kernel, or its derivative,
is required in the computation it is obtained by interpolation from an array. For this
reason the CPU time is independent of the kernel.
The expressions for the density, acceleration and viscous dissipation formally
involve a summation over all the particles. In practice this summation need only
include particles out to 3h. A substantial saving in time is achieved when N’ is large
(2150) by using cells as a bookkeeping device. If the particles are assignedto cells of
width 3h, and identified through linked lists, the calculation time is proportional to N
and not N2 as it would be if the complete summation was performed. The procedure
for carrying this out is described in detail by Hackney and Eastwood (1981) in their
discussion of short-range forces in particle simulation models. The resulting algorithm
is nearly as quick as PIC, but has the advantage that the spatial resolution is superior
since it is approximately the particle separation not the cell width as in PIC.
The results for a typical calculation with the Von Neumann-Richtmyer and bulk
viscosities are compared with the exact solution in Fig. 1. The calculations use the
super’ gaussian kernel (3.7) with h = 0.015. The numerical results show good
agreementwith the exact solution for p and p, but the shock front is broadened over
-4h. The velocity shows the irregular oscillations discussed earlier. These
oscillations, being irregular and of short wavelength, have little effect on the p and p
profiles. The oscillations are particularly large in the case of the bulk viscosity (2.25)
where the amplitude of the oscillation reaches 30 percent of the correct velocity. By
increasing GIby a factor 4 the post shock oscillations are reduced to -10 percent, but
 1.2

1.0

0.6

0.6 : 06
::
0.. L 0.1
0.2 0.2

0.0 -I I 1 . 7 . 0.04 , I 7 I I I 1 I
-0.. -0.2 0.3 0.2 -0.1 -0.2 0.0 0.2 0.
x x

“,I
1.0

0.6 :, 0.6
E
0.6

P
0.. -v
':
I.
0.2 -e--m L
+ . I
-0 a -0 2 00 02 c
x

FIG. I. Pressure, density, and velocity profiles for the shock tube problem described in Section 3 with N’ = 80 and h = 0.015. The upper frames
are for the Von Neumann-Richtmyer viscosity. The lower frames are for the bulk viscosity. The exact results are shown: ---. The SPH results arc
shown by dots and full lines. The calculation uses the super gaussian kernel.
 THE PARTICLE METHOD SPH 383

the shock is broadened over 10h. More regular oscillations appear when the Von
Neumann-Richtmyer viscosity is used. In this case a large fraction of the energy in
the oscillations is associated with a length scale kh and, as a consequence, the
oscillations show up in the p and p profiles. As before the oscillations can be removed
by increasing a sufficiently, but the broadening of the shock front becomesunaccep-
tably large.
For a given a the oscillation increase in amplitude as the magnitude of the initial
discontinuity is increased. In the next section we discuss the construction of an alter-
native artificial viscosity which removes the oscillation while retaining good
resolution of the discontinuities.

4. THE NEW ARTIFICIAL VISCOSITY

The results described in Section 3 show the need for an artificial viscosity which
acts more directly on the relative motion of the particles. The following intuitive
arguments are those that led us to an appropriate form of artificial viscosity.
We need to replace the factor

where IZ, is the new term which acts as an artificial viscosity. To keep the analysis
simple we confine our attention to motion in one dimension and look for a flij that is
similar to a bulk viscosity. We expect n, to be an approximation to

hc au
-a--,
P ax

which is sensitive to the relative motion of two particles on a scale ch. Since
vi-vj 1
)I
%i-
hi &Ii
--=- Vi- Vi+Xjiaxi+“’
xij xij xij [ i -z(’

we could try

n..=-nhC’!k
IJ
Pi xij’
However, (4.3) has two deficiencies: (a) the xij in the denominator might lead to a
384 MONAGHAN AND GINGOLD

viscous force which is too large and (b) II, # ZZ,, so that linear momentum is not
automatically conserved. To correct the first deficiency we replace

-
1 xij
by x; + ch2 ’
xij

where 0 < E< 1, and in practice we find E2 0.1 is adequate for the numerical
experiments we describe here. When the density contrast is larger a smaller E is
appropriate. The general rule is that if the initial particle separation on the high-
density side is f times that on the low-density side then E-f 2.
To correct the second deficiency we replace ci and pi by E, and pij, respectively,
where, for any function B, we define
gij := ;(B, + Bj). (4.5)
The final expression for n, is

n,, = -a hE, ‘ijxij

1, vijxij < 0,
pij x;++h2 ’
= 0, UijXij > O, (4.6)
where uijxij is used to determine whether &/ax is positive or negative. The artificial
viscosity is only neededin the latter case. In particle terms the viscosity only operates
on approaching particles. If the kernel is the gaussian (3.6) the total viscous
acceleration is

We can write (4.7) in a form which is simple to interpret by defining

cjjx; wij
tTij := (4.8)
~ij(X~ + Eh’) ’
and
ai :=x fJij. (4.9)

With the use of cij and ci we can write (4.7) in the form

2t?lUOi
h (Vi - V;.), (4. IO)
where
5, = Ci Oijvj
(4.11)
’ cjoij .
 THE PARTICLE METHOD SPH 385

Expression (4.10) shows that the viscous force attempts to drive ui to zYi,where the
latter is an average velocity taken with weight uij. Expression (4.10) bears a formal
similarity to the viscosity used by Marder (1975), who, however, uses a different
definition of Gi and a different rule for specifying the coefficient multiplying (ui - Vi).
Furthermore, since Marder uses a grid with approximately 10 particles per cell, the
resolution in his calculations is -5 times coarser than can be achieved by the SPH
calculations.
In the appendix we show that if h is sufficiently small, and summations are
replaced by integrations, the viscous acceleration is
1 ha av
yap% (PC& 1y
which establishes, as expected, that the artificial viscosity is approximately equivalent
to a bulk viscosity in the exact equations.
The energy equation we use is (2.23) with ZI, instead of (qi/pf + qj/pj).

1.2 1.2 ,

1.0
1
1.0 \ 1
0.6 0.8
i
K
: 0.6 / \
0.6
i
0.4 E -L--T,>
0.4 .I
I

0.2

0.02 0.0-1

0.8

FIG. 2. As for Fig. 1 except that the SPH calculation uses the new artificial viscosity. The kernel is
super gaussian.
386 MONAGHANANDGINGOLD

5. NUMERICAL RESULTS-NEW VISCOSITY

The numerical results for the shock tube problem described in Section 3 using the
SPH method with I7, are illustrated in Fig. 2. For this calculation the super gaussian
kernel (3.7) was used. The results are much better than those found using the
standard artificial viscosity. The shock front and contact discontinuity are now
broadened over only 2/r, which improves on the results discussed in Section 3 and
rivals the best of the methods discussed by Sod (1978). By comparison with the
results in Section 3 oscillations are negligible. The end points of the rarefaction wave
are not perfectly sharp but the error is small.
The pressure profile shows a weak blip at the contact discontinuity. It is due to the
fact that the variable A is‘discontinuous at the contact discontinuity, and when (Ap)
is formed it is inevitable that (p) is slightly lower to the left of the contact discon-
tinuity, and slightly higher to the right. This blip increases with increase in the jump
in A at the contact discontinuity, but it has a negligible effect on the motion.
In recent experiments we have found that if the pressure acceleration term Vp/p is
written in the form

Py-21wP)+ (Y- l)AVPl

‘.21 ‘.*1
1.0 -1
\

0.8
*

0.6

“\, \ ----T-,
0.4
‘: I
1
0.2 ;...
-a-___-,

i
0.04 , I . r I I I I 0.0 4 I I I
-0.4 -0.2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 c14
x x

0.8 1.4 _

FIG. 3. As for Fig. 2 except that a normal gaussian kernel was used.
 THEPARTICLEMETHOD SPH 387

the resulting equations of motion give results comparable in accuracy to those

described here without the blip in the pressure curve.
In Fig. 3 we show the results for the same problem using the gaussian kernel (3.6).
The shock front and contact discontinuity are now broadened over 3h, and the end
points of the rarefaction wave are more rounded. If the same kernel is used with the
standard viscosities the results are inferior to those shown in Fig. 3 and they show the
excessiveoscillation in the velocity already noted in Section 3. These results might be
expected becausethe standard gaussian kernel interpolates with lower accuracy than
the super gaussian kernel. However, a priori estimates of accuracy are difficult
becausein particle methods there are two sources of interpolation error. One of these
is the truncation error known to be greater for the standard gaussian kernel. The
other source of error is the estimate of integrals by sums, as in passing from (2.1) to
(2.4). We are uable to estimate how this latter error is affected by changing the
kernel. It is clear, however, that in practice the super gaussian kernel is superior.

6. GENERALIZATION

The artificial viscosity term can be easily generalized to an arbitrary number of

dimensions by writing it in the form
ahei, vii . rij
nij=- vij * rij < 0,
pij r; + ch2 ’
= 0, vij * rij > 0.

1
z

-1

-2

-3
0.0 0.2 0,4 0.6 0.6 1.0 1.2 1.4 00 o-2 04 0.6 06 10
TIME TIME

FIG. 4. The variation of the z coordinate with time for some selected particles in the three-
dimensional collapse of an isothermal self-gravitating cloud rotating about the z axis. The frame on the
left shows the collapse with no artificial viscosity and illustrates vigorous streaming through the
equatorial plane. The frame on the right shows the collapse when the new artificial viscosity is used in
the z component of the momentum equation. The streaming is now negligible.
388 MONAGHAN AND GINGOLD

In some problems (for example, the collapse of rotating, self-gravitating clouds)

where the correct transport of angular momentum is important, it is preferable to
include flii only in the component of the momentum equation which is parallel to the
angular momentum vector (Gingold and Monaghan, 1982), and replace vij . rij by
“iJ . rU, where u, is the velocity component of particle i parallel to the angular
momentum vector. In that context l7, has the further advantage that it acts very
effectively to prevent unwanted particle streaming. As an example we show in Fig. 4
the z coordinates for particles in a collapsing, rotating cloud both with and without
the artificial viscosity.
A further extension of the ideas described in Section 4 is to construct a ZZ, for a
viscosity equivalent to the Von Neumann-Richtmyer viscosity. In one dimension the
appropriate nu is
ah’v,f
pi,(x; + ch’) ’

where a is a constant.
Experiments with this n, produce results very similar to those described in
Section 5.

7. CONCLUSIONS

Our numerical experiments show that it is possible to choose an artificial viscosity

for the particle method SPH which results in negligible post shock oscillations. If a
kernel which interpolates with third-order accuracy is used, the discontinuities are
resolved with an overall accuracy that rivals the best of the methods considered by
Sod (1978), but it is not quite as accurate as Van Leer’s (1979) method. The present
method has the advantage that it is very easy to apply, and its generalization and
application to more than one dimension are straightforward.

APPENDIX

The viscous acceleration in one dimension is

a, = mah c f&i ‘+ij vi w..1J ’

J PiJ xj?J+Eh’
(A.11
where Vi is now equivalent to a/ax,. Since the presenceof Wij limits the contribution
to the summation to those points with xi-xi we expand the variables about this
point. Thus
2
v,=xuv; --v; xii + a**,
2 64.2)
 THE PARTICLE METHOD SPH 389

fj..
‘.pi $E c ‘+..., (A.31
Pij Pi 2 0Pi

where, for any function B(x), Bf E B(x,)/ax,. If (A.2) and (A.3) are substituted into
(A. 1) and the following results are used,

(i> m 5’ xpi w, + 7 vi wij + (p)[, (A.41

7 x:+&h2
and

(ii) m Y x;vi wij +mmxijViWij

7 x&+&h’ j

=ximTViWij-Vi (mCxjWij)

+ xi@); - (xp)i’

+ @)i* (A.5)

(A.l) becomes (dropping angle brackets for convenience)

The error depends on how well (A) represents A and how well (2.4) approximates
(2.1). The first source of error is a truncation error (Monaghan, 1982) which for the
kernel (3.6) is of order h* and for the kernel (3.7) is of order h4.

REFERENCES

R. A. GINGOLD AND J. J. MONAGHAN, Mon. Not. Roy. Astron. Sot. 181 (1977), 375-389.
R. A. GINGOLD AND J. J. MONAGHAN, J. Comput. Phys. 46 (1982), 429-453.
R. W. HOCKNEY AND J. W. EASTWOOD, “Computer Simulation Using Particles,” McGraw-Hill, New
York, 1981.
B. M. MARDER, M&h. Comp. 29 (1975), 434.
J. J. MONAGHAN, SIAM. J. Sci. Statist. Comput. 3 (1982), 422.
P. J. ROACHE, “Computational Fluid Dynamics,” Hermosa Pub., Albuquerque, 1975.
G. A. SOD, J. Comput. Phys. 27 (1978), 1-31.
B. VAN LEER, J. Comput. Phys. 32 (1979), 101-136.